The target proteins of loganin, ursolic acid, morroniside, paeoniflorin, and albiflorin were screened using three databases. First, targets were screened using the Traditional Chinese Medicine System Pharmacology (TCMSP) database and Analysis Platform (https://www.tcmspw.com/tcmsp.php). Based on the predicted target proteins corresponding to the chemical small molecules identified by TCMSP, the names of the target proteins were matched with the abbreviated gene symbols corresponding to the uniprotID in the uniport database (https://www.uniprot.org). Next, targets were screened using the PharmMapper database (http://www.lilab-ecust.cn/pharmmapper) based on the chemical small molecules identified in the database. Finally, targets were screened using the Swiss database (http://www.swisstargetprediction.ch) based on the chemical small molecules identified in the database. All predicted targets from TCMSP and the top ten targets from PharmMapper and Swiss were selected, and a Venn diagram was generated to identify overlapping targets between databases. In total, 132 unique gene targets of the five compounds were identified between the three databses.

GO and KEGG functional analysis of the previously screened targets was performed using the database for Annotation, Visualization, and Integrated Discovery website (https://david.ncifcrf.gov). p ≤ 0.05 was defined as the critical value of significant enrichment, and results were mapped using R software as bar plots and bubble plots for both GO and KEGG analysis. The “Pathview” package in R software was used to generate a diagram of signaling pathways associated with RA based on the results of KEGG analysis.

The interactions among the target proteins of the selected compounds was mapped by constructing a PPI network to present the direct and indirect regulatory relationship between the targets, using the STRING web site (https://string-db.org). For visual analysis of the target proteins, a visual PPI network diagram was generated using Cytoscape 3.6.1, and the Network Analyzer function in Cytoscape was utilized to evaluate the topology parameters of the network nodes. In addition, compound-target relationships were visualized by constructing a pharmacological network map using Cytoscape.

The potential targets of RA were screened in the following databases: OMIM (https://www.omim.org), GenCLiP3 (http://ci.smu.edu.cn/genclip3/analysis.php), CTD (http://ctdbase.org), and GeneCards (https://www.genecards.org). The targets of RA identified here were illustrated using a Venn diagram, revealing 71 common targets between the four databases (Supplementary Table S1). Then, the 71 targets of RA were compared with those of the five selected active compounds to reveal the common targets.

Based on the active compounds and core targets selected above, the protein structure of the corresponding target was obtained from the PDB database (http://www.rcsb.org/). Molecular docking analysis was performed using the AutoDock Vina software (http://vina.scripps.edu/) (Trott and Olson, 2010).

Commercially available Cornus officinalis (COR) and Paeonia lactiflora (PAE) were purchased from Zhejiang Chinese Medical University Medical Pieces Ltd. (Hangzhou, China), ursolic acid (UA, DX0019) and paeoniflorin (PF, DS0070) were purchased from Desite (Chengdu, China), and dexamethasone (DEX, D1756) was purchased from Sigma-Aldrich (St. Louis, MO). Each compound was dissolved in sterilized water, given that drug administration is performed at 10 ml per kg of animal body weight. All animal experiments were performed at Wuhan Myhalic Biotechnology Co., Ltd. (Wuhan, China). The experimental protocol was approved by the institutional review board of the Model Animal Research Institute at Wuhan Myhalic Biotechnology Co., Ltd. and adhered to the guidelines for animal care and use (approval number: HLK-20190418–01). Male specific-pathogen-free Sprague-Dawley rats weighing 200 ± 20 g, acquired from China Three Gorges University, were housed in a facility with 50–60% relative humidity at 25°C. Before the experiment, the rats adaptively fed for 7 days, where they were allowed free access to food and water. The CIA model was established following previously reported methods (Jia et al., 2014), with three rats in each group. Control rats were not treated in any way. To induce CIA, 10 mg of type II collagen (PAB43878, Bioswamp) was mixed with 0.01 M acetic acid and emulsified with an equal volume of Freund’s complete adjuvant to a final collagen concentration of 2 mg/ml. On day 0, experimental rats were intradermally injected with 0.1 ml of emulsified collagen II in the right hind toe. Boost immunization was induced after 7 days with an additional intradermal injection of 0.1 ml of emulsified collagen II in the right hind toe. Drug treatment began after another 7 days (14 days of CIA induction) by gavaging the rats daily with the following drug doses: COR: 3.36 g of drug per kg of body weight per day (g/kg/d); PAE: 6.27 g/kg/d; UA: 25 mg/kg/d; PF: 7.5 mg/kg/d; DEX (positive control): 0.5 mg/kg/d.

After 20 days of drug administration (34 days after initial CIA induction), the rats were sacrificed via an overdose of sodium pentobarbital. Synovial tissues were isolated by cutting open the knee joint to expose the kneecap and separating the muscles. The synovial and fibrous layers of the joint capsule were separated using surgical scalpels, and synovial tissues were extracted. The right ankle joint was fixed in 4% paraformaldehyde, decalcified, and paraffin-embedded for histological examination.

Tissue sections were heated at 65°C for 1 h and immersed in xylene for 30 min. After deparaffinization, the sections were washed in tap water for 10 min. Antigen retrieval was performed by immersing the sections in 0.01 M sodium citrate buffer at 125°C and 103 kPa for 20 min, and the sections were then washed three times with phosphate-buffered saline (PBS) after cooling. Endogenous peroxidase was performed in 3% H₂O₂ at 37°C for 10 min and the sections were thereafter washed three times with PBS. The sections were then blocked in 10% goat serum (SL038, Solarbio) at 37°C for 20 min in a humidified atmosphere and washed three times with PBS. For staining, the sections were incubated overnight at 4°C with primary antibodies against galectin-3 (LGALS3, PAB32080, Bioswamp, Wuhan, China), signal transducer and activator of transcription 3 (STAT3, PAB30641, Bioswamp), and vascular endothelial growth factor-A (VEGFA, PAB45842, Bioswamp) (all diluted at 1:50 in PBS). Thereafter, the sections were washed three times with PBS and incubated for 1 h at 37°C with secondary antibodies in the MaxVision™ HRP-Polymer anti-Mouse/Rabbit IHC Kit (KIT-5020, Maxim, Fuzhou, China). After three washes with PBS, the sections were stained with diaminobenzidine, washed with tap water, and counterstained with hematoxylin (G1140, Solarbio) for 3 min. After ethanol dehydration and xylene washes, the sections were sealed and observed under a microscope. Brown areas represent positive staining.

The intensity of positive immunohistochemical staining was quantified by ImagePro Plus using three sections from each treatment group (n = 3). The data are represented as the average integrated optical density ± standard deviation. One-way analysis of variance with Tukey’s post-hoc analysis was performed to evaluate the differences among groups. p < 0.05 was considered as statistically significant.

The molecular structures of the selected active compounds of Cornus officinalis (loganin, ursolic acid, and morroniside) and Paeonia lactiflora (paeoniflorin and albiflorin) investigated in this study are shown in Figure 1A. The “absorption, distribution, metabolism, and excretion” (ADME)-related pharmacological parameters of these compounds are listed in Table 1, and the detailed descriptions of these parameters can be found on the TCMSP web site (https://tcmspw.com/load_intro.php?id=29). The “drug-likeness” of a compound can be evaluated from these parameters by referring to Lipinski’s “Rule of 5”, which is a guide to determining whether a compound’s pharmacological activity is likely to render it orally active in humans (Lipinski et al., 2001). Among the five compounds, ursolic acid is the only one that satisfies three out of the four criteria in the Rule of 5: Hdon < 5, Hacc < 10, and MW < 500. In addition, paeoniflorin exhibits the highest OB (oral bioavailability) score of 53.87% and the highest DL (drug-likeness) score of 0.79. These results signify that among the five active compounds, ursolic acid and paeoniflorin (active compounds in Cornus officinalis and Paeonia lactiflora, respectively) exhibited the most “drug-like” pharmacological properties.

The protein targets of these compounds were predicted using three databases: Swiss, TCMSP, and PharmMapper (Figure 1B). In total, 132 unique target proteins were identified between the five compounds, with five common targets between Swiss and TCMSP. However, no target was identified in all three databases. The results of GO analysis revealed that the 132 target genes were significantly enriched in 93 GO terms, of which the top 20 terms based on adjusted p-value and gene count are shown in Figures 2A,B, respectively. The common terms found in the bar and bubble plots were identified and listed in order of gene count in Table 2. Among the eight GO terms, 52 unique target genes were identified, among which 12 (ATF2, CTSB, DDP4, FOS, IGHG1, JUN, KLK7, MMP1, MMP2, MMP3, MMP10, and PLAU) were found in three GO terms (the maximum). Similarly, the results of KEGG analysis revealed that the 132 target genes were significantly enriched in 118 signaling pathways. The top 20 KEGG terms based on adjusted p-value are shown in Figure 3A, and these terms were rearranged in order of gene count in Figure 3B. The top 10 terms in Figure 3B and their target genes are listed in Table 3. Among the ten KEGG terms, 46 unique target genes were identified, among which RELA and IKBKG were found in all ten KEGG terms. The next most enriched terms were TP53, which appeared nine times, and CASP3, IL6, and JUN, which each appeared eight times in the top ten KEGG terms.

To evaluate the relationship between the identified target proteins, a PPI network was constructed to present the direct and indirect regulatory interactions among the targets (Figure 4A). The nature of the interactions is differentiated using lines of different colors, as illustrated in Figure 4A. In addition, a visual PPI network diagram was generated to evaluate the topology parameters of the network nodes (Figure 4B). The size and darkness of the nodes represent the number of degrees of each node. The thickness and color of the edges represent the combination score, with yellow and blue indicating lower and higher combination scores, respectively. The core genes were further analyzed and ranked by the number of times that each gene appears among relationship pairs within the network diagram (Figure 4C). Among these interactions, there were 131 nodes and 913 edges in total, with an average node degree of 13.9. Nodes that are larger and darker in red are those with higher degrees of freedom (Figure 4B), suggesting that the represented genes are closely linked to other genes in the network and may be critically associated with RA pathogenesis. Among these, the top ten key genes in terms of degrees of freedom are TP53, IL6, TNF, VEGFA, STAT3, CASP3, JUN, MAPK8, IL1B, and PTGS2, with respective degrees of freedom of 62, 61, 60, 59, 53, 51, 46, 45, 42, and 42 (Figure 4C).

Next, a pharmacological network map was constructed to illustrate the interplay between compound-target relationships (Figure 5). The active compounds investigated in this study are displayed in pink hexagons and their targets are shown in yellow circles. Genes that are targeted by two or more compounds are revealed by identifying the nodes with multiple connections (Table 4). Among the identified targets, four (LGALS3, LGALS9, ADORA2A, and IGHG1) were the common targets of three compounds and 18 were the common targets of two compounds.

We then screened for the target genes of RA in four databases: OMIM, GenCLiP3, CTD, and GeneCards. Among all of the identified targets, Venn diagram analysis revealed that there were 71 common targets of RA from these databases (Figure 6A). These 71 targets were compared with the 132 previously identified unique targets of the five active compounds (loganin, ursolic acid, morroniside, paeoniflorin, and albiflorin), revealing eight common targets (Figure 6B). These are namely IL1β, VEGFA, STAT3, TP53, IL6, TNF, FOS, and LGALS3.

To elucidate the mechanism of compound-target binding, we evaluated the binding energy between the eight common targets and the five active compounds (Table 5). Lower binding energies correspond to stronger binding between compound and target (Aamir et al., 2018). Consistent with the compound-target relationships shown in Figure 5 and Table 4, the binding strengths between each identified target with their corresponding active compounds/compounds (bold values in grey cells in Table 5) are the highest among all active compounds. The exception was LGALS3, which was predicted to be targeted by morroniside, paeoniflorin, and albiflorin (Figure 5) but showed the greatest binding strength with ursolic acid, morroniside, and paeoniflorin. Using AutoDock Vina, we performed a molecular docking analysis between the five active compounds and three selected targets among the eight identified: LGALS3, STAT3, and VEGFA. The PDB IDs of the receptor proteins of LGALS3, STAT3, and VEGFA are 4BL1, 6NJS, and 1KAT, respectively. The results of molecular docking visualization (Figure 7) showed that the three targets were able to spontaneously bind to each of the five active compounds, via forces such as hydrogen bonds to form a stable conformation.

For a preliminary in vivo experimental validation of the predicted compound-target relationship, we evaluated the expression of LGALS3, STAT3, and VEGFA in a rat model of CIA (Figure 8A). Successful establishment of the CIA model was confirmed by assessing the arthritic score and swelling degree in the rats throughout the experimental period, as well as histological verification (data not shown). After the induction of CIA, Sprague-Dawley rats were subjected to daily treatment of Cornus officinalis (COR) or/and Paeonia lactiflora (PAE), their respective active compounds ursolic acid (UA) or/and paeoniflorin (PF), or dexamethasone (DEX) as a positive control (Cohen et al., 1960). UA and PF were selected for these experiments because they had previously been identified as the compounds with the most drug-like pharmacological properties. After 20 days of treatment, synovial tissues were isolated from the experimental rats and the synovial expression of LGALS3, STAT3, and VEGFA was assessed by immunohistochemistry (Figure 8B). Upon CIA induction, we noticed a significant increase in the positive synovial expression of LGALS3, STAT3, and VEGFA compared to that in the Control group (p < 0.05). The administration of COR or/and PAE, UA or/and PF, or DEX significantly reduced the expression of the three proteins in the synovial tissues of CIA-induced rats. In terms of VEGFA expression, the difference between the treatments was not significant. However, the effect of PF + UA and DEX was significantly better than that of all other treatments (except COR alone) in downregulating LGALS3 expression. Furthermore, PF + UA exerted a stronger effect in downregulating STAT3 expression than PAE and COR + PAE. Collectively, these results signify that the expression of LGALS, STAT3, and VEGFA, which can be considered as tissue biomarkers of RA, act as indicators of the effectiveness of treatment using various compounds.