An Advanced Systems Pharmacology Strategy Reveals AKR1B1, MMP2, PTGER3 as Key Genes in the Competing Endogenous RNA Network of Compound Kushen Injection Treating Gastric Carcinoma by Integrated Bioinformatics and Experimental Verification

Gastric carcinoma (GC) is a severe tumor of the digestive tract with high morbidity and mortality and poor prognosis, for which novel treatment options are urgently needed. Compound Kushen injection (CKI), a classical injection of Chinese medicine, has been widely used to treat various tumors in clinical practice for decades. In recent years, a growing number of studies have confirmed that CKI has a beneficial therapeutic effect on GC, However, there are few reports on the potential molecular mechanism of action. Here, using systems pharmacology combined with proteomics analysis as a core concept, we identified the ceRNA network, key targets and signaling pathways regulated by CKI in the treatment of GC. To further explore the role of these key targets in the development of GC, we performed a meta-analysis to compare the expression differences between GC and normal gastric mucosa tissues. Functional enrichment analysis was further used to understand the biological pathways significantly regulated by the key genes. In addition, we determined the significance of the key genes in the prognosis of GC by survival analysis and immune infiltration analysis. Finally, molecular docking simulation was performed to verify the combination of CKI components and key targets. The anti-gastric cancer effect of CKI and its key targets was verified by in vivo and in vitro experiments. The analysis of ceRNA network of CKI on GC revealed that the potential molecular mechanism of CKI can regulate PI3K/AKT and Toll-like receptor signaling pathways by interfering with hub genes such as AKR1B1, MMP2 and PTGERR3. In conclusion, this study not only partially highlighted the molecular mechanism of CKI in GC therapy but also provided a novel and advanced systems pharmacology strategy to explore the mechanisms of traditional Chinese medicine formulations.

Gastric carcinoma (GC) is a severe tumor of the digestive tract with high morbidity and mortality and poor prognosis, for which novel treatment options are urgently needed. Compound Kushen injection (CKI), a classical injection of Chinese medicine, has been widely used to treat various tumors in clinical practice for decades. In recent years, a growing number of studies have confirmed that CKI has a beneficial therapeutic effect on GC, However, there are few reports on the potential molecular mechanism of action. Here, using systems pharmacology combined with proteomics analysis as a core concept, we identified the ceRNA network, key targets and signaling pathways regulated by CKI in the treatment of GC. To further explore the role of these key targets in the development of GC, we performed a meta-analysis to compare the expression differences between GC and normal gastric mucosa tissues. Functional enrichment analysis was further used to understand the biological pathways significantly regulated by the key genes. In addition, we determined the significance of the key genes in the prognosis of GC by survival analysis and immune infiltration analysis. Finally, molecular docking simulation was performed to verify the combination of CKI components and key targets. The anti-gastric cancer effect of CKI and its key targets was verified by in vivo and in vitro experiments. The analysis of ceRNA network of CKI on GC revealed that the potential molecular mechanism of CKI can regulate PI3K/AKT and Toll-like receptor signaling pathways by interfering with hub genes such INTRODUCTION Gastric carcinoma (GC) is a malignant tumor of the digestive system that poses a serious threat to human health (Ajani et al., 2017). According to the statistics of International Agency for Research on Cancer, there were about 1.08 million new cases of gastric cancer worldwide in 2020, and about 0.77 million of people die from gastric cancer (Bray et al., 2018;Sung et al., 2021). The morbidity and mortality of GC have declined sharply in recent decades in some Western countries, while it is still relatively high in East Asia and represents a significant medical burden (Ferro et al., 2014;Torre et al., 2015). Among the factors that increase the risk of human GC, Helicobacter pylori gastric infection plays a vital role, and 75% of GC cases worldwide are caused by Helicobacter pylori infection (Plummer et al., 2015). Although early stage GC is highly treatable, the median survival time of advanced GC is only 9-10 months. Unsatisfactorily, the global 5-year survival rate of patients is still less than 30% (Verdecchia et al., 2007). The combination of different forms and different drugs of chemotherapy, radiotherapy and surgery are common treatment methods to treat GC (Coccolini, 2016). However, because of the internal metastasis and changes of the tumor, the heterogeneity of different patients and the side effects of radiotherapy and chemotherapy, patients' options in clinical practice are minimal (Ajani et al., 2017).
Compound Kushen injection (CKI), also called Yanshu injection, consists of Kushen (Radix Sophorae flavescentis) and Baituling (Rhizoma Smilacisglabrae) (Zhao et al., 2014). CKI has been adopted clinically for a decade to treat various solid tumors, including GC, liver cancer, lung cancer, breast cancer and other cancers (Xu et al., 2011;Sun et al., 2012;Wang et al., 2015). It is worth mentioning that CKI can also relieve cancer pain, regulate immunity, and improve conventional chemotherapy to reduce tumor efficacy and chemotherapy toxicity (Tu et al., 2016;Yang et al., 2020). The anti-tumor effect of CKI has been confirmed while the underlying molecular mechanism is still poorly understood.
Molecular studies have provided a vast quantity of new information for potential mechanisms to use in cancer treatment. Microarray and high-throughput sequencing technologies provide a reliable guarantee for the deciphering of significant genetic or epigenetic alterations in carcinogenesis and the discovery of potential biomarkers for cancer diagnosis, treatment, and prognosis (Cancer Genome Atlas Research Network, 2014). MicroRNA (miRNA) and long noncoding RNA (lncRNA) are the two most common subtypes of noncoding RNA (ncRNA). Their abnormality leads to the inability of mRNA to be transcribed normally and may contribute to unhindered growth, and invasion of cancer cells (Schmitt and Chang, 2016;Rupaimoole and Slack, 2017). At present, studies have confirmed that the lncRNA-miRNA-mRNA network plays a vital role in the occurrence and development of cancer, which may have substantial clinical prospects for identifying potential biomarkers and therapeutic targets of various tumors (Liu et al., , 2020. Thus, in this study, we firstly analyzed the microarray dataset in the Gene Expression Omnibus (GEO) database and The Cancer Genome Atlas (TCGA) to find miRNAs that are differentially expressed in GC compared to normal tissues. Secondly, we applied weighted gene co-expression network analysis (WGCNA) and merged with differentially expressed miRNAs (DEMs) predicted targets to identify genes associated with GC progression systematically. Afterward, we undertake a systematic study of the molecular mechanism of CKI in the treatment of GC using a network pharmacology analytical model and proteomics analysis. Drawing on the above research, we conducted a meta-analysis of key targets to verify their expression changes in GC and conduct immune infiltration to explore their prognostic impact on GC patients. At last, to better analyze and predict the molecular mechanism of CKI on GC, enrichment analysis, molecular docking and biological experiments were exploited to discover and verify the involved pathways and the binding of CKI components to key targets. Figure 1 depicts a workflow of the advanced systems pharmacology strategy used in this study. a Waters HSS T3 UPLC C18 analytical column (1.7 µm, 2.1 mm × 100 mm, Milford, MA, United States). The oven was set at 40 • C; the injection volume was 5 µL; the flow rate was set at 0.3 mL/min; the mobile phase was consisted of 0.1% formic acid in water (A) and carbinol (B). The eluting program was: 5-5% B for 0-1 min, 5-95%B for 1-12 min, 95-95%B for 12-15 min. The ion source was electrospray ionization (ESI); MS was operated in positive mode; the scan mode was Full scan/ddMS2, the scan range was 100-1200 Da; The capillary temperature was 350 • C; The spray voltage in negative mode was 3800 V; The spray voltage in positive mode was 3200 V; The sheath gas was 35 arb; The aux gas was 15 arb; three collision energies of low, medium and high were used for MS2. The positive ion mode was 30, 40, 50 eV, and the negative ion mode was 30, 50, 70 eV. The resolution of the full scan was 70000 FWHM, and the resolution of MS2 was 17500 FWHM. The reference marker compounds present in the sample were identified based on retention time, MS fragmentation and UV spectra.

Construction of Compound Kushen Injection Ingredient Prediction Target Network
For a more profound and comprehensive study, we input the 3D structure of chemical composition in Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) (Ru et al., 2014), Search Tool for Interactions of Chemicals (STITCH) (Szklarczyk et al., 2016), SuperPred (Nickel et al., 2014), SwissTargetPrediction (Gfeller et al., 2014) for target prediction. Moreover, the predicted multiple target information of the compounds and the obtained information was introduced into Cytoscape 3.6.1 (Franz et al., 2016) 1 to obtain an intermolecular interaction network and carry out complex network analyses.

Differential Expression of MicroRNAs in Gastric Carcinoma
Microarray data for gene expression GSE23739 were downloaded from the GEO database. A total of 80 samples were obtained, including 40 primary tumors and 40 normal samples. After the raw data has undergone background correction and standardization, the Limma (Ritchie et al., 2015) R package was applied to analyze the difference between cancer and normal tissues. The miRNA-seq data in TCGA contains 446 tumor samples and 45 normal samples. To verify and obtain DEMs, the edgeR (Robinson et al., 2010) package was used to analyze the difference between the groups.

Differentially Expressed MicroRNAs Target Genes Prediction
MiRWalk2.0 (Dweep and Gretz, 2015) is a comprehensive archive that fully integrates the interactions of multiple existing miRNA target prediction databases and provides predictive and experimentally verified miRNA target prediction. On the one hand, the interactions between miRNAs and genes from 12 servers were speculated and only the genes predicted by more than six of the servers were identified as target genes. On the other hand, 5 servers using miRWalk, miRanda, PITA, RNAhybrid and Targetscan were utilized to predict miRNA-lncRNA targets.
Weighted Gene Co-expression Network Analysis for Gastric Carcinoma mRNA

Data Collection and Preprocessing
The TCGA-STAD RNA-seq data includes 407 samples of its HTSeq-Counts data and associated clinical information downloaded in February 2020. After removing samples that contained incomplete analytical data and/or other malignancies, 375 samples were retained. Because, some genes without significant changes in expression between samples, we selected the top 5000 genes that are most relevant to differential expression for the following WGCNA analysis.

Weighted Gene Co-expression Network Analysis and Module Preservation
The gene co-expression networks were constructed using the WGCNA package (Langfelder and Horvath, 2008). The 1 http://www.cytoscape.org/ similarity between gene expression profiles was used to construct a similarity matrix based on pairwise Pearson correlation coefficient matrices. To improve co-expression similarity and achieve a scale-free topology, an appropriate soft threshold power β was selected by using the integration function (pickSoftThresshold) in the WGCNA software package (Jeong et al., 2001). We also reconstructed the topological overlap matrix by calculating the Topological Overlap Measure (TOM) which is a robust measure of network interconnectedness (Li and Horvath, 2007;Yip and Horvath, 2007). Finally, the Dynamic Tree-Cut algorithm method was applied to identify the module of gene co-expression with maxBlockSize of 6000, minModuleSize of 30 and mergeCutHeight of 0.2.

Identification of Clinical Significant Modules
The module eigengene (ME) is the first principal component of each gene module and the expression of ME is considered representative of all genes in one module. The Module Membership (MM) is the correlation between the ME and the gene expression profile. Gene Significance (GS) is the absolute value of the correlation between a specific gene and a clinical trait. According to ME, GS, MM, we can associate modules with clinical traits, not only to calculate the correlation between ME and clinical traits, but also to analyze clinically vital modules (Langfelder and Horvath, 2008).

Prediction of Competing Endogenous RNA Network of Compound Kushen Injection Intervention in Gastric Carcinoma
To systematically describe GC-associated underlying molecular mechanism, a competing endogenous RNA (ceRNA) network was conducted by merging the predictive correlation of DEMs and key modules in WGCNA. The target network predicted for CKI active component was combined with the ceRNA network of GC for CKI intervention in GC ceRNA network prediction, and the overlapping proteins in the two networks are likely to be the potential key gene for GC treatment by CKI active ingredients. Cytoscape constructed a potential ceRNA network for CKI treatment of GC, and the potential targets for CKI treatment of gastric cancer were systematically analyzed.

Functional Enrichment Analysis
To analyze the enrichment of key proteins, we first used the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) 10.5 2 database to construct a protein-protein interaction (PPI) network for key proteins (Szklarczyk et al., 2017). The STRING database aims to collect, score, and integrate all publicly available sources of protein-protein interaction information, and to complement these with computational predictions. Its goal is to achieve a comprehensive and objective global network, including direct (physical) as well as indirect (functional) -interactions. We performed the Gene Ontology (GO) Functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the predicted key targets of CKI compounds used in GC therapy to identify their biological functions. In addition, the R package clusterProfiler was used to perform GO and KEGG functional enrichment analysis. Particularly, the function and pathway enrichment analyses of the validated target genes of miRNAs, were used by the DIANA tool, which is based on the collaboration of the previously mentioned database (TarBase v7.0) and mirPathv3.0 (a miRNA pathway analysis web server deciphering miRNA function with experimental support) (Paraskevopoulou et al., 2016;Vlachos and Hatzigeorgiou, 2017).

Data Collection
A microarray search for hub genes was conducted in the GEO database with the following terms: ("stomach neoplasms"[MeSH Terms] OR gastric cancer[All Fields]) AND "Homo sapiens"[porgn] AND ("gse"[Filter] AND "Expression profiling by array" [Filter]) and the latest searching time was April 5, 2020. The criteria for inclusion were as follows: (1) patients diagnosed with gastric cancer were investigated; (2) cancerous and noncancerous samples were involved; (3) datasets samples were no less than 20. In addition, the following conditions caused the exclusion of a study: (1) lack of original data; (2) the patients with gastric cancer were accompanied by other tumors (3) the interventions included surgery, radiotherapy or other cancer treatments.

Statistical Analysis and Comprehensive Meta-Analysis
The expression profiling information of the datasets were exploited to calculate the mean (M) and standard deviation (SD) for each hub gene in the experimental control group. Then, the meta package of R software was brought into play the standardized meta-difference (SMD) and 95% confidential interval (CI) analysis. In addition to determining a reasonable choice of random effects and fixed effects models and evaluate heterogeneity, the chi-squared test of Q and the I2 statistic were calculated.

Survival Analysis of Hub Genes
The correlation between hub gene expression and overall survival was assessed using the Kaplan-Meier estimation method, based on the "survival" package in R. A significant difference in survival curves was assessed using a log-rank test. P value less than 0.05 was considered statistically significant.

Immune Infiltrates Analysis
TIMER 3 is a database that can comprehensively study the molecular characterization of tumor-immunity interactions. Not only can the association between immune infiltrates and a variety of factors be explored interactively, but also the dynamic analysis 3 https://cistrome.shinyapps.io/timer/ and visualization of these associations can be performed using a TIMER. In this study, we evaluated the hub gene expression in GC and its correlation with the abundance of tumor-infiltrating immune cells, via gene modules .

Molecular Docking
Molecular docking can reflect the binding energetics of drug molecules to protein receptors by calculating the binding affinity between ligands and receptors and the corresponding intermolecular interactions (Huang and Zou, 2010;Ferreira et al., 2015). The crystal structure of the key gene was downloaded from Protein Data Bank (PDB) 4 database. The 3D protein crystal structure had to be determined by X-ray crystallography and the crystal resolution was less than 3 Å. The protein receptor and ligand files were pre-processed by AutoDock Tools and then Autodock vina was used for molecular docking (Trott and Olson, 2010;Forli et al., 2016). In addition, Pymol and Ligplot were used to visualize the results to show the intermolecular interaction and docking more clearly (Laskowski and Swindells, 2011;Mooers, 2016).

Assessment of the Anticancer Effect of Compound Kushen Injection in vivo
The Institute of Hydrobiology, Chinese Academy of Sciences, provided normal AB wild-type zebrafish. The zebrafish were raised and bred under the conditions recommended in The Zebrafish Book: A guide for the laboratory use of zebrafish (Danio rerio). The experiment was in accordance with the Animal Management Rules of the Ministry of Science and Technology of the People's Republic of China for experimental care and use of animals and approved by the Animal Ethics Committee of Beijing University of Traditional Chinese Medicine.
According to the 1:1 ratio of males to females, healthy adult zebrafish were arranged for natural spawning. Then, zebrafish embryos were collected and cultured in a constant temperature incubator at 28.5 • C for the construction of the zebrafish xenograft GC model. Before tumor cell injection, the GC cells HGC-27 were labeled using Cell Plasma Membrane Staining Kit with DiO (Green Fluorescence) (Beyotime, China), and incubated at 37 • C in the dark for 20 min. Then the labeled cells were washed with PBS, and digested with trypsin. Subsequently, the cells were resuspended in fresh RPMI-1640 medium (10% FBS, 1% PS) to prepare a cell suspension of 3 × 10 6 cells/mL liquid. Larva zebrafish at 2 days post fertilization (dpf) were randomly selected for microinjection. The labeled HGC-27 cells were injected into the yolk sac of the zebrafish. After 24 h, the larva zebrafish were placed under a stereofluorescence microscope (Axio Zoom.V16, Carl Zeiss) to observe and image their tumors (63×) at 1 day post-injection (dpi).
The selected tumor-bearing zebrafish were randomly grouped and placed in a 12-well plate with 10 embryos per well. CKI 150, 50, 25, 0 µg/mL (model group) and cisplatin 100 ng/mL (positive control group) were given, respectively. Zebrafish embryo culture water was used for drug preparation and dilution. The treated zebrafish were incubated in a 31 • C incubator, and the drug was changed every 24 h. With 4 dpi as the end point of the experiment, the tumor-bearing zebrafish were photographed again with a stereoscopic fluorescence microscope. The tumor fluorescence area and fluorescence intensity in the zebrafish yolk sac area were quantified by Image pro-Plus 6.0 software. Using the tumor fluorescence area and integral optical density (IOD) at 1 dpi as a starting point, the tumor growth rate was calculated.

Cell Viability Assay
The viability of HGC-27 and BGC-823 cells was detected by Cell Counting Kit-8 (CCK-8, Dojindo, Japan) assay. Cells were blown into single cell suspension (0.6 × 10 4 cells/mL) after routine digestion and seeded in 96-wells plates (100 µL/well), and cultured for 24 h, routinely. Then the cells were cultured with a drug-containing medium for 24, 48, and 72 h, respectively. After the drug treatment, the CCK-8 solution was added into 96-well plates (10 µL/well) and incubated at 37 • C for 4 h. The optical density (OD) was detected at 450 nm by using a microplate reader (Molecular Devices, United States).

Proteomics Analysis
The BGC-823 cells after CKI intervention were used for Tandem Mass Tag (TMT) labeling proteomics and Proteome Discoverer 1.4 software was performed for identification and quantitation analysis. A detailed description of the method was shown in the Supplementary Method 1.

Western Blot Assay
Cells were collected in RIPA lysis buffer and centrifuged at 12000 rpm and 4 • C for 20 min. The supernatants were preserved and used for western blot assay. Total protein concentration was gauged by BCA Protein Assay Kit (Solarbio, China). 10 µg of total protein was mixed with 5× sample buffer, boiled at 99 • C for 5 min and loaded onto 10% SDS-PAGE gels. Then the protein bands were transferred onto NC membranes and blocked with 5% non-fat milk for 2 h at room temperature. The NC membranes with proteins were incubated with diluted primary antibodies (Proteintech, China) at 4 • C overnight, including anti-AKR1B1 (1:1000), anti-PTGER3 (1:1500), anti-MMP2 (1:1000) and anti-GAPDH (1:2000) antibodies. Then, membranes were incubated with relative sources of secondary antibodies (1:5000) at room temperature for 1.5 h. At last, the specific protein bands were recognized with immobilon western chemiluminescent HRP substrate (MilliporeSigma, United States). Image J software was used for image analysis and the signals of specific proteins were normalized to GAPDH or β-Tubulin.

Statistical Analysis
Data were presented as mean ± SD and statistical analysis was performed with the two-tailed unpaired Student's t-test using GraphPad Prism 8.0 software. In all statistical analyses, statistical significance was pinpointed by a single asterisk ( * : P < 0.05), two asterisks ( * * : P < 0.01).

Identification of Major Compounds in Compound Kushen Injection by UPLC-MS
In this study, 10 marker ingredients of CKI were identified by UPLC-MS (Figure 2 and Table 1). Simultaneously, we also supplemented the chemical ingredients by a literature research (Ma et al., 2013;Wang et al., 2015). Finally, a total of 16 active ingredients in CKI were selected for the next in-depth study, which included 9α-hydroxymatrine, adenine, baptifoline, isomatrine, lamprolobine, piscidic acid and 10 compounds detected by UPLC-QE-Orbitrap-MS chromatography, and the three-dimensional structures of these active ingredients were derived from PubChem database (Kim et al., 2016).

Compound Kushen Injection-Predicted Target Network
The active compound-predicted target network ( Figure 3A) consists of 301 nodes (16 compound points and 285 gene points) that constitute 636 active compound-predicted target linkages. Details of active compound-predicted targets can be viewed in Supplementary Table 6.

Screening of Differential MicroRNA in Gastric Carcinoma
In this study, GSE23739 and TCGA were adopted to analyze the DEMs in GC. | log 2 FC| ≥ 1, P value < 0.05 and adjust P value < 0.05 were considered statistically significant for the DEMs. For GSE23739, a total of 13 up-regulated gene and 15 down-regulated genes were found. Furthermore, there are 107 up-regulated genes and 56 down-regulated genes in the TCGA analysis (Figures 3B,C). Overlapping DEMs (hsa-miR-20a, hsa-miR-30a, hsa-miR-21, hsa-miR-145) between the GSE23739 and TCGA analysis were retained for further study.

Construction and Screening of Weighted Gene Co-expression Network Analysis Key Modules
After normalization, no outlier samples were eliminated in the present study. To build a scale-free network, the power of β = 6 (scale free R 2 = 0.85) was chosen as the softthresholding parameter (Figures 3D,F) Figure 3E, the blue, turquoise and brown modules were highly correlated with clinical traits and were identified as key modules. Figure 3G indicated the topological overlap measurement (TOM) heat map of adjacency or topological overlap. TOM plot was made up of randomly selected 400 genes. Each row and column represented a module and the genes of the module. The TOM of co-expressed RNA in the key modules was high, and the internal RNA correlation was also stronger. The network building the key modules was filtered with a weight Cutoff = 0.1 between genes. The blue module consists of 632 genes and 74485 gene linkages. The turquoise module consists of 1239 genes and 126102 gene linkages and the brown module consists of 232 genes and 11316 gene linkages (Supplementary Tables 3-5).

Construction of Competing Endogenous RNA Network of Compound Kushen Injection Intervention in Gastric Carcinoma
The intersection of the WGCNA key module network and the hub DEMs prediction target constitutes the ceRNA Network of the lncRNA-miRNA-mRNA Axis in GC. Moreover, the ceRNA network and CKI-predicted target network were merged by Cytoscape to obtain the prediction of the ceRNA network of CKI intervention in GC ( Figure 4A). As shown in Figure 4B, the prediction of the ceRNA network of CKI intervention in GC involved 73 nodes and 203 linkages between genes. The overlapping targets (AKR1B1, TLR4, ESR1, PRKCQ, PIK3CD, CTSK, MMP2, ADRB2, PDE1C, ITGB3, PDE10A, PTGFR, AR and PTGER3) were considered as the key genes for the CKI treatment of GC.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis
A total of 14 putative targets were uploaded to the STRING database to identify the functional partnerships and interactions between them. The key genes and their interacting proteins from the PPI network were showed in Supplementary Figure 1.
To further interpret the function of the key gene, KEGG and GO annotations were performed in the R software. A total of 127 GO entries were identified, including 93 biological process (BP), 25 molecular function (MF), and 9 cellular component (CC) (FDR < 0.01 and P < 0.01). The top ten GO terms were tissue homeostasis, anatomical structure homeostasis, toll-like receptor signaling pathway, bone resorption, intracellular receptor signaling pathway, multicellular organismal homeostasis, integrin complex, a protein complex involved in cell adhesion, tissue remodeling, response to ketone (Figures 4C-E). The KEGG results demonstrated that 37 entries satisfy FDR < 0.05 and P < 0.05 ( Figure 4F). These targets were significantly enriched in many pathways related to cancer and signaling pathways, such as the PI3K/AKT signaling pathway. In addition, the Toll-like receptor signaling pathway related to immunity was also significantly enriched ( Figure 4G). We also conducted a modular analysis of the lncRNA-miRNA-mRNA Axis intervened by CKI in Cytoscape by Mcode. A total of key modules were analyzed, including CTSK, ITGB3, PTGER3, hsa-miR-20a-5p and hsa-miR-30a-5p five targets ( Figure 4H). To gain insights into the pharmacological mechanisms of CKI on GC, we performed KEGG analysis for two key miRNAs. The results illustrated that the validated targets of hsa-miR-20a-5p and hsa-miR-30a-5p were associated with signaling pathways closely related to cancer and development, such as Pathways in cancer, Hippo signaling pathway and p53 signaling pathway ( Figure 4I).

The Mechanism of Compound Kushen Injection Anti-Gastric Carcinoma in Cellular Proteome
To further explore the potential mechanism underlying CKI in GC treatment, TMT labeling proteomics was utilized to detect proteome differences between CKI-treated and CKI-untreated GC cells. Quality control, identification and quantification data of proteomics were shown in Supplementary Figure 2. The number of up-regulated DEPs was 490, and down-regulated DEPs was 304 (Figures 5A,B). 210 targets were shared by the proteomic results with those predicted by main active compounds and 561 linkages were attributed to main active compounds and common targets (Supplementary Table 7), which supported the credibility of our network pharmacology results.
Subcellular organelles are micro-organs with specific morphology and function, such as mitochondria, endoplasmic reticulum, etc. They are important sites for proteins to realize various functions. Different subcellular organelles often perform different cellular functions, the analysis of subcellular localization of proteins is helpful to explore the functions of proteins in cells. We found that 416 DEPs were located in the nucleus, 205 DEPs were located in the cytoplasm, 124 DEPs are located in plasma membrane and 119 DEPs were located outside the cells (Figure 5C).
Interactions between proteins and proteins or other ligands with low molecular weight are often based on domains, so domain prediction is of great significance for studying the key functional regions of proteins and their potential bioactivity. Interproscan, a domain prediction software, was used to predict the domain of differentially expressed proteins. Figure 5E shows the first 20 protein domains. And Figure 5F exhibited the domain enrichment characteristics of DEPs, which determine the significance level of protein enrichment.
Proteins, that play a biological regulatory role can interact with other proteins to achieve their function. Systematic analysis of the interaction relationships of a large number of proteins in biological systems is important to understand how proteins work in biological systems, and the response mechanisms of biological signals and metabolism of energy substances in special physiological states, for example, the diseases, and the functional links between proteins. DEPs were brought into the STRING database for PPI annotation, and the PPI network was constructed by Cytoscape software. This network consists of 709 protein nodes and 3125 connections ( Figure 5D).
To comprehensively understand the functions, localization and biological pathways involved in DEPs, GO was used for the biological annotation of proteins. Blast2GO is a highthroughput functional annotation and data mining software. And its features are the combination of various annotation strategies and tools controlling type and intensity of annotation, and the numerous graphical features such as the interactive GOgraph visualization for gene-set function profiling or descriptive charts (Götz et al., 2008). Blast2Go software was utilized to annotate the GO function of DEPs, and counted the number of DEPs at the GO secondary function annotation level ( Figure 6A). Fisher's exact test was applied to analyze the GO functional enrichment of DEPs. The overall functional enrichment characteristics of DEPs were revealed by evaluating the significant level of protein enrichment of a GO functional item. As flashed in Figure 6B, the enrichment analysis revealed that the top five enrichment entries of DEPs were associated with multiple biological processes (BP, including mitotic karyokinesis, regulation of mitotic karyokinesis, sister chromatid separation, mitotic sister chromatid segregation, respiratory electron transport chain), cell compositions (CC, including extracellular matrix containing collagen, mitochondrial respiratory chain complex I, NADH dehydrogenase complex, respiratory chain complex I, oxidoreductase complex) and molecular functions [MF, including structural components of the extracellular matrix, ubiquitin-like protein ligase binding, NADH dehydrogenase activity, NADH dehydrogenase (panquinone) activity, NADH dehydrogenase (quinone) activity]. Figure 6C exhibited the top 25 GO entries (P value), including 18 BP entries, 6 CC entries and 1 MF entry.
Kyoto Encyclopedia of Genes and Genomes pathway annotation and enrichment analysis were utilized to clarify the protein-related pathways of DEPs. KEGG Automatic Annotation Server (KAAS) is an online annotation tool provided by KEGG database itself, which is convenient and precise (Moriya et al., 2007). KEGG pathway annotation revealed that DEPs were enriched in cancer-related pathways, cell growth and death-related pathways, such as PI3K-AKT signaling pathway, MAPK signaling pathway, etc. (Figure 6D). KEGG pathway enrichment analysis exhibited that the top 10 pathways (P < 0.05) were the retrograde endogenous cannabinin signaling pathway, complement and coagulation cascade pathway, non-alcoholic fatty liver disease pathway, oxidative phosphorylation pathway, ABC transporter pathway, glycosphingolipid biosynthesis-ganglion series pathway, ECMreceptor interaction pathway, fever pathway, aldosterone synthesis and secretion pathway, and glycosaminoglycan degradation pathway ( Figure 6E).

Meta-Analysis of Key Genes
We selected microarray datasets of gastric cancer tissues from the GEO database for meta-analysis to demonstrate the differential expression of key genes in gastric cancer tissues. A total of 8 microarrays from the GEO database met the entry criteria. The features of the included GEO datasets are depicted in Table 2. The expression data from the tumor and control groups were collected based on the GEO database. Meta-analysis was conducted based on the expression data of the 8 included microarrays. The results ( Figure 7A) demonstrated that 7 of the key genes exhibited remarkable abnormal regulation in the gastric cancer groups. Given the apparent heterogeneity, a random effects model was applied. On the one hand, ADRB2 (SMD = −1.46; 95% CI −2.02, − 0.91; P < 0.01), PDE1C (SMD = −0.75; 95% CI -1.11, −0.39; P < 0.01) and PTGER3 (SMD = −0.58; 95% CI −1.08; −0.07; P < 0.01) were downregulate in the cancer groups, which might be a tumor suppressor gene. On the other hand, AKR1B1 (SMD = 0.3; 95% CI 0.01; 0.60; P < 0.01), CTSK (SMD = 1.52; 95% CI 0.98; 2.06; P < 0.01), MMP2 (SMD = 1.02; 95% CI 0.51; 1.53; P < 0.01), TLR4 (SMD = 0.85; 95% CI 0.34; 1.37; P < 0.01) were up-regulate in cancer groups, which might be the tumor proto-oncogene. Later, a sensitivity analysis was performed to explore whether a particular microarray played a vital role in the significant heterogeneity. It was found that none of the included studies had a decisive role. A funnel plot was generated to estimate publication bias ( Figure 7B). The points in the funnel were asymmetrically distributed on both sides of the midline, indicating that the bias was mainly related to publication bias, but there could be other reasons such as the lack of included literature ( Figure 8A).

Survival Analysis of Key Genes
A Kaplan-Meier curve was later used to identify the effects of hub genes expression on survival time. As shown in Figure 8B, AKR1B1 (p = 0.000988), AR (p = 0.0102), ITGB3 (p = 0.0389), MMP2 (p = 0.0465), PTGER3 (p = 0.0449), and PTGFR (p = 0.0439) with the p values were all less than 0.05, suggesting that these genes may be the key targets affecting the survival of GC patients.

Immunoinfiltration Analysis
Analysis using TIMER showed that hub genes were negatively associated with purity, and ADRB2 (cor = −0.275) was most negatively correlated with tumor purity. In addition, key genes were strongly correlated with macrophages and dendritic cells. Where PTGER3 correlated most strongly with macrophages (cor = 0.637) and CTSK (cor = 0.624) for dendritic cell (Table 3 and Figure 8C). The Univariate Cox survival analysis showed that among the six types of immune cells, only macrophages were associated with survival of GC patients, which was an indicator of survival of GC patients.

Molecular Docking
Docking studies were carried out between CKI and hub genes. The 3D protein structures of PTGFR and PDE1C were not found in the PDB database. The molecular docking results are shown in Table 4. AR, ITGB3, AKR1B1, ADRB2, and PTGER3 were five genes with the highest predicted for interaction between each of the five protein targets and corresponding CKI components.

Anti-Gastric Carcinoma Effect of Compound Kushen Injection
Zebrafish has great advantages in real-time monitoring of tumor cell growth in vivo and is widely used in human cancer research through the tumor xenograft technique. In this study, different concentrations of CKI and the positive drug cisplatin were used to treat tumor-bearing zebrafish. After 72 h of continuous exposure to the drug, the fluorescence area and intensity of the tumor in the yolk sac area of the low-, medium-, and high-dose groups and the positive control group were decreased, compared with the model group, indicating that 25-150 µg/mL CKI had an inhibitory effect on the growth of the transplanted tumor (Figure 10). And the inhibitory effect was concentration-dependent. Figure 11A, showed that CKI treatment had no significant inhibitory effect for 24 h, except for the 2.0 mg/mL group. At 48 h, 1.0 and 2.0 mg/mL had a significant inhibitory effects on BGC-823 cells. After 72 h, concentrations above 0.5 mg/mL showed significant inhibitory effects. As shown in Figure 11D, the concentration of CKI reached above 0.25 mg/mL, it began to produce apparent inhibitory effect on cell proliferation of HGC-27 cells. The inhibitory effect on cell proliferation gradually increased with the increase of drug concentration.

Effect of Compound Kushen Injection on the Expression of AKR1B1, MMP2, PTGER3 in Gastric Carcinoma Cells
BGC-823 and HGC-27 cells were used to analyze the regulatory effects of CKI on key genes (AKR1B1, MMP2, and PTGER3) in GC cells to evaluate the mechanism of CKI in GC treatment. After 48 h of CKI treatment, compared with the control group in BGC-823 cells (Figure 11B), mRNA expression levels of AKR1B1 decreased significantly in CKI 2.0 and 4.0 mg/mL groups (P < 0.05), MMP2 declined significantly in CKI 1.0 and 2.0 mg/mL groups (P < 0.05), and PTGER3 increased significantly in all CKI groups (P < 0.05). Similarly, we obtained the same validation in HGC-27 cells (Figure 11E). To further confirm the regulatory effects of CKI on key genes (AKR1B1, MMP2, and PTGER3), a western blot assay was used to detect the extent of protein variation intervened by CKI in BGC-823 and HGC-27 cells. After 48 h of CKI treatment, protein expression levels of AKR1B1 were significantly reduced (P < 0.05) in  BGC-823 cells compared with the control group ( Figure 11C) in the 2.0 and 4.0 mg/mL CKI groups, MMP2 was significantly decreased (P < 0.05) in all CKI groups, and PTGER3 was significantly increased (P < 0.05) in the 1.0 and 2.0 mg/mL CKI groups. Approximately, we obtained the comparable validation in HGC-27 cells ( Figure 11F).

DISCUSSION
Gastric cancer is one of the most common cancers and its mortality rate remains high (Waldum et al., 2017). Since GC is difficult to detect in the early stages, there is always a delay in diagnosis in patients with gastric cancer . Therefore, due to the crisis situation of GC with low chance of cure and poor prognosis, there is an urgent need to develop new therapies. CKI is a prescribed drug approved by the Chinese Medicine Administration of China (NMPA). It has also passed the standardized Good Manufacturing Process (GMP) certification and is widely used clinically to treat gastric cancer (Zhao et al., 2014). According to several systematic reviews and meta-analysis studies, it has been found that CKI can not only improve the clinical efficacy in gastric cancer patients but also mitigate the adverse effects of radiotherapy and chemotherapy (Zhang et al., , 2018Huang and Wei, 2019).
Cancer is a complex disease that results from changes in multiple biological networks (Pache et al., 2012). Therefore, in this study, we established an advanced systems pharmacology strategy that was used to reveal the mechanism underlying the effects of CKI on GC. In detail, this initiative strategy was integrated WGCNA, molecular target prediction technology, microarray analysis, meta-analysis, molecular docking technology, proteomics and experimental verification in vivo and in vitro. The high-throughput data analysis method  was used to find miRNAs closely associated with gastric cancer for target prediction and overlapped with the key modules of WGCNA analysis in TCGA to identify ceRNA networks closely related to GC. Finally, considering network pharmacology as a core concept, the key target of CKI in gastric cancer was identified, and 14 intersection genes were identified as hub genes. To further explore the impact of CKI on gastric cancer, we conducted a meta-analysis of key targets to compare the differential expression of key genes in gastric cancer tissues and normal tissues. Second, we performed functional analysis to understand the biological regulatory pathways involved in key genes. In addition, survival analysis and immune infiltration analysis were used to analyze the relationship between key genes and the prognosis of gastric cancer. Finally, molecular docking simulation was used to verify the binding of CKI components to key targets.
Network pharmacology can explain the effects of drugs on the disruptions of biological networks from the perspective of macro or overall regulation, and explain the treatment of diseases from the perspective of multi-component-multi-target-multi-pathway (Jing et al., 2019). We discovered through network pharmacology that 14 intersection genes could be the key targets for CKI treatment of GC. After meta-analysis of GC gene expression profile chip, it was found that 7 of these genes, including AKR1B1, CTSK, MMP2, TLR4, ADRB2, PDE1C, and PTGER3, had significant differences in gastric cancer tissues. Besides, AKR1B1, MMP2 and PTGER3 were found to be important in the analysis of gastric cancer survival, so the above three genes are considered to be the most significant hub genes for CKI to treat gastric cancer and improve the prognosis of GC. AKR1B1 as a common highexpressed gene in cancer, including gastric cancer, may lead to increased proliferation, metastasis and invasion of tumor cells by driving the epithelial-tomesenchymal transition (EMT) Schwab et al., 2018). Studies have shown that ADRB2 can directly interact with and upregulate AKR1B1 in pancreatic cancer cells, promoting cell proliferation and inhibiting apoptosis through the ERK1/2 pathway (Xiao et al., 2018). In addition, the expression of MMP2 in AKR1B1 knockdown cancer cells also decreased significantly compared with the control group (Schwab et al., 2018). MMP2 has been implicated in tumor FIGURE 11 | The expressions of AKR1B1, MMP2, and PTGER3 were detected in GC cells. In BGC-823 cells, (A) Dose-inhibition curves of CKI. For 24 h, IC 50 = 1.72 ± 0.22 mg/mL; for 48 h, IC 50 = 1.20 ± 0.11 mg/mL; for 72 h, IC 50 = 0.79 ± 0.05 mg/mL. (B) mRNA level of key genes was measured by RT-qPCR after CKI intervention. (C) Protein level of key genes was determined by western blot after CKI intervention. In HGC-27 cells, (D) Dose-inhibition curves of CKI. For 24 h, IC 50 = 0.98 ± 0.17 mg/mL; for 48 h, IC 50 = 0.64 ± 0.13 mg/mL; for 72 h, IC 50 = 0.62 ± 0.22 mg/mL. (E) mRNA level of key genes was measured by RT-qPCR after CKI intervention. (F) Protein level of key genes was determined by western blot after CKI intervention. Data were presented as mean ± SD. n = 3. * P < 0.05; * * P < 0.01. development and morphogenesis (Zhao et al., 2019). It has been demonstrated that MMP-2 can regulate the extracellular matrix (ECM) degradation, which plays an important role in cancer development (Zhang et al., 2012). Previous studies have confirmed that matrine, an important component of CKI, can downregulate the abnormal expression of MMP2 and thus inhibit tumor cell invasion and metastasis (Qian et al., 2015;Huang et al., 2017;Gao et al., 2018). Whole-genome analysis showed that PTGER3 is abnormally low in gastric cancer. PTGER3 can inhibit the secretion of gastric parietal cells and gastric acid (Kim et al., 2018). The lack of PTGER3 leads to abnormal secretion of gastrin and gastric acid and accelerates the occurrence of gastric cancer (Nishio et al., 2007). In addition, PTGER3 can also upregulate the expression of related MMP2 and AR to promote the proliferation of cancer cells and the deterioration of gastric cancer (Robertson et al., 2010;Kashiwagi et al., 2013). Indeed, we also found that CKI could induce the changes of intracellular AKR1B1, MMP2 and PTGER3 at mRNA and protein levels in GC cell lines BGC-823 and HGC-27, and the change trends were consistent in these two levels. Therefore, we consider that these may be important approaches of CKI treatment of GC.
In this study, we performed GO enrichment and KEGG pathway analysis to elucidate the multiple mechanisms of CKI against GC from a systematic level. The key genes of CKI on GC were enriched in the PI3K/AKT signaling pathway, Tolllike receptor signaling pathway and other pathways, as indicated by the functional enrichment analysis. With frequent alterations identified in GC, the PI3K/AKT pathway is significantly involved in gastric carcinogenesis and progression (Matsuoka and Yashiro, 2014). The PI3K/AKT pathway can be activated by various factors, including hormone and ECM signaling pathways, thereby regulating several basic cellular activities such as cell proliferation, apoptosis and metastasis (Fresno et al., 2004). As the most common dysregulatory pathway in cancer, the PI3K/AKT pathway has received increasing attention due to its potential for targeted therapy in many malignancies (Hu et al., 2019). We proposed that CKI plays a therapeutic role in the treatment of GC mediated by the PI3K/AKT pathway, and this has been previously confirmed experimentally. Previous studies have confirmed that the active ingredients of matrine, oxymatrine and sophoridine in CKI can treat various tumors by inhibiting the PI3K/AKT signaling pathway Dai et al., 2018;. Peng et al. (2016) found that matrine can inhibit the proliferation and metastasis of gastric cancer cell SGC7901 via PI3K/Akt pathway. Accordingly, by comparing the proteomics results between the CKI group and the control group, we found that DEPs were significantly enriched in the PI3K/AKT signaling pathway, MAPK signaling pathway, ABC transporter pathway, ECM receptor interaction pathway and other signaling pathways. Excessive activation of the PI3K/AKT signaling pathway leads to over-expression of P-gp (ABC carrier), which affects the normal function of the ABC transporter signaling pathway (Sui et al., 2014). Abnormal ABC transporter pathway is common in multidrug resistance. MAPKs transduce signals through tertiary kinase cascades, and activated MAPKs maintain normal cellular functioning in the nucleus (Meister et al., 2013). MAPK subgroups, including ERK and JNK, and their signaling pathways are intimately associated with the regulation of ABC transporters, affecting multidrug resistance Ji et al., 2015). Accordingly, we speculate that PI3K/AKT, MAPK, ABC transporters, ECMreceptor interaction pathways and the crosstalk between them are the potential mechanisms of CKI that regulate the basic activities (proliferation, migration, invasion) of GC cells and reverse the drug resistance of GC cells.
From a comprehensive analysis of the above results, we also found that immunization may be an important potential influencing factor in CKI treatment of GC, therefore we performed an analysis of the immune invasion of key genes for gastric cancer. The analysis revealed that the degree of macrophage infiltration affects the prognosis of GC patients, and the expression of key genes is positively correlated with the degree of macrophage infiltration. Consistent with this, sophoridine has been shown to polarize tumor-associated macrophages (TAMs) into M1-TAMs and suppress M2-TAMs polarization through the TLR4/IRF3 axis. In addition, it can inhibit the migration ability of macrophages and reshape the immunological microenvironment of gastric cancer (Zhuang et al., 2020). Taken together, we propose that CKI can not only directly inhibit the proliferation and metastasis of gastric cancer tumor cells, but also improve the prognosis of cancer patients through immunotherapy.
In addition to the key mRNAs, we found two significant miRNAs possibly involved in CKI regulation of gastric cancer, namely hsa-miR-20a-5p and hsa-miR-30a-5p, by module analysis. KEGG pathway analysis of miRNA's demonstrated that most of the target genes were enriched in the Hippo signaling pathway and p53 signaling pathway, which are closely associated with cancer cell proliferation and metastasis. Matrine plays an important role in regulating the p53 signaling pathway to inhibit cell proliferation in liver cancer, lung cancer and esophageal cancer Xie et al., 2015;Lu et al., 2017). Besides, matrine can also promote apoptosis of colorectal cancer cells via the Hippo signaling pathway . It is worth mentioning that the study found that sophoridine can significantly activate the Hippo and p53 signaling pathways and inhibit the progression of lung cancer and enhance the effect of the anticancer drug cisplatin against lung cancer cells (Zhu et al., 2020). However, there are few studies on the direct effect of CKI and its active ingredients on miRNA in the treatment of gastric cancer, which needs to be confirmed by experiments in vivo and in vitro.

CONCLUSION
This study explored the potential molecular mechanism of CKI in the treatment of GC and established a new advanced system pharmacology strategy. Based on the traditional network pharmacology, the potential molecular mechanism of CKI in the treatment of GC was initially determined by the analysis of high-throughput chip data and WGCNA. Moreover, chip meta-analysis methods and survival analysis were used to verify the expression and prognosis of key genes in GC. Functional enrichment analysis and immune infiltration analysis focus on the functional impact of key genes. Finally, molecular docking was employed to verify the strong binding between the target and the components. By using network pharmacology combined with multiple integrated bioinformatics methods, we systematically revealed that CKI might be involved in regulating the ceRNA network for the treatment of GC. Among them, mRNA (including AKR1B1, MMP2 and PTGER3) and miRNA (including hsa-miR-20a-5p and hsa-miR-30a-5p) may play an essential role in the therapy. Our preliminary conclusion is that CKI could be used for GC therapy by activating signaling pathways such as PI3K/AKT, MAPK, ABC transporter and ECMreceptor interaction pathways to inhibit cancer cell proliferation and regulate immunity. Based on a multidisciplinary approach, this study might provide a new perspective for the profound exploration and provide a reference for multicomponentmultitarget-multipathway clinical research.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The experiment was in accordance with the Animal Management Rules of the Ministry of Science and Technology of the People's Republic of China for experimental care and use of animals and approved by the Animal Ethics Committee of Beijing University of Traditional Chinese Medicine.

AUTHOR CONTRIBUTIONS
WZ, JW, and RY designed the experiments and wrote the manuscript. WZ and XL collected and analyzed the bioinformatics data. CW, CZ, ZH, and XF performed most of the experiments. SL and ZW interpreted the experimental data and visualized the results. YT, AS, and JZ substantively revised the manuscript. All authors gave the final approval of the version to be published.

FUNDING
The design of the study and the collection, analysis, and interpretation of data were supported by the National Nature Science Foundation of China [Grant No. 82074284], Cooperation project between Beijing University of Chinese medicine and Shanxi Zhendong Pharmaceutical Co., Ltd. [Grant No. BUCM-ZD 2020001], and the Young Scientists Training Program of Beijing University of Chinese Medicine [Grant No. BUCM-QNLJ 2019001].