An Integrated Pharmacology-Based Analysis for Antidepressant Mechanism of Chinese Herbal Formula Xiao-Yao-San

Clinical studies and basic science experiments have widely demonstrated the antidepressant and anxiolytic effects of the herbal formula Xiao-Yao-San (XYS). However, the system mechanism of these effects has not been fully characterized. The present study conducted a comprehensive network pharmacological analysis of XYS and sorted all pharmacologically active components (149) through the TCMSP webserver. Then, all potential molecular targets (449) were predicted, of which there were 99 genes clearly related to depression. To further investigate the mechanism of antidepressant effects of XYS, a compound-depression targets (C-DTs) network was constructed, and Gene Ontology (GO) functional and KEGG pathway enrichment analyses were performed for the 99 targets. Enrichment results revealed that XYS could regulate multiple aspects of depression through these targets, related to metabolism, neuroendocrine function, and neuroimmunity. Prediction and analysis of protein–protein interactions resulted in selection of three hub genes (AKT1, TP53, and VEGFA). In addition, a total of seven ingredients from XYS could act on these hub genes and they were identified through ultra-high-performance liquid chromatography-quadrupole time-of-flight mass spectrometry (UPLC-Q/TOF-MS), including paeoniflorin, quercetin, luteolin, acacetin, aloe-emodin, Glyasperin C, kaempferol. Hereafter, we investigated the effects of paeoniflorin and its predicted target, the results suggest that it can reverse the neurotoxicity produced by CORT and could be a neuroprotective effect by promoting the phosphorylation of Akt. Overall, our research revealed the complicated antidepressant mechanism of XYS, and also provided a rational strategy for revealing the complex composition and function of Chinese herbal formula.


INTRODUCTION
Depression is a psychological disorder with widespread prevalence, affecting more than 350 million individuals worldwide, and the prevalence is gradually increasing (Kuo et al., 2015). Because of high morbidity and mortality, depression has received extensive attention. Currently, the main treatment strategy for depression is pharmacotherapy, but antidepressants typically do not completely alleviate depressive symptoms and may lead to drug addiction and adverse side effects for patients (Tom et al., 2014;Solem et al., 2017).
XYS has been widely prescribed as a safe and effective treatment or adjuvant therapy for depression, because psychological stress syndromes primarily are classified as LSSDS in TCM theory. Clinical studies and basic science experiments have demonstrated the antidepressant and anxiolytic effects of XYS (Dai et al., 2010;Chen and Hou, 2017), but the associated mechanisms have not been characterized.
Systems pharmacology is an emerging field that integrates bioinformatics and experimental methods to advance drug discovery and provides a method for clarification of mechanisms of action of Chinese herbs (Berger and Iyengar, 2009). A drug (active compounds)-targets-diseases (clinical symptoms) network can be constructed and analyzed using a holistic view though this method (Kim et al., 2019;Zhou et al., 2019). It can determine the enrichment of targets, and elucidate complex effects and pharmacological mechanisms of Chinese herbs (Kloft et al., 2016). The present study aimed to identify the bioactive components and mechanisms of action of the TCM formula XYS in the treatment of depression using system network pharmacological analysis.

Construction of XYS Chemical Constituent Databases
All chemical constituents of each of the herbs in XYS were obtained from the TCMSP data server (Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform, http://lsp.nwu.edu.cn/tcmsp.php), which is a unique systems-level pharmacology server for TCM (Ru et al., 2014). TCMSP can also provide pharmacology-related properties and predict targets and related diseases for each active ingredient based on critically evaluated pharmacological and clinical knowledge. This ground-breaking data server highlights the role that systems pharmacology can play in the traditional Chinese medicine discipline. Moreover, ADME properties data were also derived from TCMSP, comprising molecular weight (MW), octanol-water coefficient (ALogP), number of possible hydrogen-bond donors (nHdon), and number of hydrogen-bond acceptors (nHacc).

Evaluation the Drug-Likeness
Potential drugs in XYS were mainly identified by integrating oral bioavailability (OB) and drug-likeness (DL) properties. Comprehensive analysis of bioavailability and structural descriptors for predicting OB values for each compound were previously determined in-silico (Xue et al., 2012). [OB%]. Data server-dependent DL models were used to determine solubility and chemical stability based on the Tanimoto coefficient, using the following formula: T (A, B) = (A×B)/(|A| 2 + |B| 2 -A×B) (Tao et al., 2013). Then, according to recommendations from TCMSP, we selected screening criteria of OB ≥ 30% and DL ≥ 0.18 for determination of possible bioactive components.
Known Antidepressant Drug/Drug-Like Compounds Database Construction Known antidepressant drug/drug-like compounds from the cMap database Gene profiling of depression was performed using GSE12654 from GEO database, which contains frozen brain tissues from 11 depressed individuals and brain tissues from 15 nondepressed individuals (Iwamoto et al., 2004). Differentially expressed genes (DEGs) with the criteria |logFC|>1.2 and P < 0.05 were determined by "limma" package in R studio software. Moreover, the connectivity map (cMap) database is a collection of genome-wide transcriptional expression data from cultured human cells treated with bioactive small molecules that uses simple pattern-matching algorithms allowing for elucidation of functional connections between drugs, genes, and diseases through transitory common gene expression changes (Fang et al., 2003). After determining dysregulated genes, we converted the associated gene symbols to Affymetrix probe IDs, as the cMap database is based on these identifiers. P <0.05 was used as the threshold for statistical significance.

Known antidepressant drug/drug-like compounds from DrugBank database
DrugBank (http://www.drugbank.ca/) is a unique online database used as a bioinformatics and cheminformatics resource that compiles extensive drug-related information (Wishart et al., 2017). Therefore, DrugBank was also searched for antidepressant drugs or drug-like compounds which were approved by the United States Food and Drug Administration (FDA) for clinical trials.

Principal Component Analysis and Correlation Analysis
To evaluate the chemical distribution of antidepressant-related compounds in XYS, we used principal component analysis (PCA) for physicochemical-related parameters (MW, ALogP, nHDon, and nHAcc) though the prcomp function in R studio software. Therefore, known antidepressant drug/drug-like compounds from cMap, DrugBank, and active ingredients from XYS were also analyzed by PCA. Furthermore, the correlation of ADME properties of these compounds was also analyzed. Pearson r coefficients were calculated to evaluate the correlation between the three groups.

Target prediction
Corresponding protein targets of each active ingredient in XYS were predicted using several databases. Predicted targets from TCMSP were extracted and converted to target names using UniProtKB (http://www.uniprot.org) (Emmanuel et al., 2016). Then, PubChem, STITCH (Kuhn et al., 2012), and PharmMapper were used to obtain potential targets. The prediction results from the PharmMapper server (http://lilab.ecust.edu.cn/pharmmapper/) were based on the pharmacophore model, and this server is a large, in-house collection of pharmacophore databases extracted from all targets in TargetBank, DrugBank, BindingDB, and PDTD . Target information was set to Homo sapiens and aggregated for further analysis.

Construction of an Active Compound-Targets Network for XYS
To clarify the relationship between active compounds in XYS (C) and depression-related targets, depression-related genes were retrieved from the OMIM database (Ada et al., 2005), DisGeNET database (Piñero et al., 2015), TTD database (Zhu et al., 2010), and GSE12654 gene chip. Intersections of disease (depression)related genes (DGs) with the above predicted targets were determined to create an XYS antidepressant C-DTs network. This network was constructed and visualized using Cytoscape 3.4.0 software.

Functional enrichment analysis of XYS related targets
To analyse targets groups of active compounds in XYS, we used enrichment analysis, as this method can effectively increase the reliability of the identification of biological phenomena, resulting in meaningful annotation information. (i) GO annotation: In this study, "clusterProfiler" package in R studio was used for GO (Gene Ontology) enrichment analysis. "ClusterProfiler" is based on multiple resources, and also serves as a userfriendly enrichment tool with integrated gene cluster analysis (Yu et al., 2012). This R package provides for GO gene annotation enrichment analysis, including cellular component (CC), molecular function (MF), and biological process (BP) analyses. P < 0.05 were considered statistically significant. (ii) KEGG pathway analysis: To explore potential biological mechanisms of predicted targets, the "clusterProfiler" R package was used to annotate Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of targets genes potentially regulated by active compounds contained in XYS. In addition, significantly enriched KEGG pathways should contain at least three genes, and have P < 0.05.

Interactions between targets in the C-DTs network
Protein-protein interaction data analysis is an important mechanism to characterize the molecular basis of disease. Interactions between each target in the C-DTs network which could be regulated by XYS were searched using the STRING website, which is a comprehensive prediction server (Szklarczyk et al., 2017). Interaction scores should be more than 0.9 to indicate appropriate confidence in identified protein interactions. Using this tool, it is feasible to discover targets related to active components of XYS based on protein network.
To evaluate the antidepressant role of each active compound in XYS, hub DGs were identified using four different analysis methods in "cytoHubba" of cytoscape 3.4.0 software (Chin et al., 2014). The potential key active compounds in XYS for treatment of depression were incorporated into a Sankey diagram to reveal relationships between "hub genes -active compounds."
The information of equipment used in the UPLC-Q/TOF-MS analysis was detailed described in our previous research (He et al., 2018). The universal mixed mobile phase composed of water (solvent A) and acetonitrile (solvent B), both of which contain 0.1% formic acid (v/v). Gradient elution system used a flow rate of 0.4 mL/min, and adopted as follow: 5% B, 0-0.5 min; 5-80% B, 0.5-10 min; 80-100% B, 10-12 min; 100% B, 12-13 min; 100-5% B, 13-14 min; 5% B, 14-15 min. The injection volume was set as 2 µL and under 280 nm wavelength to detect. Mass spectral data were collected in Centroided. And the purity of all reference standards was above 98% that analyzed by HPLC (High performance liquid chromatography).

Cell viability assay
The effect of paeoniflorin reverse neurotoxicity from corticosterone, cell viability was determined by the MTT assay (Fotakis and Timbrell, 2006). Briefly, PC12 cells were seeded into 96-well plate with the concentration of 5 × 10 3 /100 µl, then incubated the cells under the conventional culture condition overnight. After 24 h treatment that with the mixture of 200µM CORT plus different concentrations of paeoniflorin (10 µM, 20 µM, 40 µM), 10 µl MTT (0.5 mg/ml) was added to each well, followed by 4 h of incubation at 37 • C. Then carefully remove the supernatant, and 100 µl DMSO was added to each well to dissolve the formazan crystal. The absorbance was recorded at 570 nm by microtiter plate reader. Cell viability was calculated as a percentage of the untreated or vehicle-treated control and averaged from three independent experiments.

Western blot analysis
Total protein was extracted from the PC-12 cells in different groups. Briefly, the collected were immediately lysed in RIPA lysis buffer containing PMSF and phosphatase protein inhibitor (Beyotime Biotechnology, Jiangsu, China). The protein lysates separated by 10% SDS-PAGE gel were electrophoretically transferred to polyvinylidene difluoride (PVDF) membranes (Millipore, Billerica, MA, USA). The following specific primary antibodies were used for incubate membranes overnight at 4 • C: anti-GAPDH (Beijing BioDee Biotechnology Co., Ltd.), anti-AKT (Cell Signaling Technology), anti-Phospho-AKT (Cell Signaling Technology). After probed with a 1:5000 dilution of horseradish peroxidase (HRP) conjugated secondary antibody (Beijing BioDee Biotechnology Co., Ltd.), the labeled proteins were detected by using enhanced chemiluminescence (ECL) detection.

Molecular docking evaluation
Furthermore, AutoDock 4.0 was used to evaluate the hub genes with active compounds, which is widely used for the study of bound conformation and binding free energy of ligand to the protein (Forli et al., 2016).

Statistical Analysis
Data are presented as means ± standard deviation. The standard two-tailed student's t-test was used for analysis and P < 0.05 were considered significant.

Bioactive Compounds in Xiaoyao San
Bioactive ingredients in XYS and corresponding ADME information were extracted from the TCMSP data server by filtering by drug-like threshold (OB > 30% and DL > 0.18). Several of the herbal medicines in XYS contain the same active ingredients, such as quercetin, a natural flavonoid, which is present in BR and in GRH. In addition, stigmasterol is present in AS, RB and ZRR. In total, 149 active compounds were identified in XYS, including 12 in PN, 2 in AS, 18 in BR, 15 in PR, 7 in AMR, 91 in GRH, 5 in ZRR, and 10 in MH. Detailed information regarding these molecules is summarized in Table S1.
To determine the molecular diversity of the constituents of each herbal medicine, the constituents were evaluated based on four significant physicochemical properties (Figures 1A-D

Identified Antidepressant Drugs and PCA
We retrieved the publicly available microarray dataset GSE12654 from GEO, which included frozen brain tissues of 11 humans with depression and 15 individuals without depression (Iwamoto et al., 2004). After applying thresholds log|FC| >1.2 and P < 0.05 using "limma" in R studio, a total of 29 significantly differentially expressed genes between depressed and normal patients were identified. Fourteen genes were down-regulated, and 15 genes were up-regulated (Table S2).
Furthermore, to evaluate the relationship between ingredients form XYS and other known-compounds, we first identified the known-compounds which contain anti-depression function based on the DEGs results from cMap data server. Finally, 12 significant compounds with adequate relative connection mean scores (absolute value >0.5) were identified. Of the twelve compounds, 2 had high positive connection scores, and the other ten compounds possess negative connectivity scores, indicating they can act with these dysregulated genes in depression. Figure 2A and Table 2 summarize these results. In addition, a search of the DrugBank database resulted in identification of 64 FDA-approved compounds for treatment of depression ( Table 2).
PCA was performed to characterize the physicochemical parameters (MW, ALogP, nHDon, and nHAcc) of these antidepressant drugs/compounds. PC1 and PC2, as shown in Figure 2B and Table S3, showed a distribution of different color nodes that flocked together, indicating extensive overlap among the three groups of ingredients in XYS, and the antidepressant compounds determined by cMap and Drugbank. Furthermore,  the significant positive correlation between these compounds of three groups ( Figure 2C, Table S3). In ADME properties, Pearson r value is 0.3404 (P < 0.0001) for MW, AlogP (r = 01969, P = 0.003), Hdon (r = 0.3018, P < 0.0001), Hacc (r = 0.351, P < 0.0001). These results indicated that bioactive compounds from XYS are highly correlated with known antidepressant compounds in molecular properties, and could as potential drug therapies for depression treatment.

Targets of Xiaoyao San
An in silico target screening approach was performed to identify the putative targets of each bioactive component of XYS. Predictive tools were used, including the TCMSP server, STITCH, and PharmMapper data servers. In total, 446 targets were identified for 149 bioactive compounds (9 compounds did not receive predictions), and many compounds were predicted to target the same proteins, such as Monoamine Oxidase B (gene symbol: MAOB). Depression-related putative targets were obtained from Therapeutic Target, DrugBank, OMIM, and DisGeNET databases, and DEGs resulting from analysis of GSE12654. Finally, 121 bioactive compounds in XYS which associated with 99 depression-related targets were identified. Multiple therapeutic targets were mediated by active ingredients of XYS, such as IL2, IL4, IL6, IL10, and STAT3, which are involved in immune and inflammatory responses closely associated with depression. For example, STAT3 (Signal transducer and activator of transcription 3), a component of the JAK/STAT signaling pathway, influences the fate and function of brain cells.

Construction of a Compound-Depression Target Network (C-DT)
A C-DT network was constructed (Figure 3) for further biological analysis. This network contained information regarding complex relationships between all bioactive compounds in XYS and their depression-related targets. There were 220 nodes in the C-DT network composed of 121 active ingredients and 99 depression targets, connected by 1583 interactions, the details were listed in Table S4. As shown in Figure 3, node size indicated the number of connections calculated using "NetworkAnalyzer, " an analytical tool in cytoscape 3.7.0.

GO and KEGG Pathway Analysis
To investigate biological functions associated with this interaction network relative to treatment of depression, we performed GO and KEGG pathway enrichment analyses using the "clusterprofile" R package. Three GO categories were analyzed, which included CC, MF, and BP. The top 15 GO enrichment results for each category are listed in Table 3. Direct outcomes are shown in Figure 4A. These targets were significantly involved in formation of neuronal receptor and cell, which are directly related to generation of brain cells and tissue. Molecular functions involved activity and binding of various molecules, neurotransmitter receptors, cytokines, transmembrane receptor protein kinases, and nuclear receptor, as well as cytokine receptor binding, phoshphatase, and steroid. Moreover, combined BP results indicated that XYS may regulate depression-related biological processes, such as "regulation of cell migration, " "reactive oxygen species metabolic process, " and "aging" to exert antidepressant effects. KEGG signaling pathway results are shown in Figure 4 and Table 4. The top enriched terms were related to several critical pathways associated with disease, such as "PI3K-Akt signaling pathway, " "MAPK signaling pathway, " and "neuroactive ligandreceptor interaction" (Figure 4B and Figure S1). These findings demonstrated complex pathological features of depression and the potential of XYS as an antidepressant.

PPI Network Analysis and Screening of Key Active Ingredients
To investigate the characteristics of protein targets in the C-DT network, protein-protein interactions were analyzed using the STRING online database, and the related PPI network was constructed and visualized using cytoscape software. Using confidence of 0.9 as the threshold, we found 90 nodes and 288 interactions in the PPI network. In addition, we used four different analysis methods (Stress, Closeness, MNC and EcCentricity) using the "cytoHubba" plugin to identify the top 10 hub nodes in the PPI network. Detailed results are shown in Figures 5A-D, and the top hub genes based on these four indices were summarized in a Venn diagram ( Figure 5E, Table S5). Thus, 3 core genes were identified as hubs: AKT1, TP53, and VEGFA. Furthermore, based on GO and KEGG analyses, these three core genes were enriched in multiple biological processes and signaling pathways. For instance, AKT1 is involved in myelination of the peripheral nervous system, peripheral nervous system axon ensheathment, and the neurotrophin signaling pathway, all of which are closely related to depression (Barros et al., 2009;Musashe et al., 2016;Liu et al., 2018). Using these core genes, a Sankey diagram (Figure 6A) was used to facilitate shown of bioactive compounds that can act on multiple targets. Nine ingredients of XYS and their hub targets were identified,  including quercetin, luteolin, acacetin, aloe-emodin, gadelaidic acid, Glyasperin C, kaempferol, naringenin, paeoniflorin.

Identification of the Chemical Characterization of XYS
As shown in Figure 6B, the chemical base peak ion (BPI) chromatogram of XYS based on the positive and negative ion modes of UPLC-Q-TOF/MS. According to the above analysis, nine key compounds of XYS were further examined. All of the key compounds from XYS, seven ingredients (paeoniflorin, quercetin, luteolin, acacetin, aloe-emodin, Glyasperin C, kaempferol) were successful identified by UPLC-Q-TOF/MS analysis, the details were listed in Table 5. Among them, we found major compounds are flavonoids, only paeoniflorin as monoterpene glycoside, and most importantly it is a principal ingredient of Paeonia lactiflora Pallas ( Figure 7A). Thus, we next performed cellular model and investigated the neuron protected effect of paeoniflorin in vitro.

Evaluations
To determine the protected effect of Paeoniflorin from CORT neurotoxicity, the cell viability assay was examined. The detailed MTT assay results as shown in Figure 7B, that 200 µM CORT reduced the cell viability of PC12 cells, whereas treatment 20 and 40 µM paeoniflorin were very significant able to moderate the cell viability (P < 0.01). Western blot analysis results showed that compared with the control cells, the expression of p-Akt in CORT stimulated cells was significantly decreased (P < 0.01). However, the down-regulated expression of p-Akt was significantly reversed by paeoniflorin treatment (Figures 7C,D). Furthermore, quercetin and luteolin could act on all hub genes, so we further validated their molecular docking to detect and calculate the binding energy between small molecules and proteins. The best bound conformation site of quercetin and luteolin binding to proteins were shown in Figure 8, respectively. The binding energy between quercetin with AKT1 is −17.6341 kcal/mol, and luteolin with AKT1 is −12.7977 kcal/mol, that initiated both two compounds had strong binding activity with AKT1. Similarly, these two compounds contain strong binding activity with other hub targets, the binding energy between quercetin with VEGFA and TP53 (−67.4121 kcal/mol), luteolin with VEGFA and TP53 (−53.3754 kcal/mol).

DISCUSSION
Traditional Chinese medicine (TCM) is being more frequently chosen by patients to treat depressive and anxiety disorders. A meta-analysis reported that use of XYS as an adjuvant therapy with antidepressant drugs in 1,837 patients was more effective in improving symptoms than antidepressants alone (Zhang et al., 2011). Furthermore, previous studies indicated that XYS treatment of depression and anxiety involves several biological aspects. For instance, some psychological disorders have been associated with increased cytokine levels resulting from immune activation, known as the "cytokine hypothesis" (Sharma, 2016). Inflammatory cytokines may act as modulators of 5-HT (serotonin) receptors, and upregulate indoleamine 2,3dioxygenase (IDO), resulting in anxiety-like disorders (Myint and Kim, 2003). The antidepressant-like effects of XYS were associated with downregulation of peripheral IL-1β, IFNγ, TNF-α, and IL-6, and increased synthesis of tryptophan hydroxylase (TPH) and IDO, which play a role in promoting 5-HT synthesis (Jiao et al., 2019). XYS significantly inhibited TNF-α release in serum and the hippocampus through activation of the TNF-α/JAK2-STAT3 pathway, as shown in our preliminary study . Moreover, the HPA axis and brain-derived neurotrophic factor (BDNF) have also been widely studied in the context of pathogenesis and treatment of depressive disorders (Jiang et al., 2016). Previous studies have demonstrated that XYS could reverse chronic stress-induced anxiety-like symptoms by regulating the apelin-APJ system in the hypothalamus and activity of the HPA axis (Yan et al., 2018).
Many studies have demonstrated the antidepressant effects of XYS, which involves multiple pathways associated with depression. Therefore, the complex bioactive ingredients and targets with potential antidepressant activities require clarification. Previous studies using High Performance Liquid Chromatography coupled with LTQ Orbitrap Mass Spectrometry (HPLC-LTQ-Orbitrap-MS) identified paeoniflorin, liquiritin, glycyrrhizic acid, ferulic acid, saikosaponins A and C, curcumin, and poria cocos alcohol were eight active compounds identified in XYS (Li et al., 2015). Yuzhi Zhou and his colleagues identified 4 components: Z-ligustilide, palmitic acid, atractylenolide I, and atractylenolide II in XYS using gas chromatographymass spectrometry (GC-MS) methods (Zhou et al., 2012). In addition, network pharmacology methods have been used to predict and analyse 13 components of XYS (Gao et al., 2015). However, no comprehensive analysis of all bioactive compounds in XYS had been previously reported. Therefore, the present study adopted a comprehensive systematic network pharmacology approach for evaluation of XYS.
In this study, we identified 149 compounds in active fractions of XYS. PCA analysis and correlation analysis was used to compare these compounds with antidepressants approved by the FDA and small molecule drugs predicted by the cMap database. This was the first comparative study of drug properties of TCM formula compounds using Drugbank and a gene chip combined with cMap. Our results indicated that these 149 bioactive compounds in XYS possessed potential antidepressant properties. Therefore, we constructed a compound-depression targets (C-DTs) network, and further enriched this network using GO KEGG pathways analyses of 99 depression-related targets. CC enrichment results indicated that a considerable number of genes were involved in formation of nerves and synapses. Synapses are the fundamental structures for information transmission between neurons, and they adapt to stimuli by continuously modifying neural connections and neurological circuits. Therefore, synaptic plasticity is the main manifestation of neuroplasticity, and depression and other psychological disorders are typically associated with decreased synaptic plasticity in the hippocampus (Wainwright and Galea, 2013). Molecular functions associated with components of XYS involved activity and binding of various molecules, such as neurotransmitter receptors, cytokines, transmembrane receptor protein kinases, nuclear receptors, cytokine receptor binding, steroids, and catecholamines. Moreover, a large number of genes were significantly enriched as BP terms in regulation of aging, reactive oxygen species metabolic processes, and inflammatory response, which are associated with pathogenesis of depression. Combined BP results indicated that XYS may act via generation of nerves and synapses, alteration of neurotransmitter receptor function, and regulation of cytokines, resulting in antidepressant effects. The PI3K/Akt signaling pathway was enriched to the greatest degree, resulting in association of 20 genes. It has been wildly regarded to modulate antidepressant-like functions and contributes to synaptic plasticity and neurotransmission formation (Ludka et al., 2016). Furthermore, several metabolic pathways directly related to the nervous system such as "neuroactive ligand-receptor interaction, " "GABAergic synapse, " and "neurotrophin signaling pathway" were also enriched. Analysis of the PPI network identified hub nodes which were most closely related to other proteins: AKT1, TP53, VEGFA. These proteins were consistent with our enrichment analysis, and have been previously reported to be related to neurotrophic factors and depression.  Bioactive ingredients in XYS were classified in the following categories: volatile oils, alkaloids, flavonoids, saponins, and polysaccharides. Among these, polysaccharides can regulate intestinal flora to moderate depression though the brain-gutmicrobiota axis (Cryan and Dinan, 2012;Schwalm et al., 2016), and flavonoids can cross the blood-brain barrier, directly acting on brain targets (Youdim et al., 2003). We concluded that XYS exerts its effects in many ways through the actions of many compounds and may confer advantages over the compounds by reducing drug resistance.
After identified hub targets, 9 ingredients of XYS could interact with them, and 7 of them were identified by UPLC-Q/TOF-MS analysis. As a result, we found that 3 hub genes predicted by the PPI network could be regulated by quercetin and luteolin, which are both flavonoids. Nevertheless, most of the reports of flavonoids are limited to in vitro experiments, and these flavonoids are common constituents in a variety of herbs. It is considered important to clarify their specific role in the treatment (Heinricha et al., 2020). Actually, many studies have verified their neuroprotective effects in vitro and in vivo. For instance, the genes regulated by luteolin could act as receptors in the central nervous system, resulting in antidepressant effects (Sasaki et al., 2013;Akinrinde and Adebiyi, 2019). Quercetin, a flavonoid, is a component of Bupleurum Chinense Dc. This flavonoid confers antioxidant activity and antiinflammatory effects in the treatment of neurological diseases. A previous study reported that quercetin relieved anxiety and depression of chronic stress rats by protecting neurons from oxidation and inflammation (Mehta et al., 2017). In LPSinduced neuroinflammation of mice, quercetin still exhibited a function of eliminating inflammation through significantly reduced LPS-induced proliferation of astrocytes in the brain and decreased the expression of inflammatory factors (Khan et al., 2018). Additionally, researchers found that quercetin could inhibit neuronal apoptosis to preventing anxiety-like behaviors in vivo via regulating the Akt1 and ASK1/JNK3/caspase-3 expression (Pei et al., 2016). Besides, in a previous network pharmacology study, quercetin was found also as a high contribution compound of Huangqi and Huanglian for treating diabetes, participated in polypharmacological and synergistic mechanisms (Yue et al., 2017). However, among these 7 ingredients, only paeoniflorin, a monoterpene glycoside compound, as main component of herb (Paeonia lactiflora Pallas) from XYS. And it is one of the quality identification standards for XYS in the Chinese Pharmacopeia. Thus, it should be a representative compound from this formula, so that the current research further analyzed its neuroprotective effect and validated its potential roles with predicted target Akt. Previous reports by several researchers had provided clues supporting the role for paeoniflorin in the protection of neuron cells and the therapeutic effects in neurological diseases (Mao et al., 2012;Cong et al., 2019). Gu et al. had demonstrated that paeoniflorin via upregulating the p-Akt expression and Bcl-2/Bax ratio to reduce neuron death in Alzheimer's disease mice (Gu et al., 2016). Recent evidence suggests that CUMS rats had improved depressive-like behavior after paeoniflorin treatment, it could be regulated by acting on the ERK-CREB signaling pathway (Zhong et al., 2018). Consist with previous findings, we found paeoniflorin can reverse the neurotoxicity produced through CORT and promote the phosphorylation of Akt. On the one hand, from the perspective of network pharmacology, the representative importance of paeoniflorin in XYS has been clarified for the first time. On the other hand, our results indicated that paeoniflorin as the XYS's representative compound which could modulate the several signaling upstream target AKT's expression, that validate the current network pharmacology prediction.

LIMITATIONS AND CONCLUSIONS
Current research based on network pharmacology which could provide a novel and systematic analysis way for the research of Chinese herbal formula. However, there are still some limitations in our article. Firstly, as we all know, components of Chinese herbal are complicated, it's doubtful whether there are synergistic or side effects caused by the interaction between compounds that is hard to clarified in the prediction and composition identification of this study. Secondly, although CORT is one of the most common methods to simulate depression in vitro, the pathogenesis of depression is complicated, which means that it is difficult to simulate whole pathogenic factors in vitro. Moreover, western blots suggested that paeoniflorin could regulate the phosphorylation of Akt, but it is only in protein expression level. The direct or indirect regulation between them would need to be further explored. Last but not least, hub genes were selected by PPI network in our study that they are usually as upstream molecules or central molecules in those important signaling pathways. However, there still exist multiple downstream proteins in pathogenesis of depression. And thus, it is important to discover and identify more compounds that interact with these downstream proteins through future works.
Overall, this study systematically summarized key active compounds in XYS and comprehensively analyzed potential mechanisms of its action and pathways. We also evaluated  for representative compound paeoniflorin in vitro. The current findings provide a basis and new insights for XYS in the treatment of depression. Certainly, we will explore widely experiments of pharmaceutical ingredients to better understand the effects of XYS, such as component analyses, identification of drugs, and metabolic studies in animal models. Densitometry data analysis were quantified using Image J software. Compare with control group, **P <0.01; compare with CORT group, # P < 0.05, ## P<0.01; n = 3. PF, paeoniflorin; CORT, corticosterone.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE12654.

AUTHOR CONTRIBUTIONS
NY and JC main contributed research conception and design. NY, KT, LH, WH, XL, and QM collected and processed data. NY, LG, and KT wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the final submitted version. JC takes primary responsibility for communication with the journal and editorial office during the submission process, throughout peer review, and during publication.

FUNDING
This work was financially supported by the National Natural Science Foundation of China (Nos. 81630104, 81803998, 81973748).

ACKNOWLEDGMENTS
The authors gratefully appreciate the reviewers and editors for their constructive suggestions and kind help. And we sincerely appreciate Dr. Huang Chuiguo for his warm help.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar. 2020.00284/full#supplementary-material Figure S1 | Distribution of partial targets of XYS on the Neuroactive ligand-receptor interaction signaling pathway. The red nodes are potential targets. The green nodes are relevant targets in the pathway.