Computational Prediction of Antiangiogenesis Synergistic Mechanisms of Total Saponins of Panax japonicus Against Rheumatoid Arthritis

Objective: To investigate the anti-angiogenesis mechanisms and key targets of total saponins of Panax japonicus (TSPJ) in the treatment of rheumatoid arthritis (RA). Methods: RStudio3.6.1 software was used to obtain differentially expressed genes (DEGs) by analyzing the differences in gene expression in the synovial tissue of RA and to predict the potential targets of active compounds from TSPJ by the PharmMapper and SwissTargetPrediction databases. We evaluated the overlapping genes by intersectional analysis of DEGs and drug targets. Based on the overlapping genes, we used Cytoscape 3.7.2 software to construct a protein–protein interactions (PPI) network and applied Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to determine the mechanisms of the treatment. Finally, the correlations with angiogenesis-related genes were explored. Collagen-induced arthritis (CIA) model was established and treated with different doses of TSPJ. The manifestations of CIA were determined by evaluation of arthritis index and histology score. Serum levels of vascular endothelial growth factor (VEGF) and the hypoxia-inducible factor 1 (HIF-1) were tested by ELISA. The mRNA levels of IL-1β and IL-17A were detected by real time-quantitative PCR. Results: Altogether, 2670 DEGs were obtained by differential analysis, and 371 drug targets were predicted for four active components (Araloside A, Chikusetsusaponin IVa, Ginsenoside Rg2, and Ginsenoside Ro). A total of 52 overlapping genes were included in the PPI network and the KEGG analysis. However, only 41 genes in the PPI network had protein interactions. The results of the KEGG enrichment analysis were all related to angiogenesis, including VEGF and HIF-1 signaling pathways. Seven genes with negative correlations and 16 genes with positive correlations were obtained by correlational analysis of DEGs in the VEGF and HIF-1 signaling pathways. SRC proto-oncogene, nonreceptor tyrosine kinase (SRC), and the signal transducer and activator of transcription 3 (STAT 3) had a higher value of degree and showed a significant correlation in the pathways; they were regarded as key targets. Compared with the model group, TSPJ significantly relieved the symptoms and decreased the expression of VEGFA, HIF-1α, IL-1β, and IL-17A in serum or spleens of CIA mice. Conclusion: In the current study, we found that antiangiogenesis is one of the effective strategies of TSPJ against RA; SRC and STAT 3 may be the key targets of TSPJ acting on the VEGF and HIF-1 signaling pathways, which will provide new insight into the treatment of RA by inhibiting inflammation and angiogenesis.


INTRODUCTION
Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease affected by both genetic and environmental factors (Croia et al., 2019). Approximately 0.5-1% of the world's population suffers from RA, and many patients have symptoms of joint damage or even disability (Smolen et al., 2016;Aletaha and Smolen, 2018). When synovitis occurs, macrophages and fibroblasts in synovial tissue are induced to secrete angiogenic factors, leading to an abnormal increase of synovial vessels. A large number of inflammatory cells enter the joint through newly formed blood vessels, causing synovitis and synovial hyperplasia. Increased angiogenesis leads to the permanence of RA (Elshabrawy et al., 2015;Leblond et al., 2017;Lu et al., 2017;Chen et al., 2018). Eventually, these mechanisms destroy the joint and lead to deformity. It is worth noting that the treatment of RA with vascular targeted therapy has gradually entered the clinic (Leblond et al., 2017).
Recently, many people have become interested in research into traditional Chinese medicine (TCM) because it plays an important role in the treatment of difficult and complicated diseases (Yuan et al., 2017). The study of Chinese herbal medicine, with the identification of the main components of herbals and the molecular mechanisms of their curative effects, has become a hot research topic. Currently, TCM is increasingly being used in clinical treatment or for developing new treatment methods . Panax japonicus (T.Nees) C.A.Mey is a kind of Chinese herbal medicine that belongs to the Panax genus, which is widely grown in the southwest of China. Studies have shown that Panax japonicus (T.Nees) C.A.Mey has analgesic, antiinflammatory, antioxidant, and joint swelling-relieving properties (Yang et al., 2014;Yang et al., 2018). Saponins are the main components of Panax japonicus (T.Nees) C.A.Mey. Total saponins of Panax japonicus (TSPJ) are extracted from Panax japonicus (T.Nees) C.A.Mey and are widely used to treat RA. Its main ingredients include Araloside A, Chikusetsusaponin IVa, Ginsenoside Rg2, and Ginsenoside Ro (Yu-min et al., 2019) (Supplementary Figure  S1 and Table S1). Ding et al. have confirmed that Araloside A has an anti-inflammatory effect in the treatment of RA . Chikusetsusaponin IVa and Ginsenoside Ro exhibit antiangiogenic effects (Zheng et al., 2019). In the treatment of RA, TSPJ can inhibit synovial hyperplasia and reduce capillary permeability (Guo et al., 2019). Yuan et al. found that TSPJ could reduce the degree of foot swelling (Ding et al., 2008) and downregulate the serum levels of TNF-α and IL-1β in CIA rats (Ding et al., 2009). Related studies have shown that TSPJ can inhibit angiogenesis during tumor therapy (Yougui, 2011). However, no previous studies have shown whether TSPJ can be used as a targeted drug against angiogenesis in RA.
With the development of biomedical big data research, the network pharmacology approach and data mining in bioinformatics is being applied to help study the potential mechanism of different actions of drugs in organisms (Boezio et al., 2017) and has been proven to effectively reveal the regulation principles of small molecules of TCM (Zhang R.-z. et al., 2017). There is a large quantity of biological data in the public database, among which there is much valuable information that has not yet been highlighted. Bioinformatic methods are widely used to extract productive data because of their advantages of being quick and inexpensive (Ramharack and Soliman, 2018). By constructing an interaction network of disease-related proteins of TCM components, network pharmacology can identify the key targets in the network, and study the interactions among the drugs, targets, and pathways (Yuan et al., 2017;Sidders et al., 2018). In this study, we first identified differentially expressed genes (DEGs) of RA related to angiogenesis by bioinformatics analysis. After predicting potential targets for the active compounds of TSPJ, we used a network pharmacology approach to identify the possible antiangiogenic mechanism of the main components of TSPJ in RA therapy. Because CIA model has many similar pathological and immunological characteristics to human RA, many animal experiments have used CIA mice as RA models (Choudhary et al., 2018). We successfully established CIA model and verified partly of our computational prediction results (Supplementary Figures  S2-S5). The findings of this study deepen our understanding of the therapeutic mechanism of TSPJ and provide a new therapeutic idea in the antiangiogenic treatment of RA. The flow diagram of this study is shown in Figure 1.

Data Collection and Preprocessing
The Gene Expression Omnibus (GEO) database (https://www. ncbi.nlm.nih.gov/geo/), a public functional genomics database built and maintained by the National Center for Biotechnology Information (NCBI), collects microarray, next-generation sequencing, and other forms of high-throughput functional genomic data. Persistent synovitis and joint destruction are the main characteristics of RA (Scott et al., 2010). In our study, we mainly studied the difference in synovial tissue expression between RA patients and healthy individuals. Microarray raw data were downloaded from the GEO database by searching with the keywords: "rheumatoid arthritis," "synovial tissue," and "Homo sapiens." All the above data sets were merged according to the same probe name by RStudio3.6.1 software to build the gene expression data matrix. The gene expression data matrix was normalized using the "Limma" R package. In the end, the gene expression matrix of synovial tissue was obtained. Frontiers in Pharmacology | www.frontiersin.org October 2020 | Volume 11 | Article 566129 3 PubChem (https://pubchem.ncbi.nlm.nih.gov/) is an openaccess database that contains chemical structures, molecular formulas and patents, and article identifiers of drugs (Sayers et al., 2020). The chemical structures of the main components of TSPJ were searched in the PubChem database. Based on the chemical structure models, drug targets were predicted via the PharmMapper database (http://www.lilab-ecust.cn/ pharmmapper/) and the SwissTargetPrediction database (http://www.swisstargetprediction.ch/).

Differentially Expressed Genes' Screening
The "Limma" R package is commonly used for variance analysis. First, RStudio3.6.1 software was used to read the gene expression matrix of synovial tissue. All samples were divided into two groups (healthy control and the RA group). Second, the "Limma" package was adopted to identify the DEGs between the healthy control and the RA group. The results also included fold change (FC) and the p values to reflect the degree of gene differences. p values were adjusted for multiple test correction using the Hochberg and Benjamini test. Finally, |log2FC| > 2 and adj.p < 0.05 were used as the cutoff criteria to identify DEGs.

Network Pharmacology
The UniProt database (http://www.uniprot.org/) provides highquality annotation for proteins (The UniProt Consortium, 2017). The UniProt database was utilized for retrieving drug prediction targets information and screening out targets for the "Homo sapiens" species. The overlapping targets between the drug targets and the DEGs were identified, and input into the STRING database (https://string-db.org/) to retrieve the protein-protein interactions (PPI). Cytoscape3.7.2 software was used to import the data to construct and visualize the PPI network. We used the "Network Analyzer" plugin to analyze the network to obtain network-related parameters.

Pathway Enrichment Analysis
The DAVID database (https://david.ncifcrf.gov/) is commonly used for annotation, visualization, and integrated discovery (Dennis et al., 2003). To illustrate the possible mechanism of drug action in the treatment of RA, the overlapping targets were adopted for the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis by the DAVID database with the "Homo sapiens" setting. Finally, the signaling pathway with p < 0. 05 was selected for the next study.

Correlation Analysis
Significant signaling pathways were acquired by KEGG enrichment analysis. To learn more about the synergistic mechanisms of these pathways, we performed correlation analysis and identified negative and positive genes in the related pathways. First, we downloaded all targets in these signaling pathways from the KEGG (https://www.kegg.jp/) database. Second, the overlapping genes between the DEGs and pathways were identified by intersection analysis. Finally, the roles of these overlapping genes in the pathway were determined by correlation analysis.

Gene Expression Matrix
After a systematic review, we found two gene expression data sets that met our requirements. The gene expression data sets GSE48780 and GSE77298 were downloaded from the GEO database (Supplementary Table S2). These data sets are based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). We used the "dplyr" R package to merge the data sets to construct the gene expression matrix of synovial tissue. This matrix contained seven healthy samples and 99 RA samples.

Identification of the Overlapping Targets
A total of four active ingredients were identified from TSPJ, including Araloside A, Chikusetsusaponin IVa, Ginsenoside Rg2, and Ginsenoside Ro. The targets of the four main ingredients ( Figure 3A) of TSPJ were predicted by using the PharmMapper database and SwissTargetPrediction database. We obtained 157 Araloside A targets, 157 Chikusetsusaponin IVa targets, 163 Ginsenoside Rg2 targets, and 127 Ginsenoside Ro targets. Then, 371 drug targets were attained by using the UniProt database to screen and remove repetition. Subsequently, the drug targets and DEGs were set as a background list to get the overlapping targets using Venn diagrams. As shown in the Venn FIGURE 2 | Volcano plot of differentially expressed genes. Note: red represents upregulated genes, blue represents downregulated genes, and gray represents background genes. X-axis represents log fold change, and Y-axis represents log 10-adjusted p-value.

Protein-Protein Interactions Network
According to the related data in the STRING databases, the PPI relationships of the overlapping targets were obtained. Using the Cytoscape3.7.2 software, we constructed the PPI network of the overlapping targets. In the network (Figure 4), nodes represent targets and edges represent interaction relationships existing between two targets. The value of degree indicates the connectivity of the node in the network. The higher the value of the node, the more important the node is in the network. We found a total of 41 nodes and 124 edges in this PPI network after applying the "Network Analyzer" plugin. According to the node degree, some hub targets were identified among the network. The top 10 (degree ≥ 7) targets in the network are SRC protooncogene, nonreceptor tyrosine kinase, epidermal growth factor receptor, signal transducer and activator of transcription 3, vascular endothelial growth factor A, jun proto-oncogene, AP-1 transcription factor subunit, Fos proto-oncogene, AP-1 transcription factor subunit, matrix metallopeptidase 2, protein tyrosine phosphatase receptor type C, annexin A1, protein kinase C beta, matrix metallopeptidase 13, and renin.

Signaling Pathway
A total of 10 significant signaling pathways (p < 0.05) were obtained by KEGG analysis of the 52 drug and DEGs overlapping targets in the DAVID database. The results of pathway enrichment analysis ( Figure 5A) included VEGF, sphingolipid, Ras, Rap1, oxytocin, HIF-1, GnRH, estrogen, ErbB, and calcium signaling pathways. Through the analysis of the pathways, we found that all 10 signaling pathways were related to angiogenesis. Antiangiogenesis may be one of the effective mechanisms of TSPJ in the treatment of RA. Recent studies showed that the VEGF (Elaimy and Mercurio, 2018; Apte et al., 2019) and the HIF-1 signaling pathways (Wang et al., 2016;Serocki et al., 2018) were classic angiogenesis signaling pathways. Next, we downloaded all the genes in the two classic signaling pathways (the VEGF signaling pathway and the HIF-1 signaling pathway) from the KEGG database, and 148 genes were obtained after filtering and deleting the duplicates. Intersection analysis between the 148 genes in the signaling pathways and the DEGs identified 28 genes related to angiogenesis, including six upregulated and 22 downregulated DEGs ( Figure 5B). Five out of the top 10 targets (SRC proto-oncogene, nonreceptor tyrosine kinase, epidermal growth factor receptor, signal transducer and activator of transcription 3, vascular endothelial growth factor A, and protein kinase C beta) in the PPI network were also in the same two pathways, and they may play an important role in the mechanism underlying the antiangiogenesis effect of TSPJ.

Correlation Analysis
To analyze the synergistic mechanisms of the antiangiogenesis effect, we performed a correlation analysis for DEGs in the pathways and calculated the correlation coefficient between the expression of each gene and the others. According to the results of the correlation analysis, the genes with a correlation coefficient greater than 0.5 are defined as positive correlation genes to promote each other, and those with correlation coefficients less than −0.5 are negative correlation genes that inhibit each other. The results of the correlation analysis are shown as a heat Frontiers in Pharmacology | www.frontiersin.org October 2020 | Volume 11 | Article 566129 map ( Figure 6A) by the "Pheatmap" R package. By clustering the results of the correlation analysis, the genes in the red region indicate a positive correlation, the blue region indicates a negative correlation, and the yellow region indicates that the correlation is not significant; seven negative and 16 positive correlation genes were found. We mapped these genes into signaling pathways  Frontiers in Pharmacology | www.frontiersin.org October 2020 | Volume 11 | Article 566129 6 FIGURE 6 | The heat map of correlation analysis (A), VEGF signaling pathway (B), and HIF-1 signaling pathway (C). Note: red or pink represents positive correlation and blue or green represents negative correlation.

DISCUSSION
In this study, 52 overlapping genes were found based on 2670 DEGs and 371 drug targets. In addition, the pathways of the KEGG analysis were all involved in angiogenesis, especially the HIF-1α and VEGF signaling pathways. Two key genes in the PPI network were also in these two pathways. Moreover, we successfully established CIA model and verified partly of our computational prediction results. These results suggested that TSPJ could inhibit angiogenesis by targeting the HIF-1α or the VEGF signaling pathway through multitargets such as SRC and STAT 3, effectively treating RA.
Angiogenesis is one of the main features of RA, and it often occurs in the early stage of inflammation (Lu et al., 2017;. RA synovial macrophages and fibroblasts produce growth factors, cytokines, and chemokines after inflammatory stimulation or hypoxia (Bartok and Firestein, 2010;Konisti et al., 2012). These include the production of HIF-1α after hypoxia and the most famous angiogenic stimulator VEGF (Świdrowska-Jaros and Smolewska, 2018). New angiogenesis aggravates the inflammation and then promotes synovial neovascularization. There is a positive feedback loop between the two (Leblond et al., 2017). Inhibition of angiogenesis has become a new choice for the treatment of RA in recent years (Jia et al., 2018). However, many drugs targeting RA angiogenesis are still under research and development (Alaarg et al., 2017;Feng and Chen, 2018), including Chinese medicine compounds or monomers. For example, sinomenine inhibits angiogenesis through the HIF-1α-VEGF-ANG-1 axis . Tripterygium wilfordii, Hook f. indirectly suppresses VEGF by inhibiting TNF-α (Zhang W. et al., 2017). It has previously been reported that some components of TSPJ have an antiangiogenic effect in tumors; however, the mechanism of the antiangiogenesis effect of TSPJ in the treatment of RA was not clear. Because of the multicomponent combination therapy, it is difficult to identify the specific interactions between the components and a single target. Network pharmacology is consistent with the overall concept and can help to predict the targets and identify the mechanisms of TCM (Lee et al., 2019;Zhang et al., 2019). We use the methods of network pharmacology and bioinformatics analysis to further understand the possible mechanism underlying the antiangiogenic effect of TSPJ.
In this study, four active components (Araloside A, Chikusetsusaponin IVa, Ginsenoside Rg2, and Ginsenoside Ro) of TSPJ and the 371 potential targets of TSPJ were predicted. DEGs in RA synovium were obtained by difference analysis and could be involved in the pathogenesis of RA. Acting on these DEGs may be very helpful for the targeted therapy of RA. The pathways containing these DEGs may also be important targets for the treatment of RA. The results of the intersection analysis showed that 52 of the targets may be related. The protein interactions of the 52 overlapping genes were obtained from the STRING database and network visualized by Cytoscape3.7.2 software. The top 10 genes' (SRC proto-oncogene, nonreceptor tyrosine kinase, epidermal growth factor receptor, signal transducer and activator of transcription 3, vascular endothelial growth factor A, jun proto-oncogene, AP-1 transcription factor subunit, Fos proto-oncogene, AP-1 transcription factor subunit, matrix metallopeptidase 2, protein tyrosine phosphatase receptor type C, annexin A1, protein kinase C beta, matrix metallopeptidase 13, and renin) degree values were not less than seven according to the network visualization analysis. A high degree value indicates that these genes play a key role in the regulatory network of RA proteins in responding to TSPJ therapy. Similarly, the important pathways of TSPJ in the treatment of RA were found by the KEGG analysis of overlapping genes. To our surprise, the 10 signaling pathways were all related to angiogenesis, and the top four genes (SRC proto-oncogene, nonreceptor tyrosine kinase, epidermal growth factor receptor, signal transducer and activator of transcription 3, and vascular endothelial growth factor A) in the PPI network are all involved in the VEGF and HIF-1 signaling pathways. We supposed that TSPJ might treat RA by inhibiting angiogenesis.
The VEGF and HIF-1 signaling pathways are two important pathways of angiogenesis. VEGF-A induces endothelial cell proliferation, migration, and survival by activating VEGFR2 and its downstream signal transduction pathway, thus promoting new angiogenesis (Karaman et al., 2018). Synovial inflammation and hyperplasia consume much oxygen, leaving the tissue in a state of local hypoxia. When hypoxia occurs, activation of the HIF-1 signal pathway leads to an increase in the expression of angiogenesisrelated factors such as VEGF Fallah and Rini, 2019). To further identify the synergistic mechanism of the TSPJ against angiogenesis, we analyzed the correlation of DEGs in the VEGF and HIF-1 signaling pathways. Seven genes showed a negative correlation, and 16 showed a positive correlation. The difference in correlation reflects the competition or mutual promotion among these genes (Xuan et al., 2019). Meanwhile, we found only two genes (SRC and STAT 3) with higher degree values in angiogenesis-related pathways according to the results of Frontiers in Pharmacology | www.frontiersin.org October 2020 | Volume 11 | Article 566129 8 the correlation analysis. STAT 3 is a marker of angiogenesis that interacts with SRC (Sp et al., 2017). Li et al. and L. Claesson-Welsh et al. have shown that the activation of SRC phosphorylation promotes angiogenesis (Claesson-Welsh and Welsh, 2013;Li P. et al., 2018). In addition, the results showed that the degree value of STAT 3 in the PPI network is higher than that of other DEGs in the HIF-1 signaling pathway. Studies have shown that STAT 3 is a potential therapeutic target of RA (Oike et al., 2017). STAT 3 directly regulates VEGF, and blocking STAT 3 can effectively inhibit angiogenesis (Wei et al., 2003;Xu et al., 2005). Therefore, we propose that when SRC and STAT 3 are inhibited by TSPJ, the genes positively related to SRC and STAT 3 are also suppressed, and thus reduce angiogenesis.
To further verify the prediction results above, we successfully established CIA mice model. The data showed that TSPJ could significantly reduce the arthritis index of CIA mice. This indicated that TSPJ can effectively relieve the symptoms of CIA mice. The results of H&E staining and histology score of knee joint sections showed that TSPJ inhibited articular cartilage destruction, synovial hyperplasia, and inflammatory cell infiltration in a dose-dependent manner. In addition, TSPJ could reduce the expression of HIF-1 α and VEGFA in serum. This is consistent with our prediction that TSPJ may act as an antiangiogenic by regulating the VEGF or the HIF-1 signaling pathway. Results also showed that TSPJ could decrease the expression of IL-1β and IL-17A in spleen (Supplementary Figures S2-S5). The above results indicated that TSPJ application might improve joint inflammation, vascular proliferation, and bone erosion.

CONCLUSION
The current study demonstrated that TSPJ may regulate the VEGF or HIF-1 signaling pathways by inhibiting two targets, STAT 3 and SRC, thus inhibiting inflammation and angiogenesis, which will provide a new strategy for the treatment of RA.

DATA AVAILABILITY STATEMENT
All data sets presented in this study are included in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
ZF and ZM conceived and designed the experiments; XG, JJ, XF, and YL analyzed the data; XH provided information support; ZF and XG drafted the manuscript; and ZF, GJ, and ZM amended the final manuscript.

ACKNOWLEDGMENTS
Special thanks to GJ for English language editing.