Construction and Analysis of ceRNA Regulatory Networks Reveal the Core Genes Associated with Rheumatoid Arthritis

Introduction

Rheumatoid arthritis (RA) is a chronic systemic autoimmune disease characterized by persistent synovitis, progressive cartilage damage, bone erosion, pain, and disability.^1–3 Despite major therapeutic advances, RA remains biologically heterogeneous, and substantial unmet need persists in early diagnosis, patient stratification, and molecular monitoring.^2,3 Because peripheral blood can be sampled repeatedly and captures systemic immune activation, blood-derived molecular signatures are of particular interest for translational biomarker studies in RA.

Non-coding RNAs have emerged as important regulators of inflammatory signaling, stromal activation, and immune-cell differentiation in arthritic disease.⁴ Among them, circRNAs are covalently closed RNA molecules generated by back-splicing and are relatively stable in cells and body fluids.⁵ One well-studied conceptual framework is the competing endogenous RNA (ceRNA) model, in which RNAs containing shared microRNA-response elements can influence each other through competition for miRNA binding.⁶ In RA, dysregulated circRNAs and miRNAs have been reported in synovial tissue, fibroblast-like synoviocytes, and peripheral blood mononuclear cells (PBMCs), supporting their relevance to disease biology and biomarker discovery.^4,7–12 For example, circPTPN22 has been proposed as a PBMC biomarker in RA,¹² and the circRNA landscape of RA PBMCs has been systematically profiled in GSE189338.¹³

Most previous network studies in RA have focused on circRNA-miRNA-mRNA relationships alone. However, transcription factors form an additional regulatory layer that can amplify, constrain, or redirect post-transcriptional effects. Incorporating TFs into a circRNA-miRNA framework may therefore yield a more biologically informative model of gene regulation. At the same time, predicted ceRNA interactions remain correlative and should be interpreted cautiously, particularly when they are derived from independent public datasets generated on different platforms and in different cohorts.

In the present study, we integrated circRNA, miRNA, and mRNA expression profiles from three RA blood-based GEO datasets and combined them with TF information from TRRUST to construct a circRNA-miRNA-TF-mRNA network. We then used RT-qPCR in a collagen-induced arthritis (CIA) rat model to validate selected candidates. The aim was not to claim a definitive mechanism, but to identify a compact, biologically plausible, and experimentally supported regulatory framework for further study.

Materials and Methods

Ethics and Study Design

All animal procedures were approved by the Laboratory Animal Ethics Committee of Fujian University of Traditional Chinese Medicine (Approval No. FJTCM IACUC2021085) and were carried out in accordance with the Guide for the Care and Use of Laboratory Animals and institutional requirements. Reporting of the animal experiments follows the ARRIVE 2.0 recommendations.¹⁴

The bioinformatics component used only publicly available, de-identified human datasets from GEO and involved no contact with human participants and no attempt at re-identification. According to Article 32, items 1 and 2, of the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects (China; effective February 18, 2023), secondary analysis of publicly available, de-identified data may be exempt from local ethical review. Fujian University of Traditional Chinese Medicine does not have a separate IRB procedure for exempt secondary analyses of public data; this part of the study was therefore exempt from ethical review under the above national legislation.

Public Datasets and Preprocessing

Three RA-related GEO datasets were analyzed: GSE189338 for circRNA expression in PBMCs from 4 patients with RA and 4 healthy controls; GSE124373 for miRNA expression in PBMCs from 28 patients with RA and 18 healthy controls; and GSE17755 for mRNA expression in peripheral blood cells from 112 patients with RA and healthy comparators.^13,15 GEO is a public archive for high-throughput functional genomics data.¹⁵ Probe identifiers were converted to gene symbols according to the corresponding platform annotation files. Expression values were log2-transformed when necessary and normalized with the normalizeBetweenArrays function of the limma package in R.¹⁶ The characteristics of the included GEO datasets and the numbers of differentially expressed candidates are summarized in Table 1.

Table 1 Summary of Datasets and Numbers of Differentially Expressed Candidates

Differential Expression Analysis

Differential expression between RA and control samples was assessed with the limma package.¹⁶ In keeping with the exploratory nature of the original analysis, candidates were screened using |log2 fold change| > 0.2 and nominal P < 0.05. Because the circRNA dataset was small, these features were not interpreted as stand-alone biomarkers; instead, they were carried forward only if they were supported by network integration and subsequent experimental validation.

Functional Enrichment Analysis

Differentially expressed mRNAs were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using clusterProfiler.¹⁷ Terms with P < 0.05 were considered significantly enriched. The results were used to determine whether the transcriptomic signature converged on pathways relevant to inflammation, leukocyte trafficking, and immune activation in RA.

Construction of the TF-miRNA-mRNA-circRNA Network

TF-mRNA regulatory relationships were retrieved from TRRUST v2.¹⁸ Experimentally validated miRNA-target interactions were obtained from miRTarBase 2022.¹⁹ Candidate circRNA-miRNA interactions were predicted with CircBank.²⁰ First, differentially expressed mRNAs were intersected with the TRRUST output to define an RA-related TF-mRNA layer. Second, differentially expressed miRNAs were linked to TFs and mRNAs through miRTarBase-supported interactions. Third, differentially expressed circRNAs were matched to the miRNA layer through CircBank predictions. The final integrated network was visualized in Cytoscape.

For biological validation, we prioritized candidates that satisfied at least two of the following criteria: (i) presence in the integrated network, (ii) topological centrality or bridge position between regulatory layers, and (iii) prior RA relevance in the Comparative Toxicogenomics Database (CTD).²¹ On this basis, 6 circRNAs, 4 miRNAs, 3 TFs, and 8 mRNAs were selected for RT-qPCR. The prioritized targets selected for RT-qPCR validation and the rationale for their inclusion are summarized in Table 2.

Table 2 Targets Prioritized for RT-qPCR Validation and Their Rationale

CTD-Based Prioritization

The relationship between candidate genes and RA was cross-checked in the Comparative Toxicogenomics Database (CTD), an expert-curated resource that integrates gene-disease and chemical-gene-disease evidence.²¹ Genes with stronger prior RA annotation in CTD were considered higher-priority candidates for discussion and downstream interpretation.

Animals and CIA Model

Twelve male Sprague-Dawley rats (180–220 g) were used for in vivo validation. Six rats were assigned to the control group and six to the CIA group. The animal experiment was planned as exploratory validation, and no formal a priori power calculation was performed.

CIA was induced by intradermal injection of 100 μg chick type II collagen emulsified in complete Freund’s adjuvant at the base of the tail, followed 7 days later by booster immunization with 50 μg chick type II collagen in incomplete Freund’s adjuvant. Control animals received saline injections. Animals were monitored daily and euthanized on day 21 after the first immunization. Blood was collected under isoflurane anesthesia, and death was confirmed by cervical dislocation.

Clinical Scoring, Histology, and Micro-CT

Arthritis severity was scored using a standard 0–4 scale for each paw, yielding a maximum arthritis index of 16 per animal. Paw thickness was recorded during model establishment. For histology, ankle or knee joints were fixed, decalcified in EDTA, embedded in paraffin, sectioned at 5 μm, and stained with hematoxylin and eosin. Histological assessment focused on synovial hyperplasia, inflammatory-cell infiltration, cartilage injury, and bone damage. For micro-CT, formalin-fixed ankle joints were scanned to evaluate bone erosion and trabecular destruction.

RNA Extraction and RT-qPCR

Peripheral blood obtained in anticoagulant tubes was processed on ice. Total RNA was extracted from 1.5 mL whole blood using RNA isolater Total RNA Extraction Reagent (Vazyme, Nanjing, China; R401-01-AA), as previously described with minor modifications.²² Briefly, samples were lysed in extraction reagent, phase-separated with chloroform, precipitated with isopropanol, washed with 75% DEPC-treated ethanol, and dissolved in RNase-free water.

For mRNA, circRNA, and TF detection, 1000 ng of total RNA was reverse-transcribed using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific; K1622). RT-qPCR was performed with ChamQ SYBR qPCR Master Mix (Vazyme; Q311-03) on a QuantStudio 3 Real-Time PCR System (Applied Biosystems). Each 10-μL reaction contained 5 μL 2× master mix, 0.2 μL each primer (10 μM), 1 μL cDNA, and RNase-free water to volume; ROX dye was added when required. Cycling conditions were 95 °C for 30s, followed by 40 cycles of 95 °C for 10s and 60 °C for 30s, with a final melt-curve analysis. Relative expression was calculated with the 2^-ΔΔCt method.²³ GAPDH was used as the internal control for circRNA, mRNA, and TF assays, and U6 was used for miRNA assays. All reactions were run in triplicate. Primer sequences are provided in Table 3. CircRNA primers were designed across the back-splice junction, and specificity was checked by melt-curve analysis.

Table 3 Primer Sequences

Statistical Analysis

Normality of continuous variables was assessed using the Shapiro–Wilk test. All variables analyzed for group comparisons showed P > 0.05 and were therefore considered approximately normally distributed. Animal-experiment data are presented as mean ± SD. Comparisons between two groups at a single time point were performed with an unpaired two-tailed Student’s t test. Paw thickness and arthritis index measured over time were analyzed by two-way repeated-measures ANOVA followed by Sidak’s multiple-comparisons test. Statistical analyses were performed with GraphPad Prism 9.0, and P < 0.05 was considered statistically significant.

Results

Overview of the Analytical Workflow

The differential-expression and enrichment results are summarized in Figure 1. After preprocessing of the three GEO datasets, differential expression analysis identified 980 circRNAs in GSE189338, 282 miRNAs in GSE124373, and 920 mRNAs in GSE17755 that met the predefined screening criteria. These candidate transcripts were then subjected to enrichment analysis, network construction, CTD-based prioritization, and experimental validation in the CIA model.

Figure 1 Differential-expression profiling and functional enrichment in rheumatoid arthritis. (A) Box plots of normalized circRNA, miRNA, and mRNA expression data. (B–D) Heat maps and volcano plots of differentially expressed circRNAs (B), miRNAs (C), and mRNAs (D). (E) GO enrichment analysis of differentially expressed mRNAs. (F) KEGG pathway enrichment analysis of differentially expressed mRNAs.

Functional Characteristics of the mRNA Signature

The differentially expressed mRNAs were enriched in biological processes and pathways closely related to RA immunopathology. Prominent GO terms included leukocyte cell-cell adhesion, regulation of apoptotic signaling, focal adhesion, and ubiquitin-related binding functions. KEGG analysis further supported the inflammatory and immune context of the dataset. These findings are consistent with the established importance of leukocyte recruitment, stromal activation, and survival signaling in RA.^1–3

Construction of the Integrated Regulatory Network

Intersection of the mRNA expression matrix with TRRUST yielded 46 TF-mRNA regulatory pairs involving 6 TFs and 42 target genes. Integration of differentially expressed miRNAs through miRTarBase then produced a miRNA-TF/mRNA network containing 12 miRNAs and their supported targets. CircBank prediction followed by overlap with the differentially expressed circRNA set identified 6 circRNAs that connected to this regulatory core. The resulting initial network comprised 6 circRNAs, 4 miRNAs, 4 TFs, and 24 mRNAs (Figure 2).

Figure 2 Integrated regulatory network of differentially expressed transcripts in rheumatoid arthritis. (A) TF–mRNA regulatory pairs identified from TRRUST. Dashed green lines indicate inhibitory regulation. (B) miRNA–TF/mRNA interactions supported by miRTarBase. (C) Integrated circRNA–miRNA–TF–mRNA ceRNA network.

Prioritization of RT-qPCR Targets

Because the full network contained substantially more transcripts than could be validated in a small animal experiment, we selected a focused set of targets for RT-qPCR. The selected panel consisted of 6 circRNAs (circ0086684, circ0001605, circ0043947, circ0009723, circ0015911, and circ0029987), 4 miRNAs (miR-4646-5p, miR-195-5p, miR-29b-3p, and miR-424-5p), 3 TFs (DNMT1, MYB, and HIF1A), and 8 mRNAs (HDAC1, RB1, GATA3, STAT4, MIF, TLR2, CXCR4, and ENO1). This prioritization was based on network position, cross-layer connectivity, and CTD support rather than on fold change alone.

Validation of the CIA Model

Rats in the CIA group developed progressive paw swelling and increasing arthritis scores after immunization, whereas control animals remained clinically stable Histological examination showed synovial hyperplasia, inflammatory-cell infiltration, and cartilage/bone injury in arthritic joints. Micro-CT demonstrated obvious bone destruction and trabecular damage in the CIA group. Together, these observations confirmed successful establishment of the inflammatory arthritis model (Figure 3).

Figure 3 Experimental validation of the collagen-induced arthritis model. (A) Paw thickness and arthritis index during model establishment. (B) Summary of representative clinical observations in control and CIA rats. (C) Representative hematoxylin and eosin staining of ankle joints in control and CIA rats (50× and 100×). (D) Representative paw photographs and three-dimensional micro-CT reconstructions; red boxes indicate areas of bone erosion.

RT-qPCR Validation and Refinement of the Network

Most selected transcripts showed RT-qPCR trends consistent with the discovery datasets (Figure 4). CircRNAs, the three TFs, and the eight mRNAs were generally increased in the CIA group relative to controls, and three of the four selected miRNAs (miR-4646-5p, miR-29b-3p, and miR-424-5p) were also significantly altered in the expected direction. In contrast, miR-195-5p did not differ significantly between groups and was therefore removed from the validated subnetwork. After exclusion of this unsupported branch, the refined network contained 5 circRNAs, 3 miRNAs, 3 TFs, and 8 mRNAs.

Figure 4 RT-qPCR validation of selected network components in peripheral blood from control and CIA rats. (A) circRNAs. (B) miRNAs. (C) mRNAs. (D) Transcription factors. (*P < 0.05; **P < 0.01; ***P < 0.001 versus control). Abbreviation: ns, not significant.

Interpretation of the Expression-Direction Pattern

An apparent point of confusion is that the overall differential-expression results included both upregulated and downregulated circRNAs and miRNAs, whereas the experimentally retained subnetwork was dominated by transcripts that changed in the same direction. This does not imply that all circRNA-miRNA or miRNA-mRNA relationships are functionally concordant. Rather, it reflects the fact that the validated network represents only a small, selected subset derived from separate discovery datasets and then tested in an animal blood model. Accordingly, the refined network should be regarded as a prioritized regulatory hypothesis rather than proof of a fully resolved ceRNA mechanism.

Discussion

The main aim of this study was to move from a broad computational network to a smaller set of experimentally supported candidates. That is where the main value of this study lies. The work does not establish a causal ceRNA mechanism, but it does narrow the field to a smaller set of blood-associated molecules worth pursuing in follow-up studies.

Several of the retained genes fit well with what is already known about RA biology. TLR2 has a well-established role in innate immune activation in arthritic tissue.²⁴ MIF is a pleiotropic inflammatory cytokine with substantial evidence for pathogenic activity in RA.²⁵ HDAC1 is increased in RA synovium and has been linked to synovial fibroblast survival and inflammatory behavior.²⁶ STAT4 is a recognized RA susceptibility locus and an important transcriptional node in T-cell differentiation.²⁷ HIF-related signaling is also relevant, given the hypoxic synovial microenvironment and its effects on immune-cell function.²⁸ CXCR4 has been associated with leukocyte retention and tissue remodeling in inflamed synovium,²⁹ while alpha-enolase is a well-described autoantigenic and inflammatory molecule in RA.³⁰ Their appearance in the refined network does not prove direct regulation, but it does make the prioritization strategy more credible.

The non-coding RNA layer is less settled. CircRNAs are attractive as blood-accessible molecules because of their relative stability and growing biomarker literature in RA.^5,11–13 By contrast, miR-4646-5p has barely been studied in this disease, so at present it should be treated as a candidate rather than a conclusion. There is more support for miR-424 dysregulation in RA synovial fibroblasts³¹ and for miR-29b in apoptosis-related inflammatory responses in RA monocytes.³² The failure of miR-195-5p to validate in our model is also worth noting: it is a reminder that network predictions look cleaner on paper than they do in biological material.

Several limitations need to be kept in view. The discovery datasets came from independent cohorts and different platforms, and the circRNA dataset was small. The screening thresholds were intentionally lenient because the study was designed for discovery, not for definitive biomarker nomination. The validation step was performed in rat peripheral blood, whereas the discovery datasets were derived from human PBMCs or peripheral blood cells, so species and cellular-composition differences may have contributed to some inconsistency. We also lacked an independent human validation cohort. Finally, we did not perform gain- or loss-of-function experiments, which means the proposed circRNA-miRNA-TF-mRNA relationships remain putative.

For that reason, we view the present work as a prioritization study rather than a mechanistic paper. Its contribution is to produce a more focused shortlist of RA-related regulatory candidates that can now be tested directly in human cells, larger clinical cohorts, and perturbation experiments.

Conclusion

By integrating public transcriptomic datasets with TF information and experimental screening in a CIA model, we identified a refined circRNA-miRNA-TF-mRNA regulatory network associated with RA. The retained molecules represent candidate peripheral blood biomarkers and a tractable framework for subsequent study. These findings provide a focused basis for future human and mechanistic validation.