Integration of transcriptomics, proteomics, metabolomics and systems pharmacology data to reveal the therapeutic mechanism underlying Chinese herbal Bufei Yishen formula for the treatment of chronic obstructive pulmonary disease.

Bufei Yishen formula (BYF) is a traditional Chinese medicine formula, which has long been used as a therapeutic agent for the treatment of chronic obstructive pulmonary disease (COPD). Systems pharmacology has previously been used to identify the potential targets of BYF, and an experimental study has demonstrated that BYF is able to prevent COPD. In addition, the transcriptomic and metabolomic profiles of lung tissues from rats with COPD and BYF‑treated rats have been characterized. The present study aimed to determine the therapeutic mechanisms underlying the effects of BYF on COPD treatment by integrating transcriptomics, proteomics and metabolomics, together with systems pharmacology datasets. Initially, the proteomic profiles of rats with COPD and BYF‑treated rats were analyzed. Subsequently, pathway and network analyses were conducted to integrate three‑omics data; the results demonstrated that the genes, proteins and metabolites were predominantly associated with oxidoreductase activity, antioxidant activity, focal adhesion and lipid metabolism. Finally, a comprehensive analysis of systems pharmacology, transcriptomic, proteomic and metabolomic datasets was performed, and numerous genes, proteins and metabolites were found to be regulated in BYF‑treated rats; the potential target proteins of BYF were involved in lipid metabolism, inflammatory response, oxidative stress and focal adhesion. In conclusion, BYF exerted beneficial effects against COPD, potentially by modulating lipid metabolism, the inflammatory response, oxidative stress and cell junction pathways at the system level.

Introduction

Chronic obstructive pulmonary disease (COPD) is a chronic inflammatory disease, which is characterized by progressive, partially reversible airflow limitation. COPD is considered the third most common life-threatening disease worldwide, and is associated with high morbidity and mortality (1,2). In addition, COPD is considered to be not only a respiratory disease, but also a systemic disorder. Traditional Chinese medicine (TCM) formulas are comprehensive medicinal compounds that may provide a systemic approach to COPD therapy (3).

Bufei Yishen formula (BYF) is a TCM formula, which is composed of 12 medicinal herbs, that has long been used as a therapeutic agent for the treatment of COPD. In our previous clinical study, BYF was reported to exert beneficial effects on measured outcomes in patients with stable COPD over a 6-month treatment period and a 12-month follow-up period (4). Subsequently, a systems pharmacological model was constructed by integrating active compounds prediction, targets prediction and network pharmacology to identify 216 bioactive ingredients from BYF and 195 potential targets. Our previous study demonstrated that BYF was effective for the treatment of rats with COPD and ventricular hypertrophy, due to its inhibitory effects on the expression of inflammatory cytokines and hypertrophic factors, protease-antiprotease imbalance and collagen deposition in vivo (5). However, the systemic mechanism of BYF in the treatment of rats with COPD remains unclear. Therefore, the present study aimed to conduct a systems-level analysis of the therapeutic mechanism of BYF.

High-throughput molecular biological techniques, including transcriptomic, proteomic and metabolomic approaches, have been used to explore complex biological processes and the function of TCM formulas in systems biology. Transcriptomic profiling is a promising approach to analyze the entire genome, which provides details regarding the biological processes underlying respiratory disease development and medical intervention (6). Proteomic profiling has been used to uncover the complexity of the therapeutic effects of TCM formulas by analyzing expressed proteins and protein function in a cellular context (7). Furthermore, metabolomic profiling provides data-rich information regarding the metabolic alterations that occur as a consequence of the transcriptome and proteome, which reflects the genetic, epigenetic, and environmental factors that influence cellular physiology (8). Therefore, combining transcriptomics, proteomics and metabolomics has the potential to provide a system-wide understanding of the complex therapeutic processes of TCM formulas (9,10).

In our previous studies, the transcriptomic and metabolomic profiles of rats with COPD and BYF-treated rats were generated (11,12). The present study aimed to further analyze the molecular mechanisms of BYF on rats with COPD using proteomic datasets. Subsequently, systems pharmacology, transcriptomics, proteomics and metabolomics datasets were integrated, with the aim of providing a system-wide understanding of the molecular mechanisms underlying the therapeutic effects of BYF on rats with COPD.

Materials and methods

Chemicals and animals

Klebsiella pneumoniae (strain ID: 46114) was obtained from the National Center for Medical Culture Collections (Beijing, China). Tobacco (Hongqi Canal® Filter tip cigarette; tobacco type: Tar, 10 mg; nicotine content, 1.0 mg; carbon monoxide, 12 mg) was purchased from China Tobacco Henan Industrial Co., Ltd. (Zhengzhou, China). A total of 32 Sprague-Dawley rats (16 male and 16 female; weight, 200±20 g; age, 6–8 weeks) were obtained from the Experimental Animal Center of Henan Province (Zhengzhou, China). The rats were housed in an animal room at a constant temperature (25±2°C) under a 12-h light/dark cycle with free access to food and water. The present study was approved by the Experimental Animal Care and Ethics Committee of The First Affiliated Hospital, Henan University of Traditional Chinese Medicine (Henan, China), and the methods were conducted in accordance with the approved guidelines of the Experimental Animal Care and Ethics Committee of The First Affiliated Hospital, Henan University of Traditional Chinese Medicine (register no. 2012HLD-0001).

COPD model and drug administration

The COPD rat model and BYF formula were prepared as previously described (13). Briefly, 22 rats (COPD group) were maintained in a closed box and were exposed to tobacco and repeated K. pneumoniae infections. The control group rats were untreated. At the end of week 8, two COPD rats were sacrificed for lung tissue collection, in order to validate that the rat model was successful. The herbal drugs contained within BYF were provided by the Department of Pharmacology, The First Affiliated Hospital, Henan University of Chinese Medicine, and were prepared in fluid extract. The components of BYF were as follows: Ginseng Radix et Rhizoma, 9 g; Astragali Radix, 15 g; Corni Fructus, 12 g; Lycii Fructus, 12 g; Schisandrae Chinensis Fructus, 9 g; Epimedii Herba, 9 g; Fritillariae Thunbergii Bulbus, 9 g; Paeoniae Rubra Radix, 9 g; Pheretima, 12 g; Perillae Fructus, 9 g; Ardisiae Japonicae Herba, 15 g; and Citri Reticulatae Pericarpium, 9 g (5). On week 9, COPD rats were divided to two groups (10 rats each group) and intragastrically treated with normal saline (2 ml) or BYF (4.44 g/kg, 0.5 g/ml) every day between weeks 9 and 20. The control group (10 rats) were also intragastrically treated with normal saline (2 ml) for the same time period. On week 20, all rats were sacrificed, and lung tissues were collected.

Protein expression analysis

Proteins were isolated from the lung tissue from each of the three experimental groups. Briefly, the lung tissues were lysed in lysis buffer [4% SDS, 0.1 M DTT, 0.1 M Tris (pH 8.0)] and homogenized using a mechanical homogenizer (Retsch Technology GmbH, Haan, Germany). The lysates were cleared by centrifugation at 12,000 × g and 4°C for 5 min, prior to storage at −80°C until further use. For proteolytic digestion, trypsin (Roche Diagnostics GmbH, Mannheim, Germany) solution was added to the proteins and incubated for 24 h at 37°C. Subsequently, each of the samples (30 µl) was individually reconstituted with 70 µl isopropanol, vortexed for 1 min at room temperature. Tryptic peptides were labeled with 8-plex isobaric tags (AB Sciex Germany GmbH, Darmstadt, Germany) for relative quantitation according to the manufacturer's protocol.

Strong cation exchange fractionation was performed on a Shimadzu Prominence liquid chromatography system (Kyoto, Japan). Buffers A [10 mM KH2PO4 in 25% acetonitrile (ACN); pH 3] and B (10 mM KH2PO4 and 2 M KCl in 25% ACN; pH 3) were used as the mobile phase. The peptide mixtures were diluted 10 times with buffer A and 100 µl was then loaded onto a PolySULFOETHYL A column (5 µm; 100 Å; 100×4.6 mm i.d.; PolyLC, Columbia, MD, USA). The following gradients was used: 0–1 min, 0–5% B; 1–21 min, 5–30% B; 21–26 min, 30–50% B; 26–31 min, 50% B; 31–36 min, 50–100% B; 36–46 min, 100% B. The flow rate was 1 ml/min and the column temperature was set at 35°C. A total of 23 fractions (2 min/fraction) were collected and desalted with a C18 SPE column (Phenomenex, Torrance, CA, USA). The dried fractions were then dissolved in 0.1% formic acid (FA) for liquid chromatography (LC)-mass spectrometry (MS) analysis.

LC-tandem MS analysis was performed on a Prominence nano LC system (Shimadzu, Kyoto, Japan) coupled on-line to a micrOTOF-Q II mass spectrometer (Bruker Daltonik, Bremen, Germany). Water containing 0.1% FA, and ACN containing 0.1% FA were used as the mobile phases. The peptide samples (~1 µg of each fraction) were loaded onto a pulled tip column (15 cmx100 µm i.d.) packed with C18 Reprosil particles (5 µm; Nikkyo Technos Co., Ltd., Tokyo, Japan). At the flow rate of 300 nl/min, the gradient was as follows: 5–34% B, 25 min; 34–60% B, 5 min; 60–80% B, 5 min; 80% B, 4 min. Mass spectrometry analysis was operated in the positive mode and the ion source settings were as follows: spray voltage, 4500 V; nebulizer pressure, 5 psi; desolvation gas temperature, 200°C. and all MS and MS/MS spectra were obtained in data-dependent mode with one MS full-scan ranging from 300–1,800 m/z followed by 20 MS/MS scans.

The reporter ion ratio for each identified peptide was analyzed by Mascot (v2.2; Matrix Science, Inc., Boston, MA, USA). The proteomics data were analyzed by loess and global median normalization, and then underwent log2 transformation. All statistical analyses were performed using Student's t-test with the SPSS 19.0 software package (IBM Corp., Armonk, NY, USA). P<0.05 was considered to indicate a statistically significant difference. A fold-change >1.0 was considered upregulation, whereas a fold-change <1.0 was considered downregulation.

Gene/protein set enrichment, network and pathway analyses

The molecular function of proteins was explored using the Cytoscape v3.1.1 plugins ClueGO and BiNGO (14,15). Pathway enrichment analysis of transcripts and proteins were analyzed using the Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) database. Regulated pathways were considered statistically significant if the P-value was ≤0.05. For correlation analyses on pathway levels, KEGG pathways, including the pathways of transcripts, proteins, target proteins of system pharmacology and metabolomics, were compared. Pathway enrichment analysis was performed using KEGG and the aforementioned pathways. In addition, Metscape was used to analyze the integrated pathway of gene, protein, target proteins of system pharmacology and metabolomics data (16).

Results and Discussion

Proteomic analysis results of BYF-treated COPD rats

In our previous study (5), systems pharmacology was used to identify the bioactive ingredients, and the potential targets of BYF. In addition, BYF treatment was confirmed to exert beneficial effects on rats with COPD, due to its inhibitory effects on inflammatory cytokine expression, protease-antiprotease imbalance and collagen deposition (5). Furthermore, we identified molecular alterations at the transcriptomic and metabolomic level (11,12). To investigate the system-wide mechanism of BYF in COPD treatment, the present study examined the effects of BYF on the proteomic profiles of lung tissues.

Using an LC-MS-based proteomic analysis, 191 and 195 proteins were revealed to be regulated in the COPD model (vs. the control) and the BYF-treated rats (vs. the COPD model), respectively. According to a further analysis, the 191 proteins regulated in rats with COPD were predominantly associated with phospholipase A2 inhibitor activity, peroxiredoxin activity, oxidoreductase activity and major histocompatibility complex class II protein complex binding, etc. (Fig. 1A). In BYF-treated rats, the 198 regulated proteins were attributed to various molecular functions, including superoxide metabolic process, removal of superoxide radicals and glycolytic process through fructose-6-phosphate (Fig. 1B). These regulated proteins were involved in numerous pathways, including focal adhesion, tight junction, leukocyte transendothelial migration and regulation of actin cytoskeleton (Tables I and II).

Figure 1.

Molecular functions of regulated proteins in lung tissues from rats with COPD and BYF-treated rats. ClueGO was applied to analyze the molecular functions of the regulated proteins in (A) COPD rats and (B) BYF-treated rats. Functionally grouped networks of enriched categories were generated for the regulated proteins. In the functionally grouped networks, terms are presented as linked nodes, and functionally related groups partially overlap; node size represents the significance of term enrichment. (C) Molecular function of the overlapping proteins between COPD rats and BYF-treated rats was analyzed using BiNGO software. Node size is proportional to the number of proteins in the test set. Red indicates the predominant functions of the proteins and the molecular functions associated with BYF-therapeutic effect. BYF, Bufei Yishen formula; COPD, chronic obstructive pulmonary disease.

Table I.

Pathways associated with the proteins regulated in lung tissue from rats with chronic obstructive pulmonary disease.

Term	Count	%	P-value
Glycolysis/Gluconeogenesis	10	0.4179	6.64×10−6
Hypertrophic cardiomyopathy	9	0.3761	5.29×10−5
Dilated cardiomyopathy	9	0.3761	8.70×10−5
Pyruvate metabolism	6	0.2507	3.62×10−4
Glyoxylate and dicarboxylate metabolism	4	0.1672	7.99×10−4
Tight junction	9	0.3761	0.001079
Citrate cycle (tricarboxylic acid cycle)	5	0.2089	0.00125
Leukocyte transendothelial migration	8	0.3343	0.002374
Focal adhesion	10	0.4179	0.003903
Tryptophan metabolism	5	0.2089	0.004803
Valine, leucine and isoleucine degradation	5	0.2089	0.006122
Adherens junction	6	0.2507	0.00651
Cardiac muscle contraction	6	0.2507	0.008115

Table II.

Pathways associated with the proteins regulated in lung tissue from Bufei Yishen formula-treated rats.

Term	Count	%	P-value
Focal adhesion	13	0.5527	8.96×10−5
Tight junction	10	0.4252	3.02×10−4
Hypertrophic cardiomyopathy	8	0.3401	4.79×10−4
Regulation of actin cytoskeleton	12	0.5102	6.92×10−4
Glycolysis/Gluconeogenesis	7	0.2976	0.002728022
Leukocyte transendothelial migration	8	0.3401	0.003040578
Dilated cardiomyopathy	7	0.2976	0.003855419
Adherens junction	6	0.2551	0.007806313
Extracellular matrix-receptor interaction	6	0.2551	0.011328586
Glyoxylate and dicarboxylate metabolism	3	0.1276	0.016548236
Metabolism of xenobiotics by cytochrome P450	5	0.2126	0.017848404
Prion diseases	4	0.1701	0.020435392
Pyruvate metabolism	4	0.1701	0.027198439
Drug metabolism	5	0.2126	0.032368833
Fatty acid metabolism	4	0.1701	0.032961744
Tryptophan metabolism	4	0.1701	0.035013439
Valine, leucine and isoleucine degradation	4	0.1701	0.041555659
Cardiac muscle contraction	5	0.2126	0.041625416

Furthermore, the COPD model group (191 proteins) shared 98 common proteins with the BYF treated-group (195 proteins). Of these 98 proteins, the alterations in the expression of 61 proteins in the COPD model were suppressed by BYF treatment (Table III). These proteins were attributed to numerous biological functions and two pathways, including oxidoreductase activity, peroxiredoxin activity, purine nucleotide binding (Fig. 1C), as well as focal adhesion and the leukocyte transendothelial migration pathway (data not shown).

Table III.

Overlapping proteins between the chronic obstructive pulmonary disease and Bufei Yishen-treated groups.

Accession number	Molecular weight (kDa)	log2(A/B)	A/B	log2(B/F)	B/F
IPI00190577	404	0.1	1.071773	−0.1	0.933033
IPI00191728	48	0.3	1.231144	−0.1	0.933033
IPI00192301	22	0.1	1.071773	−0.1	0.933033
IPI00193716	46	0.3	1.231144	−0.2	0.870551
IPI00194097	54	0.3	1.231144	−0.4	0.757858
IPI00195516	51	0.3	1.231144	−0.2	0.870551
IPI00196994	23	0.3	1.231144	−0.1	0.933033
IPI00197770	56	−0.1	0.933033	0.1	1.071773
IPI00198887	57	−0.3	0.812252	0.2	1.148698
IPI00200593	47	0.1	1.071773	−0.2	0.870551
IPI00201300	44	0.1	1.071773	−0.1	0.933033
IPI00201561	22	0.1	1.071773	−0.1	0.933033
IPI00203214	95	0.2	1.148698	−0.5	0.707107
IPI00205135	77	0.3	1.231144	−0.2	0.870551
IPI00205332	35	0.8	1.741101	−0.5	0.707107
IPI00206403	38	−0.3	0.812252	0.4	1.319508
IPI00207014	45	0.1	1.071773	−0.2	0.870551
IPI00207146	16	−0.2	0.870551	0.3	1.231144
IPI00208422	88	0.3	1.231144	−0.3	0.812252
IPI00209113	226	−0.1	0.933033	0.1	1.071773
IPI00211448	61	−0.1	0.933033	0.2	1.148698
IPI00212314	68	0.4	1.319508	−0.3	0.812252
IPI00212523	20	−0.6	0.659754	0.2	1.148698
IPI00214457	21	−0.3	0.812252	0.3	1.231144
IPI00215564	24	0.1	1.071773	−0.1	0.933033
IPI00230787	29	−0.3	0.812252	0.5	1.414214
IPI00231423	64	−0.1	0.933033	0.5	1.414214
IPI00231643	16	0.1	1.071773	−0.3	0.812252
IPI00231694	146	0.3	1.231144	−0.5	0.707107
IPI00231825	16	0.2	1.148698	−0.2	0.870551
IPI00231925	41	−0.5	0.707107	1.1	2.143547
IPI00324986	51	0.1	1.071773	−0.4	0.757858
IPI00325189	17	0.3	1.231144	−0.1	0.933033
IPI00326140	167	−0.5	0.707107	0.3	1.231144
IPI00326179	50	0.3	1.231144	−0.1	0.933033
IPI00326972	62	0.5	1.414214	−0.2	0.870551
IPI00327502	26	−0.2	0.870551	0.2	1.148698
IPI00337168	58	0.6	1.515717	−0.3	0.812252
IPI00358087	37	−0.2	0.870551	0.3	1.231144
IPI00360930	28	−0.4	0.757858	0.8	1.741101
IPI00362072	45	−0.2	0.870551	0.2	1.148698
IPI00362755	204	−0.1	0.933033	0.4	1.319508
IPI00363395	21	0.2	1.148698	−0.1	0.933033
IPI00364890	42	0.6	1.515717	−0.4	0.757858
IPI00365929	49	−0.1	0.933033	0.6	1.515717
IPI00368347	118	0.8	1.741101	−0.4	0.757858
IPI00370654	29	0.3	1.231144	−0.3	0.812252
IPI00371990	27	−0.3	0.812252	0.2	1.148698
IPI00372839	110	0.3	1.231144	−0.1	0.933033
IPI00373505	39	−0.2	0.870551	0.3	1.231144
IPI00388249	80	−0.6	0.659754	0.6	1.515717
IPI00389571	54	0.1	1.071773	−0.1	0.933033
IPI00392216	68	0.1	1.071773	−0.1	0.933033
IPI00411230	26	−0.1	0.933033	0.2	1.148698
IPI00421428	29	0.4	1.319508	−0.4	0.757858
IPI00421517	53	−0.1	0.933033	0.2	1.148698
IPI00470288	43	0.3	1.231144	−0.1	0.933033
IPI00471584	83	−0.1	0.933033	0.1	1.071773
IPI00480679	48	−0.7	0.615572	0.4	1.319508
IPI00948384	146	0.4	1.319508	−0.2	0.870551
IPI01016479	375	0.7	1.624505	−0.1	0.933033

Association between genes, proteins and metabolites

In our previous study, 18 samples were randomly chosen from three experimental groups for gene expression experiments. According to a cut-off value of P<0.05, 2,463 and 2,292 differentially expressed genes were detected between the control and COPD groups, and the COPD and BYF-treated groups, respectively (11). In addition, 49 and 31 regulated metabolites were detected in lung tissues from COPD rats and BYF treated-rats, respectively (12). The present study aimed to provide a system-wide view of the therapeutic mechanism of BYF in COPD treatment by integrating transcriptomics, proteomics and metabolomics data.

Metscape software was used to investigate the latent relationships between the gene, protein and metabolite measurements. Initially, two gene-metabolite networks were generated based on the transcriptomics and metabolomics data of COPD and BYF-treated rats. As shown in Fig. 2A and B, the metabolite-gene networks were mainly associated with lipid and purine metabolism. Subsequently, protein-metabolite networks were generated using the metabolomics and proteomics data from COPD and BYF-treated rats (Fig. 2C and D). The results demonstrated that these proteins and metabolites were also predominantly associated with lipid and purine metabolism. Furthermore, the majority of metabolites were involved in lipid metabolism in the gene/protein-metabolite networks. These findings suggested that lipid metabolism may be the critical biological process associated with COPD development and medical intervention.

Figure 2.

Correlation analysis of the metabolites, genes and proteins regulated in lung tissues of COPD rats and BYF-treated rats. Compound reaction networks of the metabolites, genes and proteins were visualized using Metscape: Metabolites (hexagons) and metabolic enzymes (circles) are presented as nodes and reactions are presented as edges. Inputted genes and proteins are shown in red, inputted metabolites are shown in blue. (A) Metabolite-gene network of the COPD model group. (B) Metabolite-gene network of the BYF-treated group. (C) Metabolite-protein network of the COPD model group. (D) Metabolite-protein network of the BYF-treated group. BYF, Bufei Yishen formula; COPD, chronic obstructive pulmonary disease.

Comprehensive analysis of systems pharmacology, transcriptomics, proteomics and metabolomics data

Systems pharmacology was previously applied to identify the active compounds and potential targets of BYF (5). To provide a more in-depth understanding of the systemic mechanism of BYF in treating COPD rats, the systems pharmacology, transcriptomics, proteomics and metabolomics data were integrated.

Initially, the direct correlation between the potential targets and transcripts were analyzed. A total of 10 overlapping proteins [aldehyde dehydrogenase, mitochondrial (ALDH2); cell division protein kinase 4; checkpoint kinase 1; dopamine receptor D2; glucose-6-phosphate 1-dehydrogenase (G6PD); glutathione S-transferase Mu (GSTM)2; interleukin (IL)-1β; vascular endothelial growth factor receptor 2; glucocorticoid receptor; prostaglandin G/H synthase (PTGS)1] between the targets of BYF and the transcripts regulated in the BYF-treated group were detected, which could be attributed to various molecular functions, including oxidoreductase activity, nuclear factor (NF)-κB binding and IL-1 receptor binding (Fig. 3A). Subsequently, 8 overlapping proteins [GSTM1; GSTM2; glutathione S-transferase P; superoxide dismutase (Cu-Zn) (SOD1); xanthine dehydrogenase/oxidase; 78 kDa glucose-regulated protein; ALDH2; ATP synthase subunit β, mitochondrial] were identified between the potential targets and the proteins regulated in BYF-treated rats (Fig. 3B). These proteins were primarily involved in oxidoreductase and antioxidant activity. Subsequently, the latent correlation between the target proteins and metabolites regulated in the BYF-treated group was determined using Metscape software (Fig. 3C). The results demonstrated that the metabolites and target proteins were predominantly involved in lipid metabolism.

Figure 3.

Molecular functions of the overlapping proteins, and the correlation between the metabolites regulated in BYF-treated rats and target proteins. The molecular functions of the overlapping proteins were analyzed using BiNGO. The area of a node was proportional to the number of proteins in the test set. (A) Molecular function of overlapping proteins between the potential targets and transcript measurements in lung tissues of BYF-treated rats. (B) Molecular function of overlapping proteins between the potential targets and proteome measurements in lung tissues of BYF-treated rats. Red indicates the predominant functions of the proteins, and the molecular functions associated with BYF-therapeutic effect. (C) Metscape software was used to build a compound reaction network: Metabolites (hexagons) and metabolic enzymes (circles) are presented as nodes and reactions are presented as edges. Inputted metabolites are shown in blue, inputted target proteins are shown in red. BYF, Bufei Yishen formula.

Based on the integrated analysis, a comprehensive image of the therapeutic mechanism of BYF in COPD treatment was generated (Fig. 4). The comprehensive image predominantly consisted of four groups: Lipid metabolism, inflammatory response, oxidative stress and focal adhesion.

Figure 4.

Comprehensive analysis of the potential targets of BYF, and the regulated genes, proteins and metabolites in the lung tissues of BYF-treated rats. The potential targets, genes, proteins and metabolites are presented as different-colored rectangles. Regulation is color-coded: Red arrow indicates upregulation, blue arrow indicates downregulation, and a gray rectangle represents no regulation. BYF, Bufei Yishen formula.

Our previous study demonstrated that BYF achieved its ameliorative effect on rats with COPD by suppressing the expression of inflammatory cytokines, including IL-1β, IL-6 and tumor necrosis factor-α (5). The systems pharmacology results indicated that extracellular signal-regulated kinase (ERK), p38, c-Jun N-terminal kinase (JNK) and NF-κB were the potential targets of the active compounds contained in BYF. These data suggested that BYF may decrease the expression levels of inflammatory cytokines, potentially by regulating the activation of ERK, p38, JNK and NF-κB.

In our metabolomics study, it was demonstrated that arachidonic acid metabolism was a significantly dysregulated pathway; with the metabolites linoleate, dihomo-γ-linoleate, lipoxin (LX)A4, LXB4, leukotriene (LT)B4, and prostaglandin (PG)E2, decreased by BYF treatment (12). Specifically, LXA4, LXB4, LTB4 and PGE2 participated in inflammatory processes in the airways of patients with COPD (17,18), and metabolic enzymes, including arachidonate 5-lipoxygenase, PTGS1/2, LTA4 hydrolase and aldo-keto reductase family 1 member C3, were potential targets of BYF. Therefore, these findings indicated that BYF achieved its anti-inflammatory activity probably through suppressing lipid metabolism, including arachidonic acid metabolism.

Oxidative stress can trigger sustained inflammatory responses and is the major contributing factor to obstructive lung disorders (19,20). In the present study, 8 overlapping proteins were identified between the targets of BYF and proteomic measurements of the BYF-treated group, which were predominantly involved in oxidative stress. For example, the potential targets (glutamate-cysteine ligase catalytic subunit; glutathione reductase, mitochondrial; G6PD; glutathione S-transferase P; glutathione S-transferase A1/2; GSTM1/2), were involved in glutathione metabolism (21,22). In addition, the levels of antioxidant proteins, including SOD1, which were involved in the pathogenesis of COPD, were increased by BYF treatment (23,24). These results suggested that regulating oxidative stress status may be one of the main causes of the anti-inflammatory activity of BYF.

In addition, focal adhesion, which is an overlapping pathway among proteomic measurements and the potential targets of BYF (Tables II and IV), was significantly regulated by BYF treatment. Specifically, the activation of mitogen-activated protein kinase 1/3 and 8, and JUN, which are potential targets of BYF, results in upregulation of the transcript levels of preinflammatory cytokines (25–27). Taken together, the present study demonstrated that BYF provided therapeutic benefits against COPD through modulating numerous biological functions, including lipid metabolism, oxidative stress, inflammatory response and focal adhesion pathway.

Table IV.

Associated pathways of the potential targets of Bufei Yishen formula.

Term	Count	%	P-value
Neuroactive ligand-receptor interaction	30	0.7413	P<0.0001
Amyotrophic lateral sclerosis	12	0.2965	P<0.0001
Pathways in cancer	27	0.6672	P<0.0001
Drug metabolism	12	0.2965	P<0.0001
Calcium signaling pathway	19	0.4695	P<0.0001
Bladder cancer	10	0.2471	P<0.0001
Metabolism of xenobiotics by cytochrome P450	11	0.2718	P<0.0001
Non-small cell lung cancer	10	0.2471	P<0.0001
Glutathione metabolism	9	0.2224	0.0001
Colorectal cancer	11	0.2718	0.0001
Small cell lung cancer	11	0.2718	0.0001
Pancreatic cancer	10	0.2471	0.0001
Prostate cancer	11	0.2718	0.0001
Vascular endothelial growth factor signaling pathway	10	0.2471	0.0002
Gap junction	10	0.2471	0.0007
Thyroid cancer	6	0.1483	0.001
Focal adhesion	15	0.3706	0.0011
Gonadotropin-releasing hormone signaling pathway	10	0.2471	0.0013
Glioma	8	0.1977	0.0015
Alzheimer's disease	13	0.3212	0.0015
Progesterone-mediated oocyte maturation	9	0.2224	0.0022
Prion diseases	6	0.1483	0.0024
T cell receptor signaling pathway	10	0.2471	0.0026
Arginine and proline metabolism	7	0.173	0.003
Melanoma	8	0.1977	0.003
Fc epsilon RI signaling pathway	8	0.1977	0.0051
Toll-like receptor signaling pathway	9	0.2224	0.0061
Neurotrophin signaling pathway	10	0.2471	0.0065
Nucleotide-binding oligomerization domain-like receptor signaling pathway	7	0.173	0.0065
Apoptosis	8	0.1977	0.0093
Adipocytokine signaling pathway	7	0.173	0.0095
Vascular smooth muscle contraction	9	0.2224	0.0111
Insulin signaling pathway	10	0.2471	0.0111
Renal cell carcinoma	7	0.173	0.0117
Endometrial cancer	6	0.1483	0.0132
Mitogen-activated protein kinase signaling pathway	15	0.3706	0.0137
Caffeine metabolism	3	0.0741	0.0142
B cell receptor signaling pathway	7	0.173	0.0161
Arachidonic acid metabolism	6	0.1483	0.0178
Graft-versus-host disease	5	0.1235	0.021
Phenylalanine metabolism	4	0.0988	0.021
Tryptophan metabolism	5	0.1235	0.0228
Oocyte meiosis	8	0.1977	0.0302
ErbB signaling pathway	7	0.173	0.031
Tyrosine metabolism	5	0.1235	0.0312
p53 signaling pathway	6	0.1483	0.0374
Epithelial cell signaling in Helicobacter pylori infection	6	0.1483	0.0374
Type II diabetes mellitus	5	0.1235	0.0385
Complement and coagulation cascades	6	0.1483	0.0395

Our previous systems pharmacology study identified the active compounds and potential targets of BYF, and demonstrated that BYF treatment was able to exert therapeutic effects against rats with COPD (5). In addition, the transcriptomic and metabolomic profiles of lung tissues derived from COPD and BYF-treated rats were generated. In the present study, a systems biology approach was used to analyze the therapeutic mechanism of BYF in treating COPD rats. Initially, proteomic profiles were obtained from the lung tissues of COPD and BYF-treated rats, and the three levels of omics data were then integrated. The gene/protein-metabolite model was generated from the regulated genes, proteins and metabolites, and was used to identify the affected pathways and to examine them according to the measured abundances of genes, proteins and metabolites. The results indicated that these transcripts, proteins and metabolites were attributed to various functions, including oxidoreductase activity, antioxidant activity and lipid metabolism. Subsequently, a comprehensive method was used to integrate the systems pharmacology and 3-omics data. The system-wide findings demonstrated that BYF potentially achieved its therapeutic effects over COPD rats by regulating lipid metabolism, the inflammatory response, oxidative stress and focal adhesion pathways at the systems level. In conclusion, the present study suggested that an integrated systems pharmacology and transcriptomics, proteomics and metabolomics approach has the potential to considerably advance understanding regarding the therapeutic mechanism of TCM.