Unveiling new therapeutic horizons in rheumatoid arthritis: an In-depth exploration of circular RNAs derived from plasma exosomes

Abstract

Rheumatoid arthritis (RA), a chronic inflammatory joint disease causing permanent disability, involves exosomes, nanosized mammalian extracellular particles. Circular RNA (circRNA) serves as a biomarker in RA blood samples. This research screened differentially expressed circRNAs in RA patient plasma exosomes for novel diagnostic biomarkers. In this study, samples of RA patients with insufficient response to methotrexate (MTX-IR), combined use of tumor necrosis factor inhibitors (TNFi) were followed up for half a year, and 56 circRNA samples of self-test data were stratified into training, testing, and external validation cohorts according to whether American College of Rheumatology 20% improvement criteria (ACR20) was achieved. A diagnostic xgboost model was developed using common hub genes identified by random forest and least absolute shrinkage and selection operator (LASSO), with intersection genes derived from overlapping machine learning-selected genes. Diagnostic performance evaluated via receiver operating characteristic (ROC) curves using pROC for area under the curve (AUC). Optimal LASSO model with 4 circRNAs determined, with AUC > 0.6 for key genes. The model validation performed well on the test set, but not significantly on the validation set. Then, circRNA screening was performed in combination with clinical data, and cross-validation identified hsa-circ0002715, hsa-circ0001946, hsa-circ0000836, and rheumatoid factor (RF) as key genes, among which hsa-circ0002715 and hsa-circ0001946 were emphasized as key markers on the training set. In addition, the morphology and size of exosomes and the expression of CD9 and CD81 verified the successful extraction of exosomes. The qPCR analysis of plasma exosomes in TNFi patients found that the expression of hsa-circ0002715 was higher than that in patients who didn’t reach ACR20, and the expression of hsa-circ0001946 was lower than that in patients who didn’t reach ACR20. The above studies suggested that hsa-circ0002715 and hsa-circ0001946 may become markers for predicting MTX-IR RA patients and TNFi precision treatment.

Introduction

Rheumatoid arthritis (RA) a common autoimmune disease with systemic impact, is characterized by synovial hyperplasia, persistent joint inflammation, and diverse extra-articular presentations [1, 2]. Females exhibit a higher predisposition to RA, with a prevalence rate of 3.6%, in contrast to males, who have an estimated risk of 1.7% [3]. RA not only results in joint deformities and a potential for lifelong disability but also markedly elevates the risk of comorbidities, notably cardiovascular disease and malignancies [4, 5]. Tumor necrosis factor inhibitors (TNFi), a representative biological disease-modifying antirheumatic drugs (bDMARDs), is the primary second-line treatment for RA [6, 7]. Nevertheless, a considerable proportion of patients experience treatment failure due to lack of efficacy [8].

Treatment failure in RA escalates the risk of drug-related complications and imposes a considerable burden on both patients and society [9]. Hence, predicting the efficacy of individualized TNFi regimens for RA patients is a research focus for RA treatment. Recently, research on the prediction of efficacy in RA patients is mainly based on clinical data or conventional biomarkers [10]. Multiple studies have leveraged machine learning approaches to forecast treatment efficacy by integrating clinical, genetic marker, and multi-omics data [11,12,13]. Previous biomarker studies focused on RA’s efficacy mechanisms at gene expression and modification levels [14, 15]. Further investigation into biomarkers at cellular and gene-cell interfaces is crucial for accurate RA diagnosis and prediction.

Exosomes are membrane-bound vesicles released into the extracellular space subsequent to the fusion of intracellular multivesicular bodies with the plasma membrane [16, 17]. Extracellular vesicles also have certain potential uses in several disease [18, 19]. In addition, exosomal circRNA represents an endogenous, highly stable non-coding RNA that exhibits aberrant expression patterns in disease states and is amenable to facile collection and detection [20, 21]. Increasing research reveals the role of noncoding RNA in musculoskeletal diseases [22,23,24,25]. Extensive research has affirmed the pivotal role of circRNAs in modulating immune and inflammatory responses in RA, underscoring their significance as a crucial biomarker for evaluating RA disease status [26, 27]. However, the predictive role of exosomal circRNAs in RA prognosis and its potential as personalized indicators of the efficacy of second-line TNFi treatment in RA patients remains unknown. This study enrolled RA patients who failed methotrexate (MTX) initial treatment scheduled for TNFi therapy, assessing exosomal circRNA levels in their peripheral blood. Patients were followed post-TNFi treatment and categorized based on their treatment outcomes. Subsequently, circRNA was screened for machine learning features, and random forest was used to construct a treatment efficacy prediction model based on circRNAs indicators and patient clinical characteristics. Finally, plasma exosomes were constructed and isolated, and circRNA detection was performed for verification. This research endeavors to predict TNFi treatment effects in RA patients, offering the clinical treatment for RA patients.

Materials and methods

Study subjects

Between October and December 2023, patients with RA who had poor response to initial MTX treatment and continued to use TNFi in combination at the Affiliated Hospital of Yangzhou University, followed up for half a year, according to the 2010 ACR/EULAR classification diagnostic criteria for RA. At baseline, peripheral blood samples were procured from these patientsand preserved in sodium citrate-containing vacuum tubes. Exclusion criteria for healthy controls encompassed severe cardiovascular disease, hepatic and renal impairment, inflammatory disorders, malignancies, and other immune-mediated conditions (e.g., systemic lupus erythematosus, ankylosing spondylitis). This study was conducted by the Declaration of Helsinki and approved by the Ethics Committee of the Affiliated Hospital of Yangzhou University [approval number: 2023-YKL10-(05)], and all subjects signed written informed consent. We collected clinical data from patients, including demographics (age, gender, etc.), smoking history, serological indicators (including RF, anti-CCP, etc.), baseline DAS28-ESR, DAS28-CRP, CDAI, SDAI, HAQ, SHARP scores, etc.

Data Processing

The self-test dataset consisted of 56 circRNA samples, which were partitioned into training set, test set and external validation set by 0.65, 0.15, and 0.2. Among them, feature screening only uses training set and test set data, and the validation set is only verified in the final model.

Screening and validation of diagnostic markers

This study employed diverse machine learning algorithms, including random forest, least absolute shrinkage and selection operator (LASSO) logistic regression, and XGBoost, to identify novel key biomarkers for RA [28, 29]. In this study, the “randomForest” package in R was utilized for random forest analysis, while the “glmnet” package facilitated LASSO analysis with 10-fold cross-validation incorporating circRNA and clinical features. Additionally, the “caret” package was employed to perform random forest and XGBoost analyses, yielding the top 20 core genes from each algorithm. Subsequently, an intersection of genes identified by these three machine learning approaches was computed to obtain the consensus gene set.

The diagnostic efficacy of the identified biomarkers was assessed through the construction of receiver operating characteristic (ROC) curves for the self-test dataset. The area under the curve (AUC) was calculated utilizing the “pROC” package to evaluate the diagnostic performance. Biomarkers with AUC values exceeding 0.6 were selected as the final screened biomarkers.

Construction of random forest diagnostic model

Based on the identified biomarkers, a joint model was constructed and validated against both the test set and an external validation set to ensure its reliability. ROC curves were plotted to evaluate the sensitivity and specificity of the models, with AUC values calculated to quantify their overall reliability.

Extraction of plasma exosomes and exosomal RNA

Blood samples from 56 RA patients underwent centrifugation at 500 g for 10 min to isolate plasma, which was subsequently stored at -80 °C. Prior to exosome isolation, plasma samples were centrifuged to eliminate cellular debris. Exosomes were then harvested via ultracentrifugation at 150,000 g for 2 h, and their concentration was determined using a BCA assay. To facilitate downstream analysis, exosomes from each group were pooled. Adhering to the manufacturer’s instructions, exosomal RNA was extracted using a rapid total RNA extraction kit.

Transmission electron microscopy (TEM) and nanoparticle tracking analysis (NTA) particle size analysis

The purified exosomes were deposited onto a carbon-coated copper grid, fixed with 2.5% glutaraldehyde, and stained with uranyl acetate prior to examination. Subsequently, 2% phosphotungstic acid was applied to the grid for 2 min, followed by transmission electron microscopy (TEM) imaging for analysis.

The Zetaview (Particle Metrix) system was employed for nanoparticle tracking analysis. Upon addition of exosomes to the instrument’s loading port, it automatically determined the hydrodynamic diameter and concentration of the nanoparticles through a combination of single-particle tracking technology, classical microelectrophoresis (zeta potential), and the principles of Brownian motion, utilizing the Stockes-Einstein equation.

Western blotting

Proteins isolated from exosomes underwent separation via 12% SDS-PAGE and were subsequently transferred onto a PVDF membrane (Merck, USA). The membrane was blocked with 5% BSA in TBST and then probed with antibodies against CD9 (ab236630, 1:1000, Abcam) and CD81 (ab79559, 1:1000, Abcam). Following incubation with HRP-conjugated secondary antibodies, chemiluminescent signals were detected using the ECL method and visualized on a Tanon 5200 integrated chemiluminescent imager (China).

The validation of diferentially expressed circRNAs

Total RNA extracted from exosomes was reverse transcribed into cDNA using the RevertAid First Strand cDNA Synthesis Kit (Thermo, USA). For qPCR, a 20 µl reaction mixture was prepared containing 10 µl of 2x ChamQ Universal SYBR qPCR MasterMix (Roche, Switzerland), cDNA, and 0.5 µl of each forward and reverse primer. Gene expression was quantified relative to the housekeeping gene GAPDH.

Statistical analysis

Statistical analyses were conducted using R version 4.3.3 and GraphPad Prism 9. Comparisons between the two sample groups were performed using either the Wilcoxon rank sum test or t-test, as appropriate. An ROC curve was generated, and the sensitivity, specificity, and AUC, along with their 95% confidence intervals, were evaluated. The p-value less than 0.05 was considered statistically significant, ns represents p > 0.05, ^* represents p < 0.05, ^** represents p < 0.01, and ^*** represents p < 0.001.

Results

Screening of biomarkers in the circRNA signature construction diagnostic model

The clinical data of this group showed that in patients who achieved ACR20 after TNFi treatment, the levels of biochemical indicators Alanine Transaminase (ALT), Aspartate Aminotransferase (AST) and lactate dehydrogenase (LDH), immunological indicators Erythrocyte sedimentation rate (ESR) and C-reactive protein (CRP), and rheumatoid arthritis disease assessment indicators Disease Activity Score for 28 joints based on the Erythrocyte Sedimentation Rate (DAS28-ESR), and C-reactive protein level (DAS28-CRP), the simplified and clinical disease activity indices (SDAI, CDAI) decreased 6 months after treatment; while in patients who did not achieve ACR20, the expressions of neutrophil (NE), platelet (PLT), ALT, D-Dimer, anticitrullinated peptides/protein antibodies (ACPA), mutated and citrullinated vimentin (MCV), RF-IgM, and RF-IgA increased significantly, and LDH, CRP, SDAI, CDAI, joint tenderness, and joint swelling decreased 6 months after treatment (Table 1; Fig. 1A-D). Then, key genes were screened by LASSO regression analysis, random forest and xgboost. When 4 circRNAs were included in the LASSO model, its performance was the best. In random forest and XGBOOST, the weights of the top 10 genes were hsa-circ0002715, hsa-circ0001946, hsa-circ0000835, etc. (Fig. 2A-C). The intersection of the three key features was taken, and it was found that the three RNAs hsa-circ0002715, hsa-circ0001946, and hsa-circ0000835 were all important, so they were selected as key genes (Fig. 2D).

Table 1 Clinical features of the participants

Fig. 1

The clinical data of the groups. (A) Immunological indicators Erythrocyte sedimentation rate (ESR), C-reactive protein (CRP), rheumatoid arthritis disease assessment indicators Disease Activity Score for 28 joints based on the Erythrocyte Sedimentation Rate (DAS28-ESR), and C-reactive protein level (DAS28-CRP); (B) Clinical disease activity index (CDAI) and simplified clinical disease activity index (SDAI); (C) Bone erosion and joint space of both hands were evaluated by Sharp-van der Heijde method (Sharp); (D) Efficacy evaluation parameters tenderness number of joint; swelling of joint; Pain assessment by patients (1–10 cm); Overall evaluation of disease activity by patients (PGA); Overall evaluation of disease activity by doctors (MDGA); Evaluation of physical function by patients (HAQ)

Fig. 2

Screening of biomarkers in the circRNA signature construction diagnostic model. (A) Mean square error in Lasso regression versus Log (λ); (B) Regression coefficient versus Log (λ); (C) Random forest and xgboost top20 score graph; (D) Venn diagram for machine learning feature selection. The numbers represent the number of genes in each interval, and the percentages come from the proportion of the number of genes to all the genes screened by the three methods

Validation of screening results of biomarkers in the circRNA signature construction diagnostic model

We then analyzed whether there were differences in the biomarkers in different outcomes and their evaluation capabilities. The analysis results showed that only two markers had significant differences, but each marker had a good ability to predict the outcome, and the AUCs of the three-feature box plot and the single-feature ROC curve are 0.757, 0.875, and 0.667, respectively, all greater than 0.6 (Fig. 3A, B).

Fig. 3

Validation of screening results of biomarkers in the circRNA signature construction diagnostic model. (A) Feature box plot; (B) Single feature ROC curve

Diagnostic modeling of circRNA signature biomarkers

Then, the random forest method was used to model the biomarkers and evaluate the constructed validation model. Both the test set and the validation set were used for evaluation. The left side is the test set and the right side is the validation set. It can be clearly seen from the molecular results that the test set has a better effect, but the effect of the validation set is not very obvious (Fig. 4A, B). Therefore, we further combined clinical characteristics for modeling.

Fig. 4

Diagnostic modeling of circRNA signature biomarkers. (A, B) ROC curve of test set and validation set

Biomarker screening based on circRNA combined with clinical characteristics

In order to screen biomarkers of circRNA combined with clinical characteristics, key genes were screened by LASSO regression analysis, random forest and xgboost. When the LASSO model contains 4 features, its performance is the best (Fig. 5A, B). The weights of the top 10 genes in random forest and XGBOOST are as follows: hsa-circ0002715, hsa-circ0001946, hsa-circ0000836, etc. The intersection of the three key features was taken, and it was found that the four features of hsa-circ0002715, hsa-circ0001946, hsa-circ0000836 and rheumatoid factors (RF) were all important, so they were selected as key genes (Fig. 5C, D).

Fig. 5

Biomarker screening based on circRNA combined with clinical characteristics. (A) Mean square error in Lasso regression versus Log (λ); (B) Regression coefficient versus Log (λ); (C) Random forest and xgboost top20 score graph; (D) Venn diagram for machine learning feature selection

Validation of biomarker screening results in diagnostic models constructed by combining circRNA with clinical characteristics

Then, we tested whether there were differences in biomarkers in different outcomes and their ability to evaluate outcomes. The four-feature box plot and the single-feature ROC curve showed that only two markers, hsa-circ0002715 and hsa-circ0001946, had significant differences, and the other features, hsa-circ0000836 and RF, had certain trends (Fig. 6A-E).

Fig. 6

Validation of biomarker screening results in diagnostic models constructed by combining circRNA with clinical characteristics. (A) Feature box plot; (B-E) single feature ROC curve

CircRNA combined with clinical characteristics of biomarkers for diagnostic modeling

The random forest method was used to model the biomarkers and evaluate the constructed validation model. Both the test set and the validation set were used for evaluation, with the test set on the left and the validation set on the right. It can be clearly seen from the ROC curves of the test set and the validation set that the performance of the model has improved, and the AUC of the training set has reached 0.778 (Fig. 7A, B). In addition, although RF cannot diagnose the outcome, it can significantly improve the predictive ability of the circRNA model.

Fig. 7

CircRNA combined with clinical characteristics of biomarkers for diagnostic modeling. (A, B) ROC curve of test set and validation set

Detection of exosomal circRNA levels in patients’ plasma

Next, we examined the characteristics of exosomes isolated from plasma. Transmission electron microscopy (TEM) and nanoparticle tracking analysis (NTA) analysis showed that the exosomes had visible double-layer capsules or crescent-shaped ultrastructures, both of which were exosome morphology (Fig. 8A). The results of NTA particle size analysis showed that the exosome concentration was 4.3E + 6 particles/ml, the exosome particle size was about 105 nm, and the size range was 30 to 200 nm, which was consistent with the previous characterization of exosomes (Fig. 8B). The presence of exosome surface markers (including CD9 and CD81) was confirmed by Western blotting, further proving that the isolated particles were exosomes (Fig. 8C). Then qPCR was used to detect the expression of exosomal RNA, and our data revealed that in the plasma exosomes of patients with TNFi reaching ACR20, the expression of hsa-circ0002715 was higher than that of patients who did not reach ACR20, and the expression of hsa-circ0001946 was lower than that of patients who did not reach ACR20, and the difference was significant. There was no significant difference in the other circRNAs (Fig. 9). The analysis flow chart of this study is shown in Fig. 10.

Fig. 8

The characteristics of exosomes isolated from plasma. (A) Transmission electron microscopy photo of exosomes; (B) NTA particle size analysis; (C) Western blot for marker detection

Fig. 9

The qPCR was used to detect the expression of exosomal RNA

Fig. 10

Analysis flowchart for this study

Discussion

RA is a multifaceted autoimmune disorder marked by the aggressive proliferation and infiltration of fibroblast-like synoviocytes, ultimately resulting in joint destruction [30, 31]. While the precise pathological mechanisms underlying RA remain elusive, exosomes are suspected to contribute to disease progression through their ability to deliver diverse biomolecules to RA-affected tissues [32]. Among the molecules present in exosomes, circRNA stands out as a potentially disease-exacerbating factor due to its high stability and previously established role in the pathogenesis of RA [33,34,35]. Furthermore, circRNAs exhibit robust expression across various tissues, with notable enrichment in the human brain and blood [36]. However, limited research exists on leveraging circRNAs in RA patient’s peripheral blood exosomes to forecast the effectiveness of TNFi second-line treatment strategies.

Recent studies have demonstrated that an array of circRNA species are implicated in the progression of RA [37]. Moreover, circRNAs have emerged as biomarkers with enhanced specificity for particular diseases and have been found to maintain superior stability across diverse organisms [38]. Therefore, targeting circRNAs might serve as a promising approach for the treatment of RA. In this study, we recruited RA patients who had failed initial MTX treatment and received TNFi treatment, and obtained a total of 56 circRNA samples. The data was divided into training set, test set, and external validation set, and then the machine learning features of circRNA were screened, and a random forest was used to build an efficacy prediction model based on circRNA combined with clinical features. Finally, plasma exosomes were constructed and isolated, and circRNA detection was performed for verification.

Previously, a research has explored the regulatory network of programmed cell death genes in rheumatoid arthritis based on blood-derived circRNA transcriptome information through LASSO analysis [39]. Here, we constructed a diagnostic model based on circRNA features, screened key genes through LASSO regression analysis, random forest and xgboost, and found that the LASSO model performed best when it included four circRNAs. In random forest and XGBOOST, the weights of the top 10 genes were hsa-circ0002715, hsa-circ0001946, hsa-circ0000835, etc. Then we took the intersection of the three key features and found that hsa-circ0002715, hsa-circ0001946, and hsa-circ0000835 were all important, so they were selected as key genes. Our earlier research endeavor uncovered that the irregular expression of circRNA_0001946 stimulates colorectal cancer cell proliferation and metastatic potential by modulating microRNA-135a-5p [40]. We subsequently analyzed biomarker differences across various outcome indicators and their predictive capabilities, finding that two markers were significantly distinct, each with good predictive power. Both the AUC of the three-feature box plot and the single-feature ROC curve exceeded 0.6. Modeling these biomarkers using the random forest method showed good performance in the test set but less so in the validation set. Therefore, we incorporated clinical characteristics into our modeling approach.

In subsequent analyses, we identified circRNA biomarkers associated with clinical manifestations, optimizing performance with a four-feature combination using the LASSO model. Among the top 10 genes prioritized by random forest and XGBOOST, hsa-circ0002715, hsa-circ0001946, and hsa-circ0000836 emerged as significant. RF analysis further confirmed their importance, selecting them as pivotal genes. Evaluation of biomarker expression across outcome indicators revealed statistically significant differences for hsa-circ0002715 and hsa-circ0001946 (via feature box plots and ROC curves), with notable trends for hsa-circ0000836 and RF. Luo et al. also identified circular RNA Hsa_circ_0002715 in peripheral blood as a novel potential biomarker for new-onset RA [41]. In addition, the ROC curves of the test set and validation set clearly show that the model performance has improved, and the AUC of the training set has reached 0.778. Although RF cannot diagnose the results, it can significantly improve the predictive ability of the circRNA model.

Subsequently, TEM and NTA confirmed the presence of exosomes with discernible double-layered or crescent-shaped ultrastructures, both characteristic of exosome morphology, aligning with previous research findings [42]. The NTA particle size analysis yielded an exosome concentration of 4.3E + 6/ml, with an average particle diameter of approximately 105 nm and a size distribution spanning from 30 to 200 nm, consistent with previously reported exosome characterization outcomes [43]. Western blotting analysis validated the expression of exosomal surface markers CD9 and CD81, reinforcing the identity of the isolated particles as exosomes [44]. Furthermore, qPCR results showed that in the plasma exosomes of patients with TNFi who achieved ACR20, the expression level of hsa-circ0002715 was higher than that of patients who did not achieve ACR20, and the expression level of hsa-circ0001946 was lower than that of patients who did not achieve ACR20, and the difference was statistically significant.

However, this study has certain limitations. First, the sample size employed was comparatively modest, necessitating the validation of these findings in larger cohorts and diverse demographic groups, encompassing various ethnicities and geographical regions. Secondly, the investigation did not assess the involvement of hsa-circ0002715 and hsa-circ0001946 in the etiopathogenesis of RA. Consequently, additional experimental endeavors are warranted to examine the potential causal link between the aberrant expression of these circRNAs and the onset or progression of RA. A deeper understanding of the role of circRNAs in RA pathogenesis is prerequisite for the routine integration of novel diagnostic and therapeutic modalities.