GSE55235 and GSE84074 microarray data were downloaded from the GEO database (ncbi.nlm.nih.gov/geo/) (14). The GSE55235 chip data was based on GPL96, [hg-u133a] Affymetrix Human Genome U133A Array (15). The Affymetrix, calif. GSE84074 chip was based on the GPL19640, agilent-062918 Human lncRNA array version 4.0 (Probe name version) microarray data in Bioconductor (version 1.46.1; bioconductor.org/) (16). The preprocessing stage included background correction, normalization and expression calculation. The Bioconductor Annotation data package was used to convert the chip data probe into gene symbols. If multiple probes were mapped to a gene symbol, the average was set for the final expression value of the gene. The DEGs of each group were analyzed using the limma package (17). The DEGs were analyzed using the volcano map script in the R studio (version 1.2.5042) (18,19). P<0.05 and |log(fold change) FC|>2 were used as the cutoff criteria. The Heatmap package was used for hierarchical clustering analysis and visualization of DEGs (20).

The Clue-GO plugin in Cytoscape (21) was used to analyze the functions and pathways of the DEGs. P<0.05 was used as the cutoff criteria. Enrichment analysis was performed and visualized using ClueGO and CluePedia with P<0.05 as the cutoff criterion (22-24). The function and pathway enrichment analysis of DEGs following LLDT-8 processing was performed using the clusterProfiler package (25).

Synoviocytes were purchased from Shanghai Institute of Biochemistry and Cell Biology. Synovial cells were grown in DMEM supplemented with 10% FBS (Gibco; Thermo Fisher Scientific, Inc.) in a humidified incubator at 37˚C with 5% CO₂.

Cells in the logarithmic growth phase were used for transfection. A total of 4x10⁵ cells/well were plated in a 6-well plate and 2 ml supplemented media without antibiotics was added. Transfection was performed when the confluence of cells reached 70%. For transfection, the negative control small interfering (si)RNA was used, and for the experimental group si-CD2 was added. Transfections were performed using Lipofectamine™ 2000 (Thermo Fisher Scientific, Inc.). The final siRNA concentration used was 40 nmol/l. The sequences of the siRNAs were: si-NC, 5'-UUCUCCGAACGUGUCACG UTT-3' and 5'-ACGUGACACGUUCGGAGAATT-3'; and for si-CD2, 5'-CCAAAGGUGCAGUCUCCAATT-3' and 5'-UUG GAGACUGCACCUUUGGTT-3'. The siRNAs were purchased from Santa Cruz Biotechnology, Inc. Follow-up experiments were performed 48 h after transfection.

Cell proliferation was detected using a Cell Counting Kit-8 (CCK-8) assay. A total of 1x10³ cells in 200 µl were added per well to a 96-well plate. Following the different treatments, 10 µl CCK-8 solution (Dojindo Molecular Technologies, Inc.) was added to each well, and the plate was further incubated for 2 h in a 5% CO₂ incubator. The absorbance of cells was measured using a microplate reader at 450 nm. All experiments were repeated 3 times.

Cells were seeded in 6-well plates containing 2 ml complete medium at a density of 1x10⁶. The 6-well plates were incubated at 37˚C, 5% CO₂ and 100% humidity for 24 h. TNF-α (10 ng/ml; cat. no. 14-8329-62; Thermo Fisher Scientific, Inc.) or IL-1β (10 mg/l; cat. no. SRP3083; Sigma-Aldrich; Merck KGaA) was added to stimulate the cells for 24 h. Total RNA was extracted using TRIzol^® (cat. no. 15596026; Thermo Fisher Scientific, Inc.). Isolated RNA was reverse transcribed into cDNA using a reverse transcription kit according to the manufacturer's protocol (cat. no. RR036A; Takara Bio, Inc.). The sequences of the primers were: CD2 forward, 5'-CCCATGATTCCTTCATATTTGCA-3' and reverse, 5'- GGTCCTTCTCCAGCCTAGT-3'; IL-6 forward, 5'-CCACTCACCTCTTCAGAACG-3' and reverse, 5'-CATCTT TGGAAGGTTCAGGTTG-3'; LDH forward, 5'-GTTGGAACT GGTGCCGTAG-3' and reverse, 5'-GAAAAGACTGCCATG CTGAAG-3'; IL-1β forward, 5'-CCTGTCCTGCGTGTTGAA AGA-3' and reverse, 5'-GGGAACTGGGCAGACTCAAA-3'; and β-actin forward, 5'-GATCTTGATCTTCATGGTGCTAG-3' and reverse, 5'-TTGTAACCACCTGGGACGATATGG-3'. β-actin was used as the internal reference. For qPCR, 10 µl cDNA was used per sample, and the thermocycling conditions were: 95˚C for 5 min; followed by 40 cycles of 95˚C for 10 sec and 60˚C for 30 sec. The 2^-ΔΔCq method was used for relative quantitative analysis (26).

LDH and IL-6 ELISA kits (Beijing Solarbio Science & Technology Co., Ltd.) were used to detect the levels of inflammatory factors in cells following the different treatments, according to the manufacturer's protocol.

The CentiScape plugin (27) was used to calculate the node degree, and the molecular complex detection (MCODE) plugin (28) was employed to find clusters in the entire PPI network. Nodal proteins may have important physiological regulatory functions and may be key candidate genes. The genes in the most significant modules were extracted for functional enrichment analysis. P<0.05 was set as the cutoff value.

The online database STRING (string-db.org) was used to construct the PPI network of DEGs (29). Cytoscape was used to construct the protein interaction network and analyze the interactive relationship of the candidate DEGs encoding proteins in RA. Hub genes were determined on the basis of the degree results.

The protein structure of CD2 was downloaded from the PDB database (rcsb.org; PDB ID, 1HNF) (30). Schrödinger (Maestro-2015-2) software was used to optimize the structure of the protein, which included assigning bond sequences, hydrogenating atoms, and removing eclipsed conformations. Finally, the energy of CD2 was minimized. The Schrodinger LigPred module was used to optimize the energy of the small molecule structure (31). First, 651 Chinese medicine molecules with the highest CD2 docking scores were screened using high throughput virtual (HTV) screening. The standard precision (SP) screen was selected from 34 Chinese medicine molecules. Finally, three candidate Chinese medicine molecules (extra precision screen) were obtained via precise docking. PyMOL (version 2.3) (32) was used to visualize the structures of triptolide derivative (5R)-5-hydroxytriptolide (LLDT-8) and CD2.

Genes that may serve a key role in RA were determined through the statistical calculation on the basis of GSEA (version 3.0) (33), which was used to compare the genes to be analyzed with those that were pre-divided into gene sets. The input file in the GCT/CLS format was created in accordance with the format requirements of the software. The CLS file was the grouping description file of the expression matrix of sequencing data, and the GCT file was the gene expression matrix files in sequencing. For data analysis, the parameters were set as: Number of permutations, 1,000; Collapse dataset to gene symbols, true; Permutation type, gene_set; Max size, exclude larger sets, 500; Min size, exclude smaller sets: 15.

SPSS version 20.0 (IBM, Corp.) was used for analysis, and measurement data are expressed as the mean ± standard deviation. Comparisons between two groups were performed using a Student's t-test, and comparisons between multiple groups were performed using a one-way ANOVA followed by a post-hoc Tukey's test. P<0.05 was considered to indicate a statistically significant difference.

Data were downloaded from NCBI-GEO, a public functional genomics database that receives data from different microarray platforms, and the gene expression profiles of RA and normal synovial tissues were obtained from the GSE55235 dataset (15). The data included synovial tissues from 10 patients with RA and 10 individuals with healthy joints. P<0.05 and |log(fold change) FC|>2 were used as the cutoff criteria. A total of 279 DEGs (166 upregulated and 113 downregulated genes) were extracted from the expression profile dataset (Fig. 1A). The DEGs are shown as a heatmap based on the |logFC| value in Fig. 1B. According to the results, the DEGs identified could be used to distinguish RA from the normal group and could be used for subsequent biological analysis. The GSEA software package was used to analyze the different biological processes of synovial tissues in patients with RA and healthy individuals, and the truncated value was P<0.05. GSEA is a computational method used to determine whether a predefined set of genes can show significant consistent differences in two biological states. In the GO analysis, the upregulated functions associated with the significantly enriched genes included the activation of immune and adaptive immune responses. The connective tissue development, cell number maintenance, and cellular response to extracellular stimulus function decreased (Fig. 2A). KEGG enrichment analysis showed that the B cell receptor, chemokine, the JAK-STAT signaling pathways and oxidative phosphorylation were activated in RA lesions (Fig. 2B). Identification of significantly rich functions and pathways may assist in further studying the role of DEGs in RA.

The study of protein interaction networks assist in deepening our understanding of cell structure and function from a systematic perspective and provide a theoretical basis for the discovery of novel drug targets and drug design. Using bioinformatics methods to study protein interaction network shows great advantages, including protein interaction network mapping and display, network topology analysis, network structure module research and network comparison. Based on the obtained DEGs, the interaction network of DEGs was further analyzed. The STRING database was used to calculate the PPI network (Fig. S1). As observed from the results, the PPI network can be divided into five major interaction subset networks.

The key nodes in the PPI network were analyzed using Cytoscape plugins, to further analyze the hub genes in RA and determine the association of the protein interactions. The results showed that the interactive network had 41 nodes, 135 edges and a MCODE score of 6.750 (Fig. 3A). The degree of each node gene was further calculated, and CD2 had the highest topological connectivity (Fig. 3B). Furthermore, the Clue-GO plugin in Cytoscape software was used to enrich and annotate the key PPI module genes. Results showed that the key module genes were enriched primarily in response to corticosterone, antigen processing, and presentation of peptide or polysaccharide antigen via MHC class II, lymphocyte co-stimulation, interferon-γ-mediated signaling pathway, T cell co-stimulation and the T cell receptor signaling pathway (Fig. 3C).

According to the results of bioinformatics analysis, CD2 was determined to be a key protein involved in the development of RA. In order to verify the role of CD2 in RA, IL-1 and TNF-α were used to induce synovial cells for establishment of an in vitro RA cell model. The experimental results showed that IL-1 and TNF-α could induce the upregulation of CD2, consistent with the results of the bioinformatics analysis (Fig. 4A). Next, siRNAs were used to knockdown the expression of CD2. The experimental results show that the siRNAs effectively reduced the expression of CD2 (Fig. 4B). The results of cell proliferation experiments showed that compared with the control group, CD2 knockdown reduced the inhibitory effect of TNF-α on proliferation in synovial cells (Fig. 4C). Furthermore, changes in the levels of cellular inflammatory factors were assessed. The results showed that compared with the control group, TNF-α stimulation increased the levels of inflammatory factors, whereas knockdown of CD2 reduced the effects of TNF-α (Fig. 4D-F).

By analyzing the PPI network, CD2 was found to be a key node gene that serves an important role in the pathogenesis of RA. CD2 was significantly upregulated in RA (logFC=2.731) and therefore was chosen as the basis for the selection of candidate compounds. A set of molecular screening strategies for TCM was designed on the basis of high-throughput virtual screening method. The schematic diagram of the screening process is shown in Fig. 5A. First, 651 Chinese medicine molecules with the highest CD2 docking scores were screened using a HTV screening. The SP screen was selected from 34 Chinese medicine molecules. Finally, three candidate Chinese medicine molecules (extra precision screen) were obtained via precise docking. Amongst these molecules, LLDT-8 and CD2 were found to have the highest docking score and the best binding mode (Fig. 5B). The binding between LLDT-8 and the target protein CD2 is shown in Fig. 5C. The key amino acid residues (Gln-15, Leu-13 and Arg-105) interacted with LLDT-8 to form a hydrogen bond interaction. The 3D conformation of LLDT-8 interacting with CD2 was used to further clarify their interaction. LLDT-8 may inhibit the function of CD2 by inhibiting its biological activity center, and this inhibition may be beneficial for the management of RA.

A total of 324 DEGs were identified following pretreatment and batch elimination of DEGs in the GSE84074 dataset using the limma software package (16). Amongst these genes, 139 upregulated and 311 downregulated genes were identified in the LLDT-8 treatment group compared with the control group (Fig. 6A). DEGs are shown in the volcano plot, and the first 100 DEGs shown based on the value of |logFC| are also illustrated in the heatmap (Fig. 6B).

The function and pathway enrichment analysis of DEGs after LLDT-8 processing was performed using the clusterProfiler package (25). DEGs were divided into three functional groups, namely, biological process, molecular function and cellular component (Figs. 6C-D and 7A). LLDT-8 primarily inhibited the cellular response to IL-1, the positive regulation of ERK1 and ERK2 cascades, the chemokine-mediated signaling pathway, inflammatory and immune responses, the positive regulation of IL-6 production function and chemokine signaling pathway, innate and adaptive immune systems, and the JAK-STAT signaling pathway in the RA process (Fig. 7B).

Based on high-throughput molecular screening results, it was found that LLDT-8 inhibited CD2. Therefore, the changes in the levels of CD2 were detected in the LLDT-8 treated synovial cells. The experimental results showed that compared with the control group, LLDT-8 reduced the expression of CD2 (Fig. 8A). The results of cell proliferation experiments showed that TNF-α could inhibit the proliferation of synovial cells, whereas LLDT-8 treatment restored the proliferation of synovial cells (Fig. 8B).

Further, changes in cellular inflammatory factors were detected. The experimental results showed that, compared with the control group, TNF-α stimulation increased the levels of inflammatory factors, whereas knockdown of CD2 reduced the stimulation of TNF-α. After knocking down CD2, the administration of LLDT-8 did not further reduce the levels of inflammatory factors, suggesting that CD2 was the target of LLDT-8 (Fig. 8C and D).