Integrating Network Pharmacology and RT-qPCR Analysis to Investigate the Mechanisms Underlying ZeXie Decoction-Mediated Treatment of Non-alcoholic Fatty Liver Disease

ZeXie Decoction (ZXD) is a traditional Chinese medicine composed of Alisma orientalis (Sam.) Juzep. and Atractylodes macrocephala Koidz. ZXD has been widely used to treat non-alcoholic fatty liver disease (NAFLD). The mechanistic basis for the pharmacological activity of ZXD, however, remains poorly understood. In this study, we used a network pharmacology approach and investigated the association between ZXD and NAFLD. We identified the active ingredients of ZXD and screened the potential targets of these ingredients, after which a database of relevant NAFLD-related targets were constructed and several enrichment analyses were performed. Furthermore, the ethanol and aqueous extracts of ZXD were prepared and experimental pharmacology validation was conducted using RT-qPCR of the non-alcoholic fatty liver disease (NAFLD) model in Sprague-Dawley (SD) rats. As a result, a herb-compound-target-pathway network model was developed, and HMGCR, SREBP-2, MAPK1, and NF-κBp65 targets were validated. The gene expression results of these four targets were consistent with those of the network pharmacology prediction. Using an integration strategy, we revealed that ZXD could treat NAFLD by targeting HMGCR, SREBP-2, MAPK1, and NF-κBp65.


INTRODUCTION
Rising obesity rate has partially been linked to the increasing prevalence of non-alcoholic fatty liver disease (NAFLD), characterized by hepatic steatosis unrelated to alcohol consumption (Than and Newsome, 2015). More than 25% of people suffer from NAFLD worldwide. A significant cause of human morbidity, NAFLD, can progress to severe conditions, including fibrosis, cirrhosis, and hepatoma (Tandon, 2013;Angulo et al., 2015). Treatment options for NAFLD are limited and are centered mainly around bariatric surgery and lifestyle changes.
TCMFs are typically designed as multi-component multitarget pharmaceutical agents, in contrast to modern drugs typically designed to target a single protein . Presently, the molecular mechanisms underlying the efficacy of many TCMFs remain uncertain. In the past, experimental model systems were leveraged for pharmacological studies of TCMFs. However, developments in systems biology have led to more widespread use of network pharmacology-based analyses of traditional medicine preparations (Lee et al., 2018;Zuo et al., 2018;Huang et al., 2020). Network pharmacology tools offer an ideal opportunity to understand the complex mechanisms through which TCMFs function in vivo.
In our previous work , we studied ethanol and aqueous extracts of ZXD using ultraperformance liquid chromatography coupled with mass spectrometry to identify the major components of ZXD. Based on this work, published literature, and database sources, we developed a herb-compound-target-pathway network model and performed enrichment analysis to investigate the relationship between ZXD and NAFLD treatment. The network pharmacology framework is depicted in Figure 1. Basing on the non-alcoholic fatty liver disease (NAFLD) model in Sprague-Dawley (SD) rats, the prediction of the mechanism was validated by RT-qPCR analysis in vivo.  The CID of compound was obtained from PubChem database.
FIGURE 2 | The venn diagram of targets in ZX, BZ and NAFLD. Totally 1575 NAFLD-related targets were identified, of which 122 targets overlapped with our list of ZXD targets.
FIGURE 3 | The herbs-compounds-targets network of 143 nodes. These 19 active ingredients and 122 targets were incorporated into a herbs-compounds-targets network wherein node size was scaled according to the degree values. The light-yellow nodes represent ZX and its compounds, the light blue nodes represent BZ and its compounds and the light green nodes represent the gene symbols. The size of node is proportional to the degree value and the edges represent interactions between the nodes.

Experimental Animals
A total of 72 male SD rats weighting 200 ± 20 g (Beijing Vital River Laboratory Animal Technology Co. Ltd.) were housed in a room at a constant temperature of 23 ± 3°C, a constant humidity of 50 ± 10%, under a 12 h light/dark cycle. The rats were fed UVdisinfected fodder and drinking water. The animal experiment was approved by the Laboratory Animal Center, Institute of Medicinal Plant Development, Chinese Academy of Medical Sciences and Peking Union Medical College. The rats were adaptively bred for 1 week in the facility before the experiments.

Reagents
The dried rhizomes of ZX (Batch number: DD6081) and BZ (Batch number: DD8061) were purchased from Beijing Huamiao Pharmaceutical Co., Ltd. (Beijing, China). Their qualities were in accordance with the Chinese Pharmacopoeia. The specimens of ZX and BZ have been preserved in the Institute of Medicinal Plant Development, Chinese Academy of Medical Sciences and Peking Union Medical College.

Sample Preparation
The samples of ZXD were prepared by combining ZX and BZ in a 5:2 ratio (w/w). The samples were crushed and refluxed thrice with 10-times the volume of 95% ethanol for 2 h. The filtrates were then pooled and dried using a rotary evaporator under vacuum, yielding 14.29% ethanol extract.
A Q-TOF analyzer in the SYNAPT G2 HDMS system (Waters, Manchester, United Kingdom) was used for tandem MS analyses. Positive ionization mode was used to acquire ESI spectra in the m/z 50-1,200 mass range with the following settings: a 120°C source temperature; 3.0 kV as the chosen capillary voltage; a desolvation gas temperature of 450°C and a gas flow rate of 800 L/ h; a 40.0 V sampling cone voltage and a 50 L/h cone gas rate; 10-40 eV as the collision energy. The Masslynx 4.1 software package was used to acquire all data.

Bioactive Ingredient Identification
The widely used Traditional Chinese Medicine System Pharmacology database and Analysis Platform (Ru et al., 2014) (TCMSP, https://tcmspw.com/tcmsp.php) was utilized to screen the active ingredients within ZX, BZ, and ZXD. In total, 101 herbal ingredients were identified using the keywords "Alisma Orientale (Sam.) Juzep." and "Atractylodes Macrocephala Koidz.". The potential active ingredients were then selected using an oral bioavailability (OB) threshold of ≥30% and a drug-likeness (DL) threshold of ≥0.18 (Leeson and SPRINGTHORPE, 2007;Xu et al., 2012). In addition, any compounds in ZXD, ZX, or BZ that have previously been reported to be highly pharmacologically active were also included in this analysis . Two-dimensional (2D) structures of all compounds were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) and used for target prediction efforts.

Target Prediction
TheSwissTargetPrediction platform (David et al., 2013) (http://www.swisstargetprediction.ch/index.php) was utilized to identify the putative targets of the compounds retrieved from the PubChem database. The Swiss library contains over 3,000 proteins and 370,000 active ingredients, making it ideal for such efforts. Files of the identified bioactive ingredients with the suffix "SDF" were used as inputs for the SwissTargetPrediction platform, and probable macromolecular targets of these compounds were identified based on structural similarity. Only the SwissTargetPrediction platform was used for this predictive analysis to ensure the clarity of results. The species was set to Homo sapiens for this analysis and all subsequent analyses.

NAFLD-Related Target Identification
Annotated and predicted NAFLD-related targets were identified using the Online Mendelian Inheritance in Man (OMIM, https://omim.org/) (Hamosh et al., 2002) and GeneCards (https://www.genecards.org/) (Stelzer et al., 2016) databases. We used "Non-alcoholic fatty liver" and "Non-alcoholic fatty liver disease" as search terms and integrated the search results to develop a NAFLD-specific disease target database. NAFLD-related targets that were also identified as putative targets of ZXD bioactive ingredients were retained for further analysis.

Network Construction and Enrichment Analyses
Cytoscape 3.7.2 (Shannon et al., 2003) (https://cytoscape.org/) was used to construct a herb-compound-target network, with nodes corresponding to herbs, compounds, and targets and edges representing between-node interactions. Node sizes were scaled based on degree values determined using Cytoscape Network Analyzer such that smaller nodes were associated with smaller degree values.
The potential functional roles of targets overlapping between ZXD and NAFLD target lists were evaluated via Metascape (Zhou et al., 2019) (https://metascape.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, and Gene Ontology (GO) enrichment analyses. All genes were used as the enrichment background for KEGG analyses, with significantly enriched terms meeting the following criteria: count ≥3, p < 0.01, enrichment factor >1.5. These terms were FIGURE 8 | The herb-compounds-targets-pathways network of ZX. The light red nodes represent gene symbols, the light-yellow nodes represent pathways, the light green nodes represent compounds and the light blue nodes represent ZX. The size of node of gene symbol is proportional to the degree value and the edges represent interactions between the nodes.
Frontiers in Pharmacology | www.frontiersin.org September 2021 | Volume 12 | Article 722016 8 grouped and clustered based on the membership similarity. For GO analyses, the molecular functions (MFs), biological processes (BPs), and cellular components (CCs) for which these targets were enriched were identified and visualized using enrichment plots.
The ZX and BZ herb-compound-target-pathway networks were then plotted based on the results of the above enrichment analyses.

Protein-Protein Interaction Analysis
Protein-protein interaction (PPI) analysis was conducted to identify which ZXD targets are most crucial in the context of disease pathology. STRING 11.0 (Szklarczyk et al., 2017) (https:// string-db.org/) database was used to conduct PPI analysis of overlapping ZXD and NAFLD targets and KEGG pathway analysis. For the PPI analysis, active interaction sources were text mining, experiments, databases, co-expression, neighborhood, gene fusion, and co-occurrence. The analysis range was only set to input proteins, and output obtained in the form of files with the suffix "TSV" was visualized using Cytoscape 3.7.2; NetworkAnalyzer was used to assess the degree and combined scores.
All macromolecular targets were retrieved from the Uniprot database (Alex, 2019) (https://www.uniprot.org/) to confirm their gene symbols. The plots described above were designed based on these symbols to ensure statistical normalization.

Experimental Design
Seventy-four male SD rats were randomly divided into two groups according to their diets, including blank (n 13) and model group (n 61). The blank group rats were fed with normal chew diets (NCD, Beijing Hfk Bioscience Co., Ltd.), and the model group rats were fed with high-fat diets containing 41% fat (HFD, Beijing Hfk Bioscience Co., Ltd.) for 12 weeks to induce NAFLD. After 12 weeks of modeling, one rat liver tissue from each group was used for histopathological evaluation to validate the NAFLD model. Afterwards, a 4-weeks intervention period with ZXD and positive control drug started. Rats from model group were randomly and equally assigned to five subgroups: vehicle, positive control, low-dose, middle-dose and high-dose groups (n 12 each). The vehicle group and positive control group rats were orally administered distilled water and 10 mg/kg/d PH, respectively. The low-dose, middle-dose and high-dose group rats were orally administered 100, 200 and 400 mg/kg/d ZXD, respectively. The doses of ZXD and PH were selected based on the previous studies (Hong et al., 2010;Gong et al., 2018;Yang et al., 2019). During the 4-week intervention period, the blank group rats were fed with NCD and the remaining five group rats were fed with HFD. Animals were given free access to feed and drink. The animal weights and food consumption were monitored daily. At the end of FIGURE 9 | The PPI network of all matched targets. The size of nodes represents their degree values and the width of edges represents interactions between the nodes. The color represents their degree values and combined scores, which changes from blue to orange represent values from low to high. intervention period, animals were fasted for 16 h and then anesthetized with an intraperitoneal injection of urethane (20%, 1 g/kg). The biological samples (blood, liver and urine) were collected and stored at -80°C to conduct the histopathological evaluation and serum biochemical assay.

Histopathological Evaluation Hematoxylin and Eosin (HE) Staining
Fresh liver tissues were fixed with 4% paraformaldehyde for 24 h. The fixed tissues were then embedded in paraffin, cut into thin slices and stained with HE. Light microscopy was used to determine the histopathological changes.

Oil Red O (ORO) Staining
Frozen liver tissues were cut into thin sections and then mounted on glass slides. The slides were fixed in pre-cold 4% paraformaldehyde and rinsed with distilled water. The rinsed slides were placed in isopropanol solution and then incubated in ORO solution for 10 min. Afterwards, the ORO solution was removed and the slides were rinsed with running tap water for 5 min. Light microscopy was used to determine the histopathological changes.

Serum Biochemical Assays
The blood samples were collected and left standing for 30 min. The samples were centrifuged for 10 min at 3,000 g to separate the serum for biochemical assays. The total cholesterol (TC), triglycerides (TG), aspartate aminotransferase (AST), alanine aminotransferase (ALT) levels were measured according to the kit instructions.

RT-qPCR
RT-qPCR was then performed to measure the hepatic expression of the predicted NAFLD-related genes. Total RNA was collected with the RNA isolation kit according to the instructions. The purity of the collected RNA was detected by the NanoDrop ND-2000 spectrophotometer, FIGURE 10 | The PPI network of insulin-related targets. The size of nodes represents their degree values and the width of edges represents interactions between the nodes. The color represents their degree values and combined scores, which changes from blue to orange represent values from low to high.  Table 1.

Statistical Analysis
Statistical analysis was performed using SPSS 23.0 software, and the results are presented as mean ± SD. One-way analysis of variance was used for comparison, and p-value＜0.05 was considered statistically significant. The data graphs were generated using GraphPad Prism 8.

Active Ingredients of ZXD
In our initial search of the TCMSP database, we identified 101 active ingredients of potential relevance. Of these, we selected six ingredients in ZX (Z1 to Z6) and four ingredients in BZ (B4 to B7) based upon their OB and DL values. In addition, we selected Frontiers in Pharmacology | www.frontiersin.org September 2021 | Volume 12 | Article 722016 11 nine ingredients (Z7 to Z12 and B1 to B3) that have previously been reported to exhibit a high degree of pharmacological activity in ZXD preparations (Lin, 2012;Xu et al., 2016a;Chao et al., 2016;Song et al., 2017;Xu et al., 2018;Ho et al., 2019;Wu et al., 2021). These ingredients were numbered Z1-Z12 and B1-B7, and compound CID was confirmed with the Pubchem database ( Table 2).

Identification of ZXD and NAFLD Targets
These 19 active compounds were next introduced into the Swiss platform to identify putative macromolecular targets. We identified 325 ZX targets and 327 BZ targets, and 186 overlapping targets (Supplementary Table 1). In addition, the OMIM and GeneCards databases were used to identify 1575 NAFLD-related targets, of which 122 targets overlapped with our list of ZXD targets (Supplementary Table 2) (Figure 2). These overlapping targets were retained, and corresponding Uniprot IDs were determined and used to identify corresponding gene symbols for each target.
These 19 active ingredients and 122 targets were incorporated into a herbs-compounds-targets network wherein node size was scaled according to the degree values (Figure 3). The targets with the highest degree values included HMGCR, MAPK1, and SREBP-2, suggesting that these targets may play critical roles in the context of NAFLD treatment.

Functional Enrichment and PPI Network Analyses
Metascape platform was used to carry out KEGG pathway and GO enrichment analyses of the 122 overlapping ZXD and NAFLD targets. Results were ranked based on ascending q-values, corresponding to the degree of correlation with these targets, and the top 20 of these enriched pathways were plotted (Supplementary Table 3) (Figures 4-6). These targets were enriched for the lipid metabolism GO BP, suggesting that this pathway may be linked to the pharmacological activity of ZXD. Other identified targets associated with this pathway were also retained for further analysis (Supplementary Table 4).
To explore the interactions among the above targets in-depth, we constructed a PPI network using the STRING database with a medium confidence interaction threshold (0.4). Targets were then scored based on these predicted interactions and ranked in ascending order, retaining the targets that exhibited a combined score of ≥0.90 for further analysis (Supplementary Table 5). The degree values, centrality, and other parameters associated with these targets were then calculated using the Cytoscape Network Analyzer tool (Supplementary Table 6; Figures 9, 10). The data revealed precise, relevant interactions between these targets and highlighted the PI3K/AKT/NF-κB pathway as an essential target worthy of further investigation in the context of NAFLD treatment using ZXD. When all 122 overlapping targets were plotted, the resultant data also suggested that mitogen-activated protein kinases (MAPKs), fatty acid-binding proteins (FABPs), and matrix metalloproteinases (MMPs) may warrant future investigation in a therapeutic context.

Experimental Validation of Network Pharmacology
Significant differences (p < 0.01) in body weights were observed between the model and blank group rats that were fed HFD and NCD for 12 weeks, respectively ( Figure 11). Furthermore, the HE-stained liver tissues of the model group rats revealed significant differences compared with the tissues of the blank group rats (Figure 12 A-B; Figures 13A,B). Specifically, liver tissue sections from the blank group rats had intact, clear boundaries; The hepatocytes had a cord-like arrangement and were radially distributed around the central vein, and the cell nuclei was clearly visible. In contrast, swelling of hepatocyte, deposition of big lipid drops and hepatic ballooning degeneration were observed in liver tissue sections from the model group rats. The differences implied that the NAFLD modeling was successful. Afterwards, we treated NAFLD rats with the ethanol extract of ZXD across a range of concentrations (100, 200 and 400 mg/kg/d, respectively) for 4 weeks, and PH (10 mg/kg/d) was used as a positive control. We observed significant differences in the histopathological evaluation ( Figure 12 C-F; Figures 13C-F) and serum biochemical assay ( Figure 14) data of different groups. The curative effect of ZXD on NAFLD was found to be dose-dependent in rats at concentrations 100, 200 and 400 mg/kg/d. The liver tissues of blank, model and high-dose group rats were further conducted RT-qPCR. The results indicated that the gene expression of SREBP-2 and HMGCR of high-dose group rats decreased significantly (p < 0.05) in comparison with the blank group. The expression of MAPK1 and NF-κBp65 demonstrated an attenuation trend, which was however not significant ( Figure 15). Figures 16, 17 show the amplification and melting curves of these four genes (Figures 16, 17).

Identification of NAFLD-Related Targets of ZXD Treatment
ZXD is one of the TCMFs that exhibits anti-inflammatory and lipid-lowering properties. ZXD has been used to treat NAFLD because of hyperlipidemia-and inflammation-related etiology of this disease. However, the multi-component multi-target nature of ZXD and other TCMFs has made it difficult to firmly establish the molecular mechanisms underlying their pharmacological activity. In our previous study, we used LC-MS to conduct ZXD quality control analysis . This analysis, together with a review of relevant pharmacologically active compounds documented in the literature, revealed that alisol A (Ho et al., 2019), alisol A 24-acetate (Li et al., 2016), alisol A 23-acetate (Lin, 2012), 25-anhydroalisol F (Bi et al., 2017, alisol F (Bi et al., 2017), and alisol M 23-acetate (Xu et al., 2018) in ZX and atractylenolide I (Chao et al., 2016;Tang et al., 2017), atractylenolide II (Chao et al., 2016), and atractylenolide III (Song et al., 2017) in BZ are highly potent despite not meeting the predictive thresholds of the TCMSP platform (OB ≥ 30% and DL ≥ 0.18). As such, OB and DL alone are not sufficient to predict the active compounds within TCMFs.
The Swiss Target Platform was then used to identify putative NAFLD-related target genes. Further, we constructed a herbcompound-target network that identified HMGCR, MAPK1, and SREBP-2 as targets with high degree values that may be relevant in the context of ZXD-mediated NAFLD treatment.

Targets Selected for Validation
HMGCR is the rate-limiting enzyme that regulates cholesterol synthesis. Patients with NAFLD are reported to have high serum FIGURE 15 | Expression of SREBP-2, HMGCR, NF-κB and MAPK1 in SD rats. *p < 0.05 and **p < 0.01 vs. model group rats; #p < 0.05 and ##p < 0.01 vs. blank group rats. levels of total cholesterol (TC) (Puri et al., 2007;Caballero et al., 2009). The gene expression of HMGCR, however, showed an abnormally elevated serum level of TC. The expression of HMGCR is regulated by SREBP, especially by SREBP-2, a transcription factor primarily involved in the synthesis of cholesterol and fatty acid (Brown and GOLDSTEIN, 1997;Wu et al., 2013). Therefore, we hypothesized that HMGCR and SREBP-2 induce the onset and development of NAFLD by affecting the serum level of TC. There are reports that ZXD suppresses the expression of HMGCR and SREBP-2 (Zheng and WU, 2010;Dan et al., 2011). Further, the major pharmacological components of ZXD -alisol B 23-acetate and alisol A 24-acetatecould bind to HMGCR competitively (Xu et al., 2016b).
MAPK signaling plays a central role in many physiological and pathological processes associated with inflammation and NAFLD (Cao et al., 2017;Jiang et al., 2019). MAPK1, MAPK8, and MAPK9 were high-scoring targets in our analysis and may thus be associated with NAFLD treatment. The pharmacologically active components of BZ are reported to modulate MAPK signaling . For instance, atractylenolide I inhibits MAPK by suppressing TLR4/MyD88 pathway activity, while atractylenolide II and III suppress inflammation by affecting MAPK phosphorylation.
PI3K/AKT/NF-κB pathways were found to be closely associated with ZXD-mediated NAFLD treatment in PPI network analyses. Since many studies Fan et al., 2017) have confirmed the relationship between NAFLD and the PI3K/AKT/NF-κB pathway, NF-κBp65, the critical target of the PI3K/AKT/NF-κB signaling pathway, was selected to validate our findings. Furthermore, HMGCR, MAPK1, SREBP-2, and NF-κBp65 were selected for further validation through RT-qPCR.

RT-qPCR Experimental Validation Using in vivo NAFLD Model in SD Rats
The NAFLD model in rats was established using HFD and the success of modeling was verified by comparing changes in body weight and liver tissue section. The curative effect of ZXD was evaluated by histopathological evaluation and serum biochemical assays, and ZXD ameliorated lipid accumulation and treated NAFLD in a dose-dependent manner.
For validation purposes, RT-qPCR was performed to test the expression of HMGCR, MAPK1, SREBP-2, and NF-κBp65. ZXD significantly attenuated the expression of HMGCR and SREBP-2 (p < 0.05). The expression of MAPK1 and NF-κBp65 demonstrated an attenuation trend, which was however not significant (p > 0.05). The results indicated that our network pharmacology-based prediction model of molecular mechanism is reliable to some extent.
In future studies, we plan to perform a more in-depth investigation of the molecular mechanisms of ZXD based on the gut-liver axis to study the relationship between gut FIGURE 16 | The amplification curve of the four major genes (1, SREBP-2; 2, HMGCR; 3, MAPK1; 4, NF-κB).

CONCLUSION
In conclusion, we constructed a herb-compound-target-pathway network. Using biological enrichment analyses, we determined that ZXD can treat NAFLD in part by modulating the key targets -HMGCR, SREBP-2, MAPK1, and NF-κBp65. The experimental validation was performed through NAFLD model in SD rats. RT-qPCR analyses revealed that ZXD could significantly attenuated the expression of the four genes. The validation results indicated that our network pharmacology-based prediction model of molecular mechanism is reliable, which can be utilized for pilot mechanism investigation of TCMFs.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
YS provided the way of thinking. JW and FZ completed the design of framework of this manuscript. JW conducted the network pharmacology study. JW, FZ and HR completed the animal experiments. XC, JW and ZL investigated the background of the present review. XC, YS and WJ were responsible for supervision, review and editing. JW wrote the manuscript and FZ was responsible for visualization.