Abstract

Objective. To investigate the pharmacological mechanism of HuangQiXiXin decoction (HQXXD) on cough variant asthma (CVA) and validate the clinical curative effect. Methods. The active compounds and target genes of HQXXD were searched using TCMSP. CVA-related target genes were obtained using the GeneCards database. The active target genes of HQXXD were compared with the CVA-related target genes to identify candidate target genes of HQXXD acting on CVA. A medicine-compound-target network was constructed using Cytoscape 3.6.0 software, and a protein-protein interaction (PPI) network was constructed using the STRING database. Gene ontology (GO) function enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using RGUI3.6.1 and Cytoscape 3.6.0. We searched the main database for randomized controlled trials of HQXXD for CVA. We assessed the quality of the included studies using the Cochrane Reviewers’ Handbook. A meta-analysis of the clinical curative effect of HQXXD for CVA was conducted using the Cochrane Collaboration’s RevMan 5.3 software. Results. We screened out 48 active compounds and 217 active target genes of HQXXD from TCMSP. The 217 active target genes of HQXXD were compared with the 1481 CVA-related target genes, and 132 candidate target genes for HQXXD acting on CVA were identified. The medicine-compound-target network and PPI network were constructed, and the key compounds and key targets were selected. GO function enrichment and KEGG pathway enrichment analysis were performed. Meta-analysis showed that the total effective rate of the clinical curative effect was significantly higher in the experimental group than the control group. Conclusion. The pharmacological mechanism of HQXXD acting on CVA has been further determined, and the clinical curative effect of HQXXD on CVA is remarkable.

1. Introduction

Cough variant asthma (CVA) has the same pathogenesis as asthma and is characterized by persistent chronic airway inflammation and airway hyperresponsiveness [1]. The clinical manifestations of CVA are chronic cough, which is not always accompanied by wheezing. Chronic cough is often the only symptom of CVA [2]. CVA is the most common cause of chronic cough in China [3] and the second most common cause of chronic cough in Korea [4].

Inhalation of corticosteroids and long-acting β2-agonists is recommended for the treatment of CVA [5, 6]. However, there is a high recurrence rate of the disease, and the long-term effect of these treatments is not ideal. Recently, Chinese herbal medicine has been used in the treatment of CVA [7, 8]. HuangQiXiXin decoction (HQXXD), added or reduced from the traditional Chinese medicine (TCM) prescriptions Yupingfeng powder and Zhisou powder, mainly consists of Huangqi, Xixin, Jingjie, Fangfeng, Huangqin, Baizhu, Fuling, Chantui, Banxia, Baibu, and Gancao. It has been reported that HQXXD is effective in the treatment of CVA; however, our previous study showed that the total effective rate of HQXXD in the treatment of CVA was not significant compared with the control group, although there was some evidence of a clinical effect [9].

In this study, we used network pharmacology and an evidence-based medicine approach to determine the pharmacological mechanism of HQXXD on CVA and validate the clinical curative effect. First, we identified the active compounds and target genes of the monarch and minister herbs in HQXXD. Next, we used the CVA-related target genes and the active target genes of HQXXD to construct the medicine-compound-target network and protein-protein interaction (PPI) network. Subsequently, gene ontology (GO) function enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment were performed. Finally, we screened randomized controlled trials (RCTs) that investigated the clinical curative effect of HQXXD on CVA and performed a meta-analysis to validate the clinical curative effect.

2. Materials and Methods

2.1. Screening of Active Compounds in HQXXD

Huangqi, Xixin, Jingjie, and Fangfeng are monarch and minister herbs, which are regarded as the main active herbs in HQXXD, and have the Latin names Radix Astragali (RA), Herba cum Radix Asari (HRA), Herba Schizonepetae (HS), and Radix Saposhnikoviae (RS), respectively. We have determined that these four herbs are the main components in HQXXD. The compounds in these four herbs were obtained using the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) (http://tcmspw.com/tcmsp.php) [10]. We used the parameters of oral bioavailability (OB) ≥ 30% and drug-likeness (DL) ≥ 0.18 to screen for active compounds, and OB and DL are commonly used in network pharmacology. OB relates to the rate and extent of which the active component or components of the drug are absorbed into the body. DL indicates the degree to which a compound contains specific functional groups or has the same or similar physical characteristics to most drugs.

2.2. Screening of the Target Genes of Active Compounds

We retrieved the corresponding target genes of active compounds from TCMSP. After importing the target genes into the UniProt knowledgebase (http://www.uniprot.org/) and setting the search format as “Homo sapiens,” the human official gene symbols were identified, and the active target genes of HQXXD were screened out.

2.3. Acquisition of CVA-Related Target Genes and Identification of Candidate Target Genes of HQXXD Acting on CVA

Using “cough variant asthma” as the keyword in the GeneCards database (https://www.genecards.org/), we searched and acquired CVA-related target genes. Then, the active target genes of HQXXD were compared with the CVA-related target genes, and the intersecting genes were defined as the candidate target genes for HQXXD acting on CVA.

2.4. Construction of the Medicine-Compound-Target Network and Selection of Key Compounds

We used Cytoscape 3.6.0 software and its Network Analyzer tool function to construct and analyze a medicine-compound-target network [11]. Nodes represented compounds and target genes, and the relationships between them were represented as edges. The key compounds in HQXXD acting on CVA were selected according to the degree of connection between the compound and the target gene.

2.5. Construction of the PPI Network and Selection of Key Targets

The candidate target genes were introduced into the STRING database (https://string-db.org/) to construct a PPI network of HQXXD acting on CVA. The research species was defined as “Homo sapiens,” the lowest interaction score was set to 0.4, and the rest of the parameters were set to the default settings to obtain the PPI network. Using Cytoscape 3.6.0 software and its Network Analyzer tool, the PPI network was subjected to topology analysis, and the key targets of HQXXD acting on CVA were then selected according to the degree values of each target gene in the PPI network.

2.6. GO Function Enrichment and KEGG Pathway Enrichment Analysis

RGUI3.6.1 and org.Hs.eg.db package were applied to obtain the Entrez IDs of candidate target genes. Then, we used RGUI, the clusterProfiler package, and Cytoscape to perform GO function enrichment analysis, which included biological process (BP), molecular function (MF), and cellular component (CC) analysis, and KEGG pathway enrichment analysis [12]. Adjusted was considered as statistically significant.

2.7. Validation of the Clinical Curative Effect of HQXXD on CVA
2.7.1. Screening of RCTs Investigating the Clinical Curative Effect of HQXXD on CVA

We searched using the terms “HuangQiXiXin decoction” and “cough variant asthma” in PubMed, the Cochrane Central Register of Controlled Trials, the Chinese Biomedical Literature database, Chinese National Knowledge Infrastructure (CNKI), the Chongqing VIP database (CQVIP), and Wanfang Data from the establishment of each database to January 15, 2020. When the terms were not included in the titles and abstracts, if the terms were present in the full text, the full-text paper was also screened.

The inclusion criteria included that the study was designed as a RCT, the participants had symptoms that were in accordance with a diagnosis of CVA, HQXXD was used in the experimental group, the control group used conventional therapy without TCM therapy, and the clinical curative effect was the main outcome measure. Exclusion criteria included studies where there were incomplete data on the clinical curative effect. For reports of the same study, the earliest published report was included and the rest were excluded.

2.7.2. Data Extraction and Quality Assessment

Two reviewers independently extracted the information from the included studies. The main information included the first author, year of publication, number of patients in each group, methods of intervention in the experimental and control groups, and outcome data of the clinical curative effect.

We used the Cochrane Reviewers’ Handbook for guidelines to assess the risk of bias. Seven criteria were used: random sequence generation; allocation concealment; patient blinding; assessor blinding; incomplete outcome data; selective outcome reporting; and other risks of bias [13].

These main data were input into the Cochrane Collaboration’s RevMan 5.3 software for meta-analysis to analyze the clinical efficacy of HQXXD on CVA.

2.8. Statistical Analysis

Some of the statistical analyses were performed automatically using the bioinformatic tools of the platforms and software mentioned above. In GO function enrichment and KEGG pathway enrichment, adjusted was considered as statistically significant. When RevMan 5.3 software was used for meta-analysis, the total efficiency rates were expressed as the odds ratio (OR) with 95% confidence interval (CI). Heterogeneity was assessed by the Q test ( value and I2), and indicated heterogeneity across studies. Studies with I2 < 50% were considered to have no heterogeneity and those with I2 ≥ 50% were considered to have heterogeneity. If no heterogeneity was detected, the fixed effects model was used as the pooling method; otherwise, the random effect model was considered to be the appropriate choice [14, 15]. was considered as statistically significant.

3. Results

3.1. Screened Active Compounds in HQXXD

We obtained 87 compounds from RA, 192 compounds from HRA, 159 compounds from HS, and 173 compounds from RS in TCMSP (Supplementary files 1A, B, C, and D). Setting the conditions as OB ≥ 30% and DL ≥ 0.18, 20 active compounds from RA, 8 active compounds from HRA, 11 active compounds from HS, and 18 active compounds from RS were selected. Finally, with exclusion of duplications, 53 active compounds from HQXXD were identified. The basic information on the active compounds from HQXXD is shown in Table 1.

3.2. Screened Active Target Genes of HQXXD

The corresponding target genes of the 53 active compounds were also obtained in TCMSP (Supplementary file 2A). Using the UniProt knowledgebase, the corresponding gene symbols were screened by setting the format as “Homo sapiens” (Supplementary file 2B (available here)). Four compounds (MOL000374, MOL000398, MOL000438, and MOL001506) did not have targets, and one compound (MOL000439) did not have a compatible gene symbol. Finally, 217 active target genes of 48 active compounds from HQXXD were identified (Supplementary file 2C).

3.3. Acquired CVA-Related Target Genes and Identified Candidate Target Genes of HQXXD Acting on CVA

Using “cough variant asthma” as the keyword in the GeneCards database, 1481 CVA-related target genes were acquired (Supplementary files 3A, B).

The 217 active target genes of HQXXD were compared with the 1481 CVA-related target genes, and 132 candidate target genes of HQXXD acting on CVA were identified (Figure 1, Table 2).

3.4. Constructed Medicine-Compound-Target Network and Selected Key Compounds

The 132 overlapping candidate target genes corresponded to 47 active compounds (MOL000433 had correspondence to an overlapping candidate target gene), including 15 compounds from RA, 8 compounds from HRA, 10 compounds from HS, and 16 compounds from RS (Supplementary file 4).

Then, a network of medicine-compound-target was constructed using the Cytoscape software and analyzed by the Network Analyzer tool. There were 184 nodes (4 medicine nodes, 47 compound nodes, 132 target gene nodes, and 1 CVA node) and 653 edges in the network (Figure 2). The top 14 compounds ranked by degree in the network were MOL000098, MOL000006, MOL000422, MOL000173, MOL000378, MOL000358, MOL000392, MOL000371, MOL000449, MOL000354, MOL000380, MOL000417, MOL001460, and MOL002962 (Table 3), which could be considered the key compounds in HQXXD acting on CVA.

3.5. Constructed PPI Network and Selected Key Targets

To further study the mechanism of HQXXD acting on CVA, 132 overlapping target genes were mapped into the STRING database, and the PPI network was obtained. When we set the lowest interaction score to 0.40, 132 target proteins in the network had interactions, and 2373 edges represented the interactions between the proteins (Supplementary file 5, Figure 3). The top 20 target genes ranked by degree in the network are shown in Table 4, which can be considered the key targets of HQXXD acting on CVA.

3.6. GO Function Enrichment and KEGG Pathway Enrichment Analysis

The Entrez IDs of candidate target genes were obtained using RGUI and org.Hs.eg.db (Table 2). Then, GO function enrichment and KEGG pathway enrichment analysis were performed using RGUI and clusterProfiler.

GO BP function enrichment analysis showed that the candidate target genes of HQXXD acting on CVA were significantly enriched in response to toxic substance, response to molecule of bacterial origin, response to lipopolysaccharide, response to extracellular stimulus, response to nutrient levels, response to metal ion, response to reactive oxygen species, response to oxidative stress, cellular response to oxidative stress, response to antibiotic, and so on (Supplementary file 6A). The top 20 GO BP enrichments ranked by the adjusted value are shown in Figure 4(a).

GO MF function enrichment analysis showed that the candidate target genes of HQXXD acting on CVA that were significantly enriched in protein heterodimerization activity, cytokine receptor binding, G protein-coupled amine receptor activity, nuclear receptor activity, transcription factor activity direct ligand-regulated sequence-specific DNA binding, oxidoreductase activity acting on paired donors with incorporation or reduction of molecular oxygen, heme binding, tetrapyrrole binding, cytokine activity, ubiquitin-like protein ligase binding, and so on (Supplementary file 6B). The top 20 GO MF enrichments ranked by the adjusted value are shown in Figure 4(b).

GO CC function enrichment analysis showed that candidate target genes of HQXXD acting on CVA were significantly enriched in membrane raft, membrane microdomain, membrane region, cyclin-dependent protein kinase holoenzyme complex, mitochondrial outer membrane, serine/threonine protein kinase complex, integral component of presynaptic membrane, organelle outer membrane, outer membrane, protein kinase complex, and so on (Supplementary file 6C). The top 20 GO CC enrichments ranked by the adjusted value are shown in Figure 4(c).

KEGG pathway enrichment analysis was further conducted for determining candidate target genes of HQXXD acting on CVA. Candidate target genes were significantly enriched in advanced glycation end products- (AGE-) receptor for the AGE (RAGE) signaling pathway in diabetic complications, fluid shear stress and atherosclerosis, hepatitis B, interleukin- (IL-) 17 signaling pathway, tumor necrosis factor (TNF) signaling pathway, Kaposi sarcoma-associated herpesvirus infection, bladder cancer, prostate cancer, endocrine resistance, pancreatic cancer, and so on (Supplementary file 6D). The top 20 KEGG pathway enrichments ranked by the adjusted value are shown in Figure 4(d).

Using the ClueGO plugin app and CluePedia plugin app of the Cytoscape software, the GO function enrichment and KEGG pathway enrichment from ClueGO can be displayed more intuitively as shown in Figure 5.

3.7. Validation of the Clinical Curative Effect of HQXXD on CVA
3.7.1. Screened RCTs Investigating the Clinical Curative Effect of HQXXD on CVA

A total of 21 studies were retrieved through database searching. After removing duplication, nine studies were retained. A total of four irrelevant studies were excluded after reading the title, abstract, and full text. Five RCTs [1620] were included for further evaluation. The literature screening process and results are shown in Figure 6.

3.7.2. Description of Included RCTs and Assessment of the Methodological Quality

Five eligible RCTs were identified. The five RCTs were all conducted in China and included 342 patients. The five studies were all single-center studies. The basic features of the included studies are outlined in Table 5.

Two RCTs [16, 19] employed adequate methods of random sequence generation; none of the RCTs introduced allocation concealment; none of the RCTs described blindness; all the RCTs had complete outcome data; and for all studies, we were unable to determine if these were selective reports (Figure 7).

3.7.3. Meta-Analysis

The five studies that compared the total effective rate of the clinical curative effect included a total of 342 participants, 171 in the experimental groups and 171 in the control groups. The five studies had homogeneity (heterozygosity test, Chi2 = 1.35, , I2 = 0%). When the fixed effect model was used to merge OR values, the pooled OR was 2.86 (95% CI 1.37–5.95, Z = 2.80, ). This indicated that the total effective rate of the clinical curative effect was significantly higher in the experimental group compared with the control group (Figure 8).

4. Discussion

CVA is characterized by wheezing cough, wind cough, and stubborn cough in TCM, which attributes qi deficiency and vigorous wind as the main causes of this disease and emphasizes that the treatment of patients should be based on this aspect [21].

HQXXD is added or reduced from the TCM prescription Yupingfeng powder and Zhisou powder [9]. RA, HRA, HS, and RS are monarch and minister herbs in HQXXD, which are regarded as the main active herbs in HQXXD. RA is believed to tonify lung qi; HRA, HS, and RS are believed to tonify the spleen and stomach and dispel wind dampness, which strengthen the effect of RA.

Chinese herbal medicines are complex and contain many compounds. At present, most studies of TCM are limited to explaining the mechanism of a certain compound and lack a holistic view. Currently, network pharmacology is widely used in the study of TCM and advocates whole and group analysis, which is consistent with the multicomponent, multitarget, and multipathway characteristics of TCM [2224].

We identified the active compounds and target genes of HQXXD from TCMSP. Then, 132 candidate target genes of HQXXD acting on CVA were identified, and the key compounds of HQXXD acting on CVA were selected. The main key compounds were quercetin, luteolin, kaempferol, wogonin, 7-O-methylisomucronulatol, beta-sitosterol, formononetin, 3,9-di-O-methylnissolin, stigmasterol, and isorhamnetin. Quercetin has a potential protective and preventive effect on the frequency of attacks and in controlling the most common symptoms of asthma [25]; quercetin has both immunomodulatory and bronchodilatory properties [26]. Luteolin has been shown to have an antiallergic effect in a murine model of allergic asthma and rhinitis [27]. Kaempferol may be a potent antiallergic compound targeting the allergic asthma typical of airway hyperplasia and hypertrophy [28]; kaempferol could suppress eosinophil infiltration and airway inflammation in airway epithelial cells and in mice with allergic asthma [29]. Wogonin could attenuate ovalbumin-induced airway inflammation in a mouse model of asthma [30]. Beta-sitosterol possesses antiasthmatic actions by inhibiting the cellular responses and subsequent release/synthesis of Th2 cytokines [31]. Stigmasterol has antiasthmatic properties and suppressive effects on the key features of allergen-induced asthma [32]. The results of our study are in accordance with these previous studies, which indicate that HQXXD has an effect on CVA because of the action of these multiple compounds.

PPI network analysis suggested that the action of HQXXD on CVA is related to multiple targets. The key targets were AKT1, IL6, VEGFA, JUN, EGFR, CXCL8, CASP3, MMP9, PTGS2, and MAPK8. Previous research has indicated that AKT1 is related to pathway mechanisms in children with severe asthma [33]; IL6 polymorphisms are related to the effects of smoking on the risk of adult asthma [34]; VEGFA rs833069 SNPs are associated with asthma [35]; VEGFA is associated with the response to inhaled corticosteroids in children with asthma [36]; smoking aggravates airway inflammation by activating the c-Jun amino terminal kinase pathway in asthma [37]; EGFR activation-induced decreases in claudin 1 promote MUC5AC expression and exacerbated asthma in mice [38]; the EGRF-dependent signaling pathway is associated with diseases, such as asthma [39]; neutrophils release CXCL8, NE, and MMP-9 in response to viral surrogates with R848-induced CXCL8 release being specifically enhanced in neutrophils in asthma [40]. CASP3 may play a role in allergic asthma [41]; the PTGS2 gene is associated with diisocyanate-induced asthma [42]. The relationships between these key targets and CVA, some of which have been well studied, could be further studied, and possible mechanisms that have not been previously studied could be explored.

GO functional enrichment analysis suggested that GO BP for HQXXD in the treatment of CVA were enriched in response to toxic substances, molecules of bacterial origin, lipopolysaccharides, extracellular stimuli, nutrient levels, metal ions, reactive oxygen species, oxidative stress, antibiotics, and cellular responses to oxidative stress. GO MF analysis indicated that HQXXD in the treatment of CVA involved in protein heterodimerization activity, cytokine receptor binding, G protein-coupled amine receptor activity, nuclear receptor activity, transcription factor activity, direct ligand-regulated sequence-specific DNA binding, oxidoreductase activity acting on paired donors with incorporation or reduction of molecular oxygen, heme binding, tetrapyrrole binding, cytokine activity, and ubiquitin-like protein ligase binding. GO CC analysis indicated that HQXXD in the treatment of CVA were enriched in the membrane raft, membrane microdomain, membrane region, cyclin-dependent protein kinase holoenzyme complex, mitochondrial outer membrane, serine/threonine protein kinase complex, integral component of presynaptic membrane, organelle outer membrane, outer membrane, and protein kinase complex. These regions are closely related to airway inflammation, airway epithelial cell injury, airway hyperresponsiveness, and airway remodeling, which are closely related to the pathological mechanisms of CVA [4345].

KEGG enrichment pathway analysis showed that many pathways were closely related to the pathogenesis of CVA. The main pathways included the AGE-RAGE signaling pathway in diabetic complications, fluid shear stress and atherosclerosis, the IL-17 signaling pathway, TNF signaling pathway, and pathways involved in hepatitis B, Kaposi sarcoma-associated herpes virus infection, bladder cancer, prostate cancer, endocrine resistance, and pancreatic cancer. A relationship between these pathways and CVA has also been found or preliminarily indicated. These results indicated that HQXXD acts on CVA through multiple pathways. The main pathways and main target genes in the pathways are worthy of further study.

The pharmacological mechanisms of HQXXD acting on CVA have been investigated through a network pharmacology approach, but there are some limitations with this approach. We identified active compounds and target genes of HQXXD from the TCMSP database and acquired CVA-related target genes from the GeneCards database. Although these databases are the most comprehensive databases currently available, the data collection is still not necessarily comprehensive, and the establishment of the screening criteria for active compounds cannot be completely accurate. In addition, we constructed a PPI network and performed GO function enrichment and the KEGG pathway enrichment analysis to investigate the pharmacological mechanisms of HQXXD acting on CVA. Further investigation of the key target genes and pathways through experimental analysis is required, which will be the focus of future research in our laboratory.

In addition, we used an evidence-based medicine approach to perform a meta-analysis of the clinical curative effect of HQXXD on CVA. Our previous meta-analysis showed that the total effective rate of HQXXD in the treatment of CVA was not significant compared with the control group [9]. We searched the main database and screened RCTs investigating the clinical curative effect of HQXXD in CVA to perform a meta-analysis, which indicated that the total effective rate of the clinical curative effect was significantly higher in the experimental group than in the control group. Thus, an evidence-based medicine approach was used to indicate the efficacy of HQXXD in CVA from a clinical level, which was more meaningful than experimental validation.

5. Conclusion

In conclusion, based on a network pharmacology approach, this study investigated the relationships between the active compounds, target genes, and pathways, and further studied the multicomponent, multitarget, and multipathway characteristics of HQXXD acting on CVA. HQXXD acts on CVA via the effect of multiple compounds related to multiple targets and through multiple pathways. Finally, we investigated the clinical curative effect of HQXXD on CVA using an evidence-based medicine approach, which indicated that the total effective rate of the clinical curative effect was significantly higher in the experimental group compared with the control group. In brief, the pharmacological mechanism of HQXXD acting on CVA has been investigated, and the clinical effect of HQXXD on CVA is remarkable.

Data Availability

The data used to support this study are included within the supplementary files.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Qingqing Xia, Mingtao Liu, and Hui Li contributed equally to this work. Qingqing Xia, Mingtao Liu, Hui Li, and Yufeng Zhang conceived and designed the study and wrote the manuscript. Qingqing Xia, Mingtao Liu, Hui Li, and Lijun Tian were responsible for data collation and extraction and performed the data analysis. Lijun Tian, Jia Qi, and Yufeng Zhang performed supervision and project administration. All authors read and approved the final manuscript.

Acknowledgments

The authors would like to acknowledge the web database platform and software for data analysis. Victoria Muir, PhD, from Liwen Bianji, Edanz Group China (http://www.liwenbianji.cn/ac), edited the English text of a draft of this manuscript. This study was partially supported by the Scientific Research Project from the Health Commission of Wuxi (grant Q202055 to Xia), the Nantong Science and Technology Plan Project (grant MS12017004-2 to Tian), the National Natural Science Foundation of China (grant 81302768 to Qi), and Research Grant of Jiangyin Hospital of Traditional Chinese Medicine (grant to Zhang).

Supplementary Materials

Supplementary file 1: compounds of HQXXD in TCMSP. (A) 87 compounds of RA in TCMSP; (B) 192 compounds of HRA in TCMSP; (C) 159 compounds of HS in TCMSP; (D) 173 compounds of RS in TCMSP. Supplementary file 2: target genes of HQXXD: (A) corresponding target genes of active compounds; (B) corresponding target gene symbols of active compounds; (C) 217 active target genes of HQXXD. Supplementary file 3: CVA-related target genes. (A) CVA-related target genes from GeneCards; (B) 1481 CVA-related target genes. Supplementary file 4: data of the medicine-compound-target network. Supplementary file 5: data of the PPI network. Supplementary file 6. GO function enrichment and KEGG pathway enrichment analysis. (A) BP function enrichment analysis; (B) MF function enrichment analysis; (C) CC function enrichment analysis; (D) KEGG pathway enrichment analysis. (Supplementary Materials)