A Longitudinal Study Reveals Metabolomic Markers for Individuals at Risk, Disease Severity, and Treatment Response in Rheumatoid Arthritis

Abstract

Rheumatoid arthritis (RA) is a systemic inflammatory joint disease characterized by heterogeneous clinical manifestations, which requires deeper exploration in identifying reliable biomarkers for early diagnosis, monitoring, and treatment assessment. The aim is to discover plasma metabolomic markers to predict RA onset, assess disease activity, and forecast treatment efficacy. The study includes 209 established RA patients who are disease-modifying antirheumatic drugs-free for six months prior to enrollment, with 197 of them followed for 3-6 months to assess treatment response. Additionally, 56 individuals at risk are recruited, with 34 completing a 5-7-year follow-up. Analysis reveals that metabolites related to methylation and redox imbalance, such as S-adenosylmethionine, sarcosine, nicotinamide adenine dinucleotide, glutathione, etc., are associated with RA development and severity, and contribute to its heterogeneity across age, sex, and anti-citrullinated protein autoantibody status. Ridge regression models are constructed using metabolite and clinical features for the response to methotrexate (MTX) plus leflunomide, achieving an average receiver operating characteristic (ROC) score of 0.83, and for the MTX plus hydroxychloroquine, achieving an average ROC score of 0.92. In conclusion, our findings reveal RA metabolomic alterations, aiding early diagnosis and treatment response.

Keywords: longitudinal profiling; machine learning; metabolomics; rheumatoid arthritis; therapy response.

Figure 1

Metabolomic analysis workflow and quality control analysis. A) Schematic illustration of the plasma metabolomics workflow, including sample collection, targeted metabolomic detection, and data analyses. B) Spearman's correlation coefficients of quality controls (QCs) samples in both positive (red) and negative (blue) ion modes. The upper panels show color‐coded pairwise correlation values between QCs; the lower panels display pairwise scatter plots. Density curves of each QCs are shown along the diagonal. C) The cumulative count of metabolites identified across health (blue, n = 99), individuals at risk (IAR) (yellow, n = 56), and rheumatoid arthritis (RA) groups (red, n = 209). D) Coefficient of variation (CV) values for 185 metabolites across QCs samples. E) Classification and distribution of the 185 metabolites by superclass. Classification is annotated using the Human Metabolome Database.

Figure 2

Plasma metabolomic differences among RA, IAR, and healthy individuals. A) Boxplots illustrating the median and interquartile range of normalized metabolites’ abundance (n = 185). B) The partial least squares discriminant analysis (PLS‐DA) scatter plot displaying the distribution of samples across baseline RA (red), IAR (yellow), and health (blue) on the first two principal components. C, D, E) Volcano plot showing metabolite associations across clinical groups (RA vs. health, RA vs. IAR, and IAR vs. health) via generalized estimating equation (GEE) model, adjusted for age, sex, body mass index (BMI), smoking status, and comorbidities. Metabolites with significant positive associations (coefficient > 0, p < 0.05, two‐sided t‐test) are shown in red; significant negative associations (coefficient < 0, p < 0.05, two‐sided t‐test) are shown in blue. The top 10 most significant metabolites are highlighted in gold. The gray shaded region represents metabolites with p‐values smaller than the minimum detectable threshold in R (2.2E‐16). F) Forest plot of regression coefficients (95% CI) for metabolites significantly altered in RA versus health and IAR versus health comparisons, involved in inflammation, methylation, or redox processes. G) Heatmap showing the median gene set variation analysis (GSVA) scores of differential pathways across three clinical groups. Significance was determined using GEE analysis (p < 0.05, two‐sided t‐test), adjusting for age, sex, BMI, smoking status, and comorbidities.

Figure 3

Influences of sex, age, and ACPA status on metabolomic profiles in RA. A) Heatmap presenting baseline clinical data for 209 RA patients, including disease activity groups (clinical remission: DAS28‐CRP < 2.6; low disease activity: 2.6 ≤ DAS28‐CRP ≤ 3.2; moderate disease activity: 3.2 < DAS28‐CRP ≤ 5.1; high disease activity: DAS28‐CRP > 5.1), the percentage composition of DAS28‐CRP parameters: shared_TJC28 = $\frac{0.56\times \sqrt{\mathrm{T}JC28}}{\mathrm{D}AS28\_CRP-0.96}$, shared_SJC28 = $\frac{0.28\times \sqrt{\mathrm{S}JC28}}{\mathrm{D}AS28\_CRP-0.96}$, shared_CRP = $\frac{0.36\times \ln (\mathrm{CRP}+1)}{\mathrm{D}AS28\_CRP-0.96}$, shared_VAS =$\frac{0.014\times (\mathrm{VAS})}{\mathrm{D}AS28\_CRP-0.96}$, anti‐citrullinated protein autoantibodies (ACPA), age, sex, BMI, smoking and commodities status, with blanks indicating missing data. B–D) Volcano plots illustrating sex‐associated metabolite differences within RA (B), IAR (C), and health (D). Red dots indicate metabolites significantly higher in males, whereas blue dots represent metabolites significantly lower in males (|fold change|>1.2, p < 0.05, Wilcoxon test). E) Venn diagram showing overlapping significant metabolites with increased levels in males compared to females across RA, IAR, and health. F) Volcano plots showing metabolites with linear age associations based on GAMs model in RA (left), IAR (middle), and health (right). Red dots indicate metabolites significantly increased with age (coefficient > 0, p < 0.05, two‐sided t‐test). G) Venn diagram showing overlapping metabolites significantly increased with age between RA and IAR. H) Line plots of age‐associated metabolites with non‐linear trends, clustered by their fitted trajectories from GAMs. I) Density plot from 100 iterations of GEE model regression (adjusted for age, sex, and CRP) comparing 14 randomly selected ACPA positive RA patients with 14 ACPA negative RA patients, significance of regression coefficients was determined by two‐sided t‐test. J) Violin plots with embedded boxplots (the median and interquartile range) displaying the results of the Wilcoxon test for differential metabolites between 183 ACPA‐positive RA and 14 ACPA‐negative RA, with ^* p < 0.05 and ^** p < 0.01, “ns” indicates not significant.

Figure 4

Investigation into the metabolites linked to development and disease activity in RA. A) Left: Analytical framework for assessing metabolite differences between converters (developed into RA) and non‐converters (not developed into RA) in IAR. Middle: Bar plot showing the significant counts of metabolites (p < 0.05, Wilcoxon test) between 4 converters and 4 randomly selected non‐converters across 100 iterations, with 7 top metabolites exceeding 10 significant counts highlighted. Right: Boxplots (the median and interquartile range) showing the distribution of log10(fold change) across significant iterations for the 7 top metabolites, with the red dashed line indicating log10(0). B) Left: Violin plots with embedded boxplots (median and interquartile range) depicting cluster‐specific metabolite level variations, with significance (p < 0.05, Wilcoxon test) determined via pairwise disease activity group comparisons. Middle: Heatmap displaying the median levels (scaled) of metabolites in distinct clusters. Right: Representative metabolites for each cluster. C) GEE model analysis of DAS28‐CRP and metabolites (adjusted for age, sex, BMI, smoking status, and comorbidities), with red indicating positively related (coefficient > 0, p < 0.05, two‐sided t‐test) and blue indicating negatively related (coefficient < 0, p < 0.05, two‐sided t‐test). D) Venn diagram showing the overlap of significant metabolites identified through the Wilcoxon test and GEE model analysis. E) The dot plot displays the correlation of significant metabolites from comparisons between disease activity groups with four DAS28‐CRP parameters based on Spearman analysis. Gray dashed lines denote the statistical significance threshold (p = 0.05).

Figure 5

Predicting drug response and investigating drug effects. A) Volcano plots showing the differential metabolites between responders (n = 30) and non‐responders (n = 43) in the MTX+LEF group identified by the Wilcoxon test (p < 0.05). B) Volcano plots illustrating the differential metabolites between responders (n = 20) and non‐responders (n = 27) in the MTX+HCQ group by Wilcoxon test (p < 0.05). C) Flowchart of constructing a machine learning model to predict drug response. D) Performance evaluation of the machine learning model for predicting response to MTX+LEF based on clinical and metabolic features across converged iterations. (Left) Distribution of area under the receiver operating characteristic curve (AUROC) values for train (red) and test (blue) sets. (Middle) ROC curves for the train set. (Right) ROC curves for the test set. (Far Right) Coefficient contributions of features to the predictive model. E) Performance evaluation of the machine learning model for predicting response to MTX+HCQ based on clinical and metabolic features across converged iterations. (Left) Distribution of AUROC values for train (red) and test (blue) sets. (Middle) ROC curves for the train set. (Right) ROC curves for the test set. (Far Right) Coefficient contributions of features to the predictive model. F) Volcano plots showing differential metabolites between baseline and follow‐up samples of responders in the MTX+LEF (n = 30) by paired Wilcoxon test (p < 0.05). G) Volcano plots showing differential metabolites between baseline and follow‐up samples of responders in the MTX+HCQ (n = 20) by paired Wilcoxon test (p < 0.05). H) Violin plots with embedded boxplots (the median and interquartile range) depicting significantly enriched pathways among health (blue, n = 99), baseline (red, n = 30), and follow‐up samples (yellow, n = 30) of responders in MTX+LEF groups. Significance was calculated by the Wilcoxon test using GSVA scores. Significance levels: ^* p < 0.05, ^** p < 0.01, ^*** p < 0.001, “ns” indicates not significant.