A flowchart of the overall experimental design is shown in Fig. 1. The mRNA datasets were downloaded from the GEO database (www.ncbi.nlm.nih.gov/geo/). To gain an improved understanding of the underlying mechanism of RA, the independent datasets of GSE55457(6) and GSE77298(7) were selected, which included synovial tissues from 13 and 16 patients with RA and 10 and 7 healthy subjects, respectively.

To identify differentially expressed mRNAs (DE-mRNAs) using the GSE55457 and GSE77298 datasets, the mRNA expression levels of the tissue samples from patients with RA and healthy subjects were compared. DE-mRNAs were defined using the ‘limma’ package (version 3.52.2) (9), based on the criteria of log₂ fold change >1 and P<0.05(8).

Database for Annotation, Visualization and Integrated Discover (DAVID) (https://david.ncifcrf.gov/) analysis was performed on Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) terms (10). In the present study, the ‘enrichGO’ function of the clusterProfiler package (10) was used to perform GO pathway enrichment analysis based on the differentially expressed genes (DEGs). Based on their function, these genes were divided into three categories: Biological process (BP), molecular function (MF) and cellular component (CC). A significant enrichment of GO terms was observed for DEGs based on the false discovery rate-adjusted P<0.05. The DAVID webtools were used to analyze KEGG enrichment of the integrated DEGs in the present study. Adjusted P<0.05 was considered to indicate a statistically significant difference. The software packages R (version 4.2.2; https://www.R-project.org/) (11), MCODE plug-in (version 2.0.2) (12,13) and CytoHubba plug-in (https://apps.cytoscape.org/apps/cytohubba) in Cytoscape (version 3.9.1) (14,15) were used to obtain the results.

Based on the STRING database (version 11.5; using confidence score >0.9) (16), a set of high-confidence PPI networks for humans were constructed. To screen for hub genes, topological algorithms of the cytoHubba plug-in application in Cytoscape (version 3.9.1) were applied to predict and explore important nodes and hubs. The core module was extracted from the PPI network using the MCODE plug-in. The parameters were set as follows: Degree cut-off=2, K-Core=2 and node score cut-off=0.2(17). It was hypothesized that the potential hub genes could be found in the core module. Thus, the genes selected were considered the hub genes.

A total of 30 male Sprague Dawley rats (6 weeks old; 180±20 g body weight) were obtained from Chengdu Dashuo Laboratory Animal Co., Ltd. (permit no. SCXK 2015-030). The animals were fed and maintained under laboratory conditions of 22-24˚C and 50-65% humidity, with natural light and dark cycles, and had free access to food and drinking water. The rats were randomized into the control and RA groups (n=15/group). All animal procedures were approved by The Ethics Committee of Chengdu University of Traditional Chinese Medicine (Chengdu, China; approval no. 2018-11).

On the first day of the experiment, to establish the RA model, the rats were in the RA group subjected to the induction process, wherein they received an injection of 0.5 ml/kg of FCA (Sigma-Aldrich; Merck KGaA) in the footpad. The control group was injected with the same amount of normal saline. RA progression was evaluated based on the swelling and redness of the foot pad. The overall condition and bipedal lesions of all rats were observed daily. A waterproof marker was used to mark the middle of the footpad of each rat and all future measurements were taken from this point. Calipers were used to measure the footpad thickness before and 1, 2 and 3 weeks after establishment of the RA model (18).

At the end of the experiment, the rats were anesthetized with 7% chloral hydrate (400 mg/kg, intraperitoneally; Sinopharm Chemical Reagent Co., Ltd.). It was confirmed that the rats showed no signs of peritonitis, pain or discomfort following this anesthesia method. Blood was subsequently collected from the abdominal aorta of the live anesthetized rats. After the 5 ml of blood was collected, the rats were subsequently sacrificed by cervical dislocation. Ankle joints, synovial tissue and the spleen were also collected from each animal.

The ankle joint was fixed with 4% paraformaldehyde (Wuhan Servicebio Technology Co., Ltd.) for >24 h at room temperature. After fixation, the paraffin-embedded tissue sections were sliced into 5-mm sections using a paraffin microtome (Leica Microsystems, Ltd.). The sections were then stained with the H&E staining solution kit (Wuhan Servicebio Technology Co., Ltd.), which included staining with hematoxylin and eosin individually (5 min at room temperature for each staining step). The stained sections were dehydrated with ethanol, cleared with xylene, mounted onto slides with neutral balsam and covered with coverslips. Histopathological examination was performed on the mounted sections and images were captured under a light microscope (Nikon Corporation). The slides were graded according to the degree of inflammatory changes. The rat synovial pathological score included synovial tissue hyperplasia, inflammatory cell infiltration and macrophage hyperplasia. The tissue was scored based on the following: 0, normal; 1, sparse or scattered; 2, diffuse increase; 3, large increase with respect to inflammatory cell infiltration, macrophage hyperplasia and fibroplasia.

The blood samples collected from each group were centrifuged at 1,200 x g for 20 min at room temperature (Heal Force Bio-meditech Holdings Ltd.) and the supernatant was collected. The expression levels of tumor necrosis factor α (TNF-α) (cat. no. JM-01587R1) and interleukin (IL-) 12 (cat. no. JM10372R1) in the supernatant were detected using an ELISA kit (all Jiangsu Jingmei Biotechnology Co., Ltd.), according to the manufacturer's guidelines. Standards (50 µl) and samples (10 µl of serum and 40 µl of diluent, supplied with the kit) were added to predetermined microplate wells. Horseradish peroxidase (HRP)-labeled conjugate (100 µl; supplied with the kit) was added to the wells, which were then sealed and incubated at 37˚C for 1 h. The plate was washed five times with PBS and substrates A (50 µl) and B (50 µl) were added to each well. After incubation for 15 min at 37˚C, stop solution (50 µl) was added. A microplate reader (BioTek Instruments, Inc.) was used to measure the optical density (OD) values in all wells at a wavelength of 450 nm. TNF-α and IL-12 concentrations were calculated using the obtained OD values of each group of samples, according to the linear regression equation of the standard curve.

Spleen tissue (100 mg) from the RA animal model and control group was homogenized using a homogenization tube. Following the manufacturer's instructions, total RNA was extracted using the TRIzol™ method (Thermo Fisher Scientific, Inc.). The concentration and purity of the isolated RNA were determined using a microvolume spectrophotometer (Thermo Fisher Scientific, Inc.). According to the manufacturer's protocol, RNA was subsequently converted into cDNA using a RT kit (Thermo Fisher Scientific, Inc.) and a PCR machine (Eastwin Life Sciences Inc.). For the qPCR step, the reaction mixture constituted the following components: 2 µl cDNA, 2.5 µl forward and reverse primers, 7.5 µl 2X SYBR Green qPCR Master Mix (Wuhan Servicebio Technology Co., Ltd.) and ultrapure water to reach a 15 µl total volume. Amplification was performed using a fluorescent qPCR instrument (Bio-Rad Laboratories, Inc.). The reaction parameters were as follows: Initial denaturation (95˚C for 30 sec), 40 cycles of denaturation (95˚C for 15 sec) and annealing/extension (60˚C for 30 sec). Data were collected and analyzed using SDS software (version 1.4; Applied Biosystems; Thermo Fisher Scientific, Inc.). GAPDH was used as the internal reference gene. Furthermore, 2^-ΔΔCq (19) was used to calculate the relative expression level of discs large homolog-associated protein 5 (DLGAP5), kinesin family member 20A (KIF20A), maternal embryonic leucine zipper kinase (MELK) and nuclear division cycle 80 (NDC80). The primer sequences used were as follows: GAPDH-forward (F), 5'-CTGGAGAAACCTGCCAAGTATG-3'; GAPDH-reverse (R), 5'-GGTGGAAGAATGGGAGTTGCT-3'; MELK-F, 5'-TCCCATTGAGTGGCAAAGCA-3'; MELK-R, 5'-CCTGCTGTTGCGGTGATGTA-3'; NDC80-F, 5'-CAAAGACCTGGAAGCCGAACAA-3'; NDC80-R, 5'-CTAGCCAACTTGTGGTACTCTGC-3'; DLGAP5-F, 5'-GAAGTGCCAACGCTTTCCTGA-3'; DLGAP5-R, 5'-CGTTCTTTCCTAGAGTCTGTTTGTG-3'; KIF20A-F, 5'-GAACCACCTAGCCATCAGCATA-3'; KIF20A-R, 5'-ACTGTGACTGCGACTAGAGTGCT-3'.

The synovial tissue was washed 2-3 times with cold PBS to remove blood, after which the tissue was cut into small pieces and placed in a homogenizer tube. Protease inhibitors and 10 times the tissue volume (10 ml buffer/g tissue) of radio immunoprecipitation assay lysis buffer (Wuhan Servicebio Technology, Co., Ltd.) were added to the tube and the contents were homogenized thoroughly on ice. The homogenate was transferred to a 1.5 ml centrifuge tube and the contents were shaken well. The homogenate was pipetted repeatedly for 30 min in an ice bath to ensure complete cell lysis, after which the contents were centrifuged at 12,000 x g for 10 min at 4˚C and the supernatant (containing total protein) was collected. Protein concentration was determined using the bicinchoninic acid method. Subsequently, sodium dodecyl sulfate-polyacrylamide gel electrophoresis was performed using 10% gels to separate the isolated proteins (30 µg/lane). The separated proteins were then transferred onto polyvinylidene difluoride membranes. After blocking the membranes with tris-buffered saline containing 5% skimmed milk powder at room temperature for 1 h, the membranes were incubated with the following primary antibodies overnight at 4˚C: KIF20A (cat. no. ab70791; 1:1,000; Abcam), NDC80 (cat. no. ab109496; 1:1,000; Abcam), MELK (cat. no. ab108529; 1:1,000; Abcam), DLGAP5 (cat. no. ab70744; 1:1,000; Abcam) and GAPDH (cat. no. gb15002; 1:2,000; Wuhan Servicebio Technology, Co., Ltd.). After repeated washing, the membranes were incubated with HRP-conjugated goat anti-rabbit (cat. no. gb23303; 1:5,000; Wuhan Servicebio Technology, Co., Ltd.) secondary antibody for 2 h at 37˚C. Images were collected using a chemiluminescent imaging system (Wuhan Servicebio Technology, Co., Ltd.) and grayscale analysis was performed using the system's ChemiScope analysis software (version ChemiScope 6100; Clinx Science Instrument Co., Ltd.).

All data had to pass the Shapiro-Wilk normality test and the Levene's test for homogeneity of variance. Then the results were analyzed by pairwise test. For correlation analyses, when the data distribution was normal, the Pearson correlation analysis was used, whereas when the distribution was not normal, the Spearman correlation analysis was used. Data were presented as the mean ± standard deviation. P<0.05 was considered to indicate a statistically significant difference. All data were assessed using SPSS 26.0 software (IBM Corp.). The relevant boxplots were compiled using Xiantao Academic (http://www.xiantao.love).

The mRNA expression levels from the GSE55457 and GSE77298 datasets from the GEO database were investigated and the results were determined through boxplots, volcano plots and heatmaps (Fig. 2). There were 13 samples from patients with RA and 10 samples from healthy individuals in the GSE55457 dataset, whereas there were 16 samples from patients with RA and 7 samples from healthy individuals in the GSE77298 dataset. In the GSE55457 dataset, 664 DE-mRNAs were identified (429 with upregulated expression and 235 with downregulated expression), while in the GSE77298 dataset, 2,150 DE-mRNAs were detected (1,170 with upregulated expression and 980 with downregulated expression). According to Venn diagram analyses, 165 DE-mRNAs that had upregulated expression and 30 DE-mRNAs that had downregulated expression were shared between the two datasets (Fig. 3 and Table SI).

To explore the potential functions of the identified DE-mRNAs, functional enrichment analysis was performed using GO and KEGG analyses. Among them, the top five BP, top five CC and top five MF GO terms are shown in Fig. 4A and Table SII. In BP, ‘leukocyte cell-cell adhesion,’ ‘mononuclear cell differentiation,’ ‘lymphocyte differentiation,’ ‘leukocyte migration’ and ‘lymphocyte mediated immunity’ were the top five most enriched terms. Additionally, ‘external side of plasma membrane,’ ‘endocytic vesicle membrane,’ ‘clathrin-coated endocytic vesicle,’ ‘clathrin-coated endocytic vesicle membrane’ and ‘alpha-β T-cell receptor complex’ were the top five GO terms on CC. Furthermore, ‘antigen binding,’ ‘immune receptor activity,’ ‘cytokine receptor activity,’ ‘chemokine receptor binding’ and ‘chemokine activity’ were the top five GO terms on MF. Based on the KEGG analysis, significant enrichments were observed in nine KEGG pathways for the apparent DE-mRNAs, such as ‘hematopoietic cell lineage,’ ‘cytokine-cytokine receptor interaction,’ ‘viral protein interaction with cytokine and cytokine receptor,’ ‘chemokine signaling pathway’ and ‘Th17 cell differentiation’ (Fig. 4B and Table SII).

To further explore the functions of these mRNAs, the STRING database was analyzed using a PPI network analysis and mRNAs with the highest confidence levels were selected. Overall, 195 DE-mRNAs were filtered into a PPI network complex containing 181 nodes and 781 edges. Then, MCODE was used to identify the most connected hub mRNAs in the PPI network (Fig. 5A and B). Then, cytoHubba was used to analyze the top 10 most interacting central mRNAs: Assembly factor for spindle microtubules (ASPM), DLGAP5, DNA topoisomerase IIα (TOP2A), KIF20A, BUB1, KIF11, CDC20, MELK, NDC80 and ribonucleotide reductase regulatory subunit M2 (RRM2) (Fig. 5C).

A schematic diagram of the flowchart after the establishment of the FCA-induced arthritis rat model is shown in Fig. S1A. The rat footpad thickness was measured and the thickness in the model group was significantly increased compared with the control group (P<0.01; Fig. S1B). H&E staining of the joint tissue demonstrated that, compared with the control group, a greater amount of inflammatory cell infiltration, macrophage hyperplasia and fibrous tissue hyperplasia was present in the model group (Fig. S1C). The groups were then scored blindly with all criteria tested by at least two investigators and the average scores were obtained. The scores represent the degree of tissue damage with a higher score indicating a higher degree of tissue damage. Compared with the control group, the joints of the rats in the model group showed apparent inflammatory cell infiltration, fibrous tissue hyperplasia and macrophage proliferation (P<0.01; Fig. S1D-F). Therefore, the results demonstrated that, compared with the control group, the RA group had a higher degree of bone tissue destruction and inflammation, indicating that the FCA-induced RA rat model was successfully established.

The mRNA expression levels of four previously unreported genes (DLGAP5, KIF20A, MELK and NDC80) among the top 10 central DE-mRNAs (or hub genes) from the PPI network analysis were examined in RA rat spleen tissues using qPCR, to validate the reliability of the dataset. Compared with the control group, the expression levels of DLGAP5, KIF20A, MELK and NDC80 in the RA group were significantly increased (P<0.01; Fig. 6A-D). In addition, western blotting was performed to detect the protein expression levels of DLGAP5, KIF20A, MELK and NDC80. Compared with the control group, the protein expression levels of DLGAP5, KIF20A, MELK and NDC80 in the RA group were significantly increased (P<0.01; Fig. S2). Furthermore, ELISA results demonstrated that, compared with the control group, IL-12 and TNF-α levels in the RA group were significantly increased (P<0.01; Fig. 6E and F). In addition, correlation analysis was performed to analyze the correlation between inflammatory factors and hub genes, and the results demonstrated that the inflammatory factors were highly correlated with the hub genes (Fig. S3).