The Osteogenic Function of Danggui Buxue Tang, a Herbal Decoction Containing Astragali Radix and Angelicae Sinensis Radix, Is Optimized by Boiling the Two Herbs Together: Signaling Analyses Revealed by Systems Biology

: The therapeutic efﬁcacy of a herbal mixture, being multi-target, multi-function and multi-pathway, is the niche of traditional Chinese medicine (TCM). Systems biology can dissect the network of signaling mechanisms in a complex biological system. In preparing TCM decoctions, the boiling of herbs together in water is a common practice; however, the rationale of this speciﬁc preparation has not been fully revealed. An approach of mass-spectrometry-based multi-omics was employed to examine the proﬁles of the cellular pathways, so as to understand the pharmacological efﬁcacy of Danggui Buxue Tang (DBT), a Chinese herbal mixture containing Astragali Radix and Angelicae Sinensis Radix, in cultured rat osteoblasts and mesenchymal stem cells. The results, generated from omics analyses, were compared from DBT-treated osteoblasts to those of treating the herbal extract by simple mixing of extracts from Astragali Radix and Angelicae Sinensis Radix, i.e., herbal mixture without boiling together. The signaling pathways responsible for energy metabolism and amino acid metabolism showed distinct activation, as triggered by DBT, in contrast to simple mixing of two herbal extracts. The result supports that boiling the herbs together is designed to maximize the osteoblastic function of DBT, such as in energy and lipid metabolism. This harmony of TCM formulation, by having interactive functions of two herbs during preparation, is being illustrated. The systems biology approach provides new and essential insights into the synergy of herbal preparation. Well-deﬁned multiple targets and multiple pathways in different levels of omics are the key to modernizing TCM.


Introduction
The preparative procedure in making a herbal extract in traditional Chinese medicine (TCM) is crucial for its pharmacological efficacy, which determines the final chemical composition within the extract. For an example, the Chinese Pharmacopeia has described four different herbal decoctions that could be derived from Coptis chinensis roots. These decoctions have the same herbal materials; however, they show distinct pharmacological activities and modes of action. The distinction is known to be triggered by the variation in preparation methods [1]. In general, TCM herbal formulae are typically composed of two

Preparation and Quality Control of Herbal Medicines
The preparation of DBT was described in detail previously [4,6,7,10]. Roots of threeyear-old A. membranaceus var. mongholicus (AR) from Shanxi Province and two-year-old A. sinensis roots (ASR) from Minxian of Gansu Province were utilized. Herbal medicine was examined by Dr. Tina Dong (herbalist) based on Standards of Materia Medica in Hong Kong. To extract herbal materials, DBT (AR and ASR in a weight ratio of 5:1) was extracted in eight volumes of water (v/w) at boiling for two hours: this process was performed twice. The mixture was dehydrated under vacuum. Identification and quantification of active compounds of DBT was conducted on an Agilent high-performance liquid chromatography (HPLC) system (Agilent, Waldbronn, Germany), coupled with a diode array detector (DAD) and an evaporative light scattering detector (ELSD). The herbal extracts were isolated by an Agilent C 18

Rat Osteoblast and Micromass Culture
Animal protocols were revised and approved by the Animal Experimentation Ethics Committee of the Department of Health, Hong Kong (No. 17-283 for Animal Ethics Approval), under the instructions of "Principles of Laboratory Animal Care" (NIH publication No. DH/HA&P/8/2/3). The postnatal day 1 SD rat was dissected to obtain calvarias. Tissues were digested by 1% trypsin for 10 min, and 0.2% collagenase for 65 min, respectively [11]. Afterwards, the supernatant was obtained via centrifugation at 1500 rpm for 5 min. Osteoblasts were incubated in MEM-α, supplemented with 10% FBS and 1% penicillin/streptomycin. Proliferation of osteoblasts was conducted by 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetra-zolium bromide (MTT) assay. For micromass culture, rat embryos were utilized to cultivate mesenchymal progenitor cells from limb buds. The ectoderms were removed after enzymatic digestion of limb buds with 0.1% trypsin and 2.4 U of dispase ll for 20 min. The cell number was adjusted to 25 × 10 6 cells/mL. The micromass was maintained with the CMRL 1066 medium, 10% FBS and 1% penicillin/streptomycin. The micromass was fixed with 4% paraformaldehyde for 15 min, incubated in Alcian blue 8 GS solution (0.1 mg/mL) for 1 h, then washed and immersed with glycerol for Alizarin Red S (40 mM, pH 4.2) staining.

Alkaline Phosphatase (ALP) and Collagen Quantification
ALP was extracted in lysis buffer (10 mM HEPES, pH 7.5 and 0.2% Triton X-100). ALP enzymatic reaction was examined by mixing the cell lysate with 10 mM p-nitrophenyl phosphate (as a substrate in enzymatic reaction) in a buffer (pH 10.4) containing 0.1 M glycine, 1 mM MgCl 2 and 1 mM ZnCl 2 at 37 • C, and colorimetric analysis was conducted at 405 nm [11]. Dexamethasone and vitamin C were used as a positive control [12]. For determination of collagen, osteoblasts were cultured for 72 h with drug treatment. The cells were incubated in methanol overnight and placed in 0.1% picrosirius red staining solution (100 µL/well) for 3 h. The staining was recognized by a microscope.

Adhesion Assay
Osteoblasts were cultured for 72 h with drug treatment. Cells were trypsinized for 1-2 min and cultured in fresh cultured medium. The cell suspension was then counted to a cell concentration of 4 × 10 5 cells/mL, and was placed onto a cell culture plate. The cells were cultured for 30 min for attachment, and the suspended cells were collected. The unattached cells were counted via phase-contrast microscopy.

Shotgun Proteomics
For protein extraction, osteoblasts were subjected to freeze-thaw cycles repeated three times, and sonication in 8M urea buffer (0.1% SDS, pH = 7.4) for 5 min. Next, the cell lysate was precipitated in cold acetone. The proteins were resuspended in 4 M urea (pH = 6.5). Forty milligrams of protein was reduced by DTT and alkylated by iodoacetamide (IAA). The modified proteins were then cleaved by trypsin (1: 50 w/w) for 18 h at 37 • C. Next, the peptide sample was desalted by C 18 ZipTip (Millipore, Darmstadt, Germany) and dried by vacuum. Peptides were dissolved in 0.1% formic acid (FA) and were directly loaded onto a C 18 capillary column (75 µm × 25 cm; 2 µm, 100 Å). The solvent elution was optimized using an Ultimate nanoLC system (Thermo Fisher Scientific, Waltham, MA) at a flow speed of 300 nL/min, and a 120 min LC gradient of 2-90% acetonitrile (ACN) in 0.1% FA was utilized to isolate peptides. The eluted peptides were detected by an Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher Scientific). The electrospray ionization voltage was set at +2.3 kV, and the ion transfer tube temperature was set at 300 • C. MS machine setting was as follows: 1 microscan for MS 1 scan at 60 K resolution, MS 2 at 30 K resolution; full MS mass range: m/z 400-1500; and MS/MS mass range: m/z 100-2000. Automatic gain control targeting for MS 2 was 40 K; maximum injection time was 20 ms; higher-energy C-trap dissociation energy was 35% and dynamic exclusion duration was 4 s.

Protein Identification and Quantification
MS data were analyzed via Thermo Scientific™ Proteome Discoverer™ 2.2. MS 2 data were searched with SEQUEST ® HT against a database of Rattus norvegicus. Carbamidomethylation (+57.021 Da) of cysteine residues was set as a fixed modification, and oxidation of methionine residues (+15.9949 Da) and acetylation of the protein N-terminus (+42.0106) were considered as variable modifications. MS 1 tolerance was set as 20 ppm, and MS 2 was set at 0.8 Da tolerance. Peptide spectral matches were validated using a percolator algorithm, based on q-values at a 1% false discovery rate. The proteomics data collection was conducted using an Rt-Aligner and feature mapper nodes, created for untargeted label-free quantification workflow in Proteome Discoverer.

Lipidomics Analysis
The cellular lipids were obtained by homogenizing 1 × 10 7 osteoblasts with 300 µL LC-MS water, 600 µL methanol and 450 µL chloroform. The mixture was subjected to water-organic two-layer separation by centrifugation at 12,000 g for 15 min. The lower layer was obtained and dried under nitrogen. Aliquots of 60 µL lipid samples were pooled as a QC sample. The cellular lipids were measured on a Waters Acquity UPLC connected to a high-resolution mass spectrometer (TripleTOF 4600, AB SCIEX). The cellular lipids were analyzed on a Waters UPLC Acquity instrument, coupled with a high-resolution MS (TripleTOF 4600, AB SCIEX). A Waters C 18 column (1.7 µm, 2.1× 100 mm) was utilized for the separation. The final acquisition method was illustrated as below. The mobile phase for lipidomics was A = ACN/water (60:40) with 10 mM NH 4  13.1-20 min: 0% B. The solvent speed was optimized at 0.2 mL/min. The MS conditions were set as follows: ion source gas 1 at 45, ion source gas 2 at 45, curtain gas at 30, temperature at 450 • C, ion-spray voltage floating at ± 5 kV, de-clustering potential (DP) at 100 V and collision energy (CE) at 10 V. For TOF-MS scan, the accumulation time was set as 0.1 s per spectra, and TOF masses were acquired from 200 to 1200 Da. For product ion experiment, the accumulation time was at 0.05 s per spectra, and masses were acquired from 100 to 1160 Da, DP was set at 100 V and CE was at 35 V with ± 15 V collision energy spread. Information-dependent acquisition (IDA) was selected to perform MS/MS scan, and its parameters were set as follows: exclude isotopes within 4 Da, mass tolerance at 10 ppm and the maximum number of candidate ions per cycle at 20. The "for" and "after" were chosen in the "exclude former target ions" part and were set at "for 15 s after two occurrences" (the time was usually set at half-length of a signal). In addition, "dynamic background subtracts" was chosen in the IDA advance module.

Data Extraction and Processing
Chromatographic peak identification and alignment were performed by using Progenesis QI 2.3 (Nonlinear Dynamics, Newcastle upon Tyne, UK). Those unstable metabolites were filtered out by applying a cut-off on the coefficient of variation at >30% in QC samples. Matrix of normalized ion abundance was exported to SIMCA ® (version 13, Umetrics AB) for multivariate data analysis. The potential candidates were selected from S-plots of OPLS-DA. The biomarkers were identified with mass fragmentation and matched with Human Metabolome Databases, Kyoto Encyclopedia of Genes and Genomes, METLIN, LipidMap (www.lipidmaps.org, accessed on 16 August 2020) as well as LipidBlast based on their mass fragmentation, retention time and mass accuracy.

Pathway and Statistical Analysis
Signaling analysis was conducted by KOBAS (KEGG Orthology-Based Annotation System), a software interpreting sequences with KEGG orthology terms and identifing the enriched pathways in the queried sequences, as compared to the background. The pathway datebsese employed here were KEGG, Reactome, Panther, and GO analysis. Pathway analysis was conducted using Ingenuity Pathway Analysis (Qiagen). The multivariate analysis was conducted using SIMCA ® . Data were represented as mean ± standard error of mean (SEM). Statistical analysis was performed with student t-test and Dunnett's test (SPSS, version 13). Statistically, difference was classified as significant: * p < 0.05, ** p < 0.01 and *** p < 0.001.

The Osteoblastic Function of DBT Requires Boiling Herbs Together
The herbal extracts at different conditions were prepared by optimized methods [4], and were chemically standardized by HPLC. The herbal extracts of DBT (authentic decoction), AR and ASR were subjected to HPLC spectrum by using UV and ELSD detectors ( Figure S1A,B). In addition, the amounts of key chemicals within the herbal extracts were measured ( Figure S1C): the results were in accordance with previous reports [4,6,7,10].
In cultured osteoblasts, application of DBT induced cell proliferation, as well as the differentiation biomarker ALP, in a dose-dependent manner ( Figure 1A,B). In both scenarios, the DBT-induced osteoblastic growth and differentiation were markedly higher than that of the herbal extract derived from AR + ASR (simple mixing of two herbal extracts), suggesting an uniqueness of DBT formulation, i.e., boiling two herbs together. To detect the amount of collagen being deposited in osteoblasts, Sirius Red stain was used. In cultured osteoblasts, applied DBT was able to enhance the staining density of collagen, significantly, by at least~5-fold, as indicated by dark red clusters of collagens distributed throughout the cultures (Figure 2A). The effect of DBT was more robust than that of simple mixing of AR + ASR. Moreover, the herbal-extract-treated osteoblasts were subjected to determination of extracellular matrix, i.e., identifying the adhesion property. The adhesion property of cultured cells was much less in AR + ASR extract, while DBT showed robust induction as compared to the positive control of dexamethasone and vitamin C ( Figure 2B).

Proteomics Analysis of Herbal-Extract-Treated Osteoblasts
The proteomics profile of cultured osteoblasts, treated with DBT or AR+ASR, was generated by LC-MS, i.e., comparing the profiles of two herbs boiled together (DBT) against that of two herbal extracts simply mixed together (AR+ASR). The data acquired were subjected to sequence searching and label-free quantitative analysis. Label Free Quantitation (LFQ) proteomics, having a calculation of the intensity of precursor ion, was used to reveal the differential expression protein (DEP) after treatment with DBT, as compared with that of AR+ASR. The overall analysis identified 19,466 peptides, mapping to 2800 proteins ( Figure 3A). The "volcano plot" of change against the p-value of LFQ results was constructed: the average spectra number for each comparison was >5 ( Figure 3B). A significantly differentially expressed protein was considered when the p-value of the ttest was <0.05, as well as when the fold of change was higher/lower than ±1.5-fold. From the results, 318 proteins (168 up-regulated and 150 down-regulated) were differentially expressed significantly in DBT-treated osteoblasts, as compared to the AR+ASR group ( Figure 3B). To reveal the difference between DBT and AR+ASR treatment groups, we subjected the normalized protein abundances to principal components analysis (PCA) and heat map clustering. From PCA projection, the maximum variability was identified between DBT vs. AR+ASR ( Figure 3C), the first component covering 41-56% of data variance. In parallel, the heat map showed similar results as that of PCA, where two major clusters separating the protein abundance were identified ( Figure 3D).  Figure 2C). The herbal extract from AR+ASR showed a weak induction of chondrocyte and osteogenic differentiation. In comparison to authentic DBT, the treatment of AR+ASR showed much weaker induction of osteoblastic differentiation in all parameters, suggesting a distinct requirement of boiling the two herbs together during the preparation.

Proteomics Analysis of Herbal-Extract-Treated Osteoblasts
The proteomics profile of cultured osteoblasts, treated with DBT or AR+ASR, was generated by LC-MS, i.e., comparing the profiles of two herbs boiled together (DBT) against that of two herbal extracts simply mixed together (AR+ASR). The data acquired were subjected to sequence searching and label-free quantitative analysis. Label Free Quantitation (LFQ) proteomics, having a calculation of the intensity of precursor ion, was used to reveal the differential expression protein (DEP) after treatment with DBT, as compared with that of AR+ASR. The overall analysis identified 19,466 peptides, mapping to 2800 proteins ( Figure 3A). The "volcano plot" of change against the p-value of LFQ results was constructed: the average spectra number for each comparison was >5 ( Figure 3B). A significantly differentially expressed protein was considered when the p-value of the t-test was <0.05, as well as when the fold of change was higher/lower than ±1.5-fold. From the results, 318 proteins (168 up-regulated and 150 down-regulated) were differentially expressed significantly in DBT-treated osteoblasts, as compared to the AR+ASR group ( Figure 3B). To reveal the difference between DBT and AR+ASR treatment groups, we subjected the normalized protein abundances to principal components analysis (PCA) and heat map clustering. From PCA projection, the maximum variability was identified between DBT vs. AR+ASR ( Figure 3C), the first component covering 41-56% of data variance. In parallel, the heat map showed similar results as that of PCA, where two major clusters separating the protein abundance were identified ( Figure 3D). After revealing all the DEPs, triggered by DBT, the expected and novel signaling pathways were analyzed by KOBAS. The signalling pathways critical for bone development, e.g., cytoskeleton (p = 1.9 × 10 −5 ), tissue development (p = 4.9 × 10 −7 ), osteoblast differentiation (p = 0.05) and the Wnt signaling pathway (p = 0.048), were up-regulated after DBT treatment ( Figure 4A). This observation was consistent with our previous reports of enhancement of osteoblastic function, triggered by DBT [13]. In addition, the DBT-triggered pathways were identified to be related with lipid metabolism, e.g., response to lipid (p = 1.5× 10 −7 ) and fatty acid metabolism (p = 0.02) ( Figure 4B), and glucose and amino acid metabolism, e.g., cellular amide metabolic process (p = 2.14 × 10 −7 ) and cellular component biogenesis (p = 2.3 × 10 −10 ) ( Figure 4C). Moreover, DBT induced a robust RNA metabolism, e.g., proteomes of the structural constituent of ribosomes (p = 4.5 × 10 −6 ) and rRNA processing (p = 8.87 × 10 −8 ) (Figure 4D), and energy metabolism, e.g., oxidoreductase activity (p = 3.6 × 10 −6 ), proton-transporting ATP synthase activity (p = 1.27 × 10 −4 ) and mitochondrial membrane (p = 1.02 × 10 −5 ) ( Figure 4E). Furthermore, other signaling pathways were up-regulated, e.g., calcium ion homeostasis (p = 2.8 × 10 −3 ), the PI3K-Akt signaling pathway (p = 5 × 10 −3 ) and the HIF signaling pathway (p = 1.6 × 10 −2 ) ( Figure 4F). In comparison to DBT, the induced events were not revealed in the treatment of AR+ASR ( Figure 4A-F). After revealing all the DEPs, triggered by DBT, the expected and novel signaling pathways were analyzed by KOBAS. The signalling pathways critical for bone development, e.g., cytoskeleton (p = 1.9 × 10 −5 ), tissue development (p = 4.9 × 10 −7 ), osteoblast differentiation (p = 0.05) and the Wnt signaling pathway (p = 0.048), were up-regulated after DBT treatment ( Figure 4A). This observation was consistent with our previous reports of enhancement of osteoblastic function, triggered by DBT [13]. In addition, the DBT-triggered pathways were identified to be related with lipid metabolism, e.g., response to lipid (p = 1.5× 10 −7 ) and fatty acid metabolism (p = 0.02) ( Figure 4B), and glucose and amino acid metabolism, e.g., cellular amide metabolic process (p = 2.14 × 10 −7 ) and cellular component biogenesis (p = 2.3 × 10 −10 ) ( Figure 4C). Moreover, DBT induced a robust RNA metabolism, e.g., proteomes of the structural constituent of ribosomes (p = 4.5 × 10 −6 ) and rRNA processing (p = 8.87 × 10 −8 ) (Figure 4D), and energy metabolism, e.g., oxidoreductase activity (p = 3.6 × 10 −6 ), proton-transporting ATP synthase activity (p = 1.27 × 10 −4 ) and mitochondrial membrane (p = 1.02 × 10 −5 ) ( Figure 4E). Furthermore, other signaling pathways were up-regulated, e.g., calcium ion homeostasis (p = 2.8 × 10 −3 ), the PI3K-Akt signaling pathway (p = 5 × 10 −3 ) and the HIF signaling pathway (p = 1.6 × 10 −2 ) ( Figure 4F). In comparison to DBT, the induced events were not revealed in the treatment of AR+ASR ( Figure 4A-F). To explore the molecular dynamics of protein-protein interaction from the DEPs, the knowledge-based ingenuity pathway analysis and activated molecular prediction in silico were constructed. Proteome trajectories were categorized into four significant clusters with up-and down-regulated DEPs ( Figure 5). The correlation network was centered in activating multiple collagen synthesis ( Figure 5A), activating β-estradiol ( Figure 5B) and activating ERK1/2 ( Figure 5C). These networks are well recorded to be related to bone metabolism, as predicted by treatment with DBT. To explore the molecular dynamics of protein-protein interaction from the DEPs, the knowledge-based ingenuity pathway analysis and activated molecular prediction in silico were constructed. Proteome trajectories were categorized into four significant clusters with up-and down-regulated DEPs ( Figure 5). The correlation network was centered in activating multiple collagen synthesis ( Figure 5A), activating β-estradiol ( Figure 5B) and activating ERK1/2 ( Figure 5C). These networks are well recorded to be related to bone metabolism, as predicted by treatment with DBT.

Lipidomics Analysis of Herbal-Extract-Treated Osteoblasts
In LC-MS analysis, the dispersed location of MS scans along with retention time and m/z illustrated the successfully optimized LC gradient and efficient separation of the acquisition method, both +ve and −ve modes ( Figure 6A). Total number of features, detected by the centWave method, an algorithm identifying peak density and the waveletbased chromatographic peak from LC/MS data, was used as a measure of metabolome coverage of this combined strategy. After feature alignment, a total number of 5611 and 6690 grouped features were obtained from +ve and −ve ionization profiles, respectively. In addition, the statistical analysis was used to find the probability of a real difference in Processes 2021, 9, 1119 9 of 14 metabolites between experimental groups. After QC filtering, 90% of metabolites had >80% confidence; these values reflected the difference between the two groups.
Processes 2021, 9, x FOR PEER REVIEW 10 of 16 Figure 5. Network analysis of proteomics data. The treatment of osteoblasts was as that in Figure  1. From the proteomics data, the regulated proteins in responding to treatments of DBT and AR+ASR were identified by molecular network analysis. The network was obtained by analyzing the DEPs using ingenuity IPA. The correlation network was centered in (A) collagen synthesis, (B) β-estradiol signaling and (C) ERK1/2 signaling.

Lipidomics Analysis of Herbal-Extract-Treated Osteoblasts
In LC-MS analysis, the dispersed location of MS scans along with retention time and m/z illustrated the successfully optimized LC gradient and efficient separation of the acquisition method, both +ve and -ve modes ( Figure 6A). Total number of features, detected by the centWave method, an algorithm identifying peak density and the wavelet-based chromatographic peak from LC/MS data, was used as a measure of metabolome coverage of this combined strategy. After feature alignment, a total number of 5611 and 6690 grouped features were obtained from +ve and -ve ionization profiles, respectively. In addition, the statistical analysis was used to find the probability of a real difference in metabolites between experimental groups. After QC filtering, 90% of metabolites had >80% confidence; these values reflected the difference between the two groups.
As a quality control procedure, the following steps were routinely performed to ensure system stability and reproducibility of column performance. A mixture of lipid standards was injected as a QC sample before each experiment. The column pressure was recorded and compared to previous runs. The instrument resolution was optimized for different batches. As shown in Figure S2, the QC pool samples were clustered in the center of the PCA plot, which suggested that the differentiation of samples resulted from metabolome differences but not from systematic variance or technical issues. The cumulative value of Q 2 estimates the predictive accuracy of a multivariate statistical model with a Figure 5. Network analysis of proteomics data. The treatment of osteoblasts was as that in Figure 1. From the proteomics data, the regulated proteins in responding to treatments of DBT and AR+ASR were identified by molecular network analysis. The network was obtained by analyzing the DEPs using ingenuity IPA. The correlation network was centered in (A) collagen synthesis, (B) β-estradiol signaling and (C) ERK1/2 signaling.
As a quality control procedure, the following steps were routinely performed to ensure system stability and reproducibility of column performance. A mixture of lipid standards was injected as a QC sample before each experiment. The column pressure was recorded and compared to previous runs. The instrument resolution was optimized for different batches. As shown in Figure S2, the QC pool samples were clustered in the center of the PCA plot, which suggested that the differentiation of samples resulted from metabolome differences but not from systematic variance or technical issues. The cumulative value of Q 2 estimates the predictive accuracy of a multivariate statistical model with a threshold of >0.5, and the difference between Q 2 and R 2 Y should not be larger than 0.3. From multivariate statistical analysis, the fit goodness for PLS-DA and OPLS-DA showed acceptable internal cross-validation results, i.e., PLS-DA: R 2 X cum = 0.820, R 2 Y cum = 0.99, Q 2 cum = 0.998 in +ve mode; R 2 X cum = 0.491, R 2 Y cum = 0.998, Q 2 cum = 0.962 in −ve mode ( Figure 6B); and OPLS-DA: R 2 X cum = 0.591, R 2 Y cum = 0.99, Q 2 cum = 0.99 in +ve mode; R 2 X cum = 0.612, R 2 Y cum = 0.997, Q 2 cum = 0.978 in −ve mode ( Figure 6C). In the score plot of OPLS-DA, the clustering of DBT-treated group was well separated from that of AR+ASR group ( Figure 6C). showed acceptable internal cross-validation results, i.e., PLS-DA: R 2 Xcum = 0.820, R 2 Ycum = 0.99, Q 2 cum = 0.998 in +ve mode; R 2 Xcum = 0.491, R 2 Ycum = 0.998, Q 2 cum = 0.962 in -ve mode ( Figure 6B); and OPLS-DA: R 2 Xcum = 0.591, R 2 Ycum = 0.99, Q 2 cum = 0.99 in +ve mode; R 2 Xcum = 0.612, R 2 Ycum = 0.997, Q 2 cum = 0.978 in -ve mode ( Figure 6C). In the score plot of OPLS-DA, the clustering of DBT-treated group was well separated from that of AR+ASR group ( Figure 6C). The differential identified lipids between DBT vs. AR + ASR groups were selected from the S-plot in terms of the threshold of variable influence on projection (VIP) for OPLS models having VIP values ≥ 1 ( Figure 7A). This indicated that the lipid profiles of osteoblasts at different treatments were very different. The differential expression metabolites (DEMs) were listed in a heat map ( Figure 7B). Compared to AR+ASR, the treatment of DBT led to an increase in ethanolamine, fatty amide, phosphosphingolipid, quinone and sphingolipid, as well as a reduction in phosphatidyl-choline (PC), phosphatidyl-ethanolamine (PE), phosphatidyl-serine (PS) and steroid, significantly. The differential identified lipids between DBT vs. AR + ASR groups were selected from the S-plot in terms of the threshold of variable influence on projection (VIP) for OPLS models having VIP values ≥ 1 ( Figure 7A). This indicated that the lipid profiles of osteoblasts at different treatments were very different. The differential expression metabolites (DEMs) were listed in a heat map ( Figure 7B). Compared to AR+ASR, the treatment of DBT led to an increase in ethanolamine, fatty amide, phosphosphingolipid, quinone and sphingolipid, as well as a reduction in phosphatidyl-choline (PC), phosphatidyl-ethanolamine (PE), phosphatidyl-serine (PS) and steroid, significantly.

Discussion
Bone disease, e.g., osteoporosis, is a threat to human health; this problem is dominant in women aged over 50. Unfortunately, osteoporosis remains practically incurable fully at this moment [14]. One of the most crucial reasons for this unfavorable situation is poor efficiency of drug(s). Under this scenario, combined herbal therapies have been proposed to have enhanced effectiveness and minimized side effects [15]; however, the action mechanism and rationale of herbal combination are unresolved. DBT, a Chinese herbal decoction containing two herbs (AR and ASR), is an excellent example to illustrate the synergy of two herbs. The boiling together of two herbs provides better efficacy in inducing osteoblastic proliferation and differentiation. Furthermore, a complete profile of data and an experimental systems biology approach in revealing the functionality of a herbal mixture were presented here. This systems biology approach can overcome the key barriers in revealing the action mechanism of TCM having a mixture of herbs.
Bone formation is composed of the processes of cellular proliferation, cellular differentiation and organelle development. As revealed here, the proteomic data suggest that the up-regulated proteins, induced by DBT but not by AR+ASR, are enriched in collagen formation, integrin signaling and cytoskeleton and cartilage condensation. The regulation of these proteins is in accord with rapid cytoskeleton dynamics and bone regeneration [16,17]. In addition, the identified networks relating to bone differentiation and proliferation are enriched by DBT, consistent with our previous report [18]. In acting on cultured osteoblasts, there are plenty of differences in the signaling/metabolism of AR+ASR compared with boiling two herbs together, i.e., DBT, as identified by systems biology. The central metabolism, a process essential for bone formation [18,19], was triggered by DBT but not by AR+ASR. In line with this notion, Wnt-induced signaling has been known to trigger glycolysis flux during osteoblastic differentiation [20]. Supporting this notion, an increase in grip strength and swimming time was revealed in DBT-treated rats, i.e., DBT had a positive influence in promoting energy metabolism [21]. Moreover, the glycolysis flux and capacity were markedly enhanced by DBT in cardiomyocytes [22]. Meanwhile, DBT could induce the proton-transporting ATP synthase activity and the morphology of the mitochondrial membrane in cultured osteoblasts, compared to AR+ASR [22]. As a hypothesis, DBT is a potential anabolic agent for clinical disorders of substrate-availability-based osteoporosis. Moreover, the current identified pathways have never been identified before under herbal treatment. For example, the metabolism of amino acids, as induced by DBT but not by AR + ASR, is known to catabolize skeletal muscle cells to produce energy during muscle regeneration [23].
From 1 H-NMR metabolic profiling, asparagine was found in a significantly higher amount in DBT-treated culture than that in AR-or ASR-treated cultures [24]. Asparagine is a key metabolite in the TCA cycle [25]. Therefore, we believe that asparagine may be a key compound to activate energy metabolism in the case of DBT. In parallel, bone cells respond to change in the microenvironment by altering the profile of gene expression. HIF 1-α is one of the critical proteins in this situation, which binds to the hypoxia response element in regulating the expressions of EPO, VEGF and PDGF; these growth factors are related to bone differentiation [26]. In line with this notion, DBT was shown to regulate HIF signaling in cultured HEK293T cells [3].
The importance of the membrane bilayer in promoting the growth and differentiation of osteoblasts has been proposed [27,28]. The lipid profile of culture after DBT treatment was distinct from that of AR+ASR. Interestingly, we found that the amount of acetylcarnitine was increased after DBT treatment. Carnitine and fatty acids are essential molecules for energy metabolism, which convey carbon sources via acetyl-CoA to the citric acid cycle, happening more robustly during inflammation [29]. Besides, acetylcarnitine, triggered by Wnt signaling, was utilized in β-oxidation in osteoblasts [30]. Thus, fatty acid oxidation in osteoblasts is required for bone acquisition in a sex-and diet-dependent manner [30].
In TCM, preparing a herbal decoction in boiling water is the most simple and common practice [31], and this approach has been utilized for thousands of years without change.
DBT is an excellent example to show the requirement of boiling two herbs together to achieve maximal pharmacological efficacy. Supporting this synergy of herbal mixture, other TCM formulae have been demonstrated to have such characteristics. Kai-Xin-San, an ancient herbal decoction containing Ginseng Radix et Rhizoma, Polygalae Radix, Acori Tatarinowii Rhizoma and Poria, can trigger the secretion of neurotrophic factors in cultured astrocytes; this activity is more robust in herbal mixture as compared to a single herb [32].
In TCM, water as a means of herbal processing could be due to the following reasons: (i) water is an excellent solvent for most phytochemicals; and (ii) dissolved phytochemical is easier to be absorbed in the gut. Despite the fact that boiling could increase the solubility of phytochemicals, many unknown chemical changes could happen during the process [33]. In DBT, the AR-derived compounds, including formononetin and calycosin, were demonstrated to play crucial roles in DBT's functions [34]; these chemicals were shown to have higher amounts in DBT, as compared to that from a single herbal extract of AR or ASR. One hypothesis is that calycosin-7-O-β-D-glucoside and ononin can transform to formononetin and calycosin during the boiling process [32]. In line with our results, calycosin is known to stimulate osteogenic differentiation via activation of IGF1R and PI3K/Akt signaling [34], and formononetin can stimulate osteoblastic differentiation via the p38 MAPK pathway [35]. Thus, their signaling pathways are being matched in our DBT-mediated omics analyses.

Conclusions
Here, we have demonstrated an in-depth and integrated pipeline to reveal the requirement of boiling two herbs together, to enhance the osteoblastic function of DBT. First, the integration of multi-omics (proteomics and lipidomics) enables an exploration of molecular perturbation in different levels of osteoblastic development. Second, the osteoblastic function of DBT is related to glycolysis, energy metabolism and lipid metabolism. Third, boiling together the herbs collaborates different compounds within DBT, as to maximize their functions. As a result, the boiling of herbs together has a crucial role in controlling the interactions of DBT components. Our omics approach could serve as a crucial backdrop for future studies that characterize the impact of a broad array of different factors in herbal medicine. In addition, the investigation of action mechanisms in applying multitarget drugs, or combinational therapeutics, should be important for the modernization of herbal medicine.