Progression to rheumatoid arthritis in at-risk individuals is defined by systemic inflammation and by T and B cell dysregulation

Abstract

Rheumatoid arthritis (RA) is preceded by an at-risk stage of disease that can be marked by the presence of anticitrullinated protein antibodies (ACPAs) but the absence of clinically apparent synovitis (clinical RA). Preemptive intervention in at-risk individuals could prevent or delay future tissue damage; however, the immunobiology of this stage is unclear. Using integrative multiomics, we longitudinally profiled at-risk individuals, where one-third of participants developed clinical RA on study. We found evidence of systemic inflammation and signatures of activation in naïve T and B cells of at-risk individuals. During progression to clinical RA, proinflammatory skewing of atypical B cells and expansion of memory CD4 T cells with signatures of activation and B cell help were present without elevations in circulating ACPA titers. Epigenetic changes in naïve CD4 T cells suggested a predisposition to differentiate into effector cells capable of B cell help. These findings characterize pathogenesis of the ACPA⁺ at-risk stage and support the concept that the disease begins much earlier than clinical RA. Additionally, an extensive immune resource of the at-risk stage and progression to clinical RA with interactive tools was developed to enable further investigation.

Fig. 1:. Active inflammation is observed in ARI prior to disease onset.

(A) Overview of study and multimodal workflow. HC1, controls; ARI, at-risk individuals; ERA, early RA; CONV, converters. (B) First sample (baseline) ACPA (anti-CCP3) measurements from HC1, NONC (non-converters), CONV, and ERA. (C) k-means clustering (k=6) of z-scored normalized protein expression (NPX) values from differential proteins in ARI vs. HC1 and ERA vs. HC1 (FDR<0.1). Rows denote proteins, columns denote baseline samples. ARI are annotated as CONV or NONC. (D) Abundance of select inflammatory plasma proteins elevated in ARI. Dots represent participant samples in each cluster. (E) Absolute concentration of select plasma proteins from participants in Prot-C1 and Prot-C6 clusters as assayed by Meso Scale Discovery or LegendPlex. (F) Number of differentially expressed genes (DEGs; FDR<0.1 and absolute log2 fold change≥0.1) per immune cell type, elevated in ARI (above 0) or HC1 (below 0). Cell types are based on (24). CM, central memory; EM, effector memory; ISG, interferon-stimulated gene. Boxplots show median (centerline) and first and third quartiles (lower and upper bound of the box); whiskers show the 1.5x interquartile range of data. Effect sizes and P values were determined by linear regression models (C), Wald test (F), Kruskal-Wallis test with Dunn’s post-hoc testing (B), or Wilcoxon rank-sum test (D and E). FDR values are indicated for all panels and values below 0.1 were considered significant.

Fig. 2:. Longitudinal changes in naïve and CM CD4 T cells dominate progression to clinical RA.

(A) Overview of longitudinal comparison of converters (CONV) from ‘at-risk’ to clinical RA. (B) Number of genes per cell type with higher average intra-donor coefficients of variation (CVs) over time in CONV during progression to clinical RA (orange) or in HC2 (green). cDC, conventional dendritic cell (DC); NK, natural killer cell; pDC, plasmacytoid DC. (C) Comparison of the number of differentially expressed genes (DEGs) (y-axis) with the change in frequency over time (x-axis; centered log-ratio (CLR) transformed) as CONV progress to clinical RA. Bubble size corresponds to the number of DEGs. All cell type abundance changes were above FDR of 0.1. (D) Number of DEGs from longitudinal model (FDR<0.1) per cell type, elevated (above 0) or diminished (below 0) in CONV progressing to clinical RA. CM, central memory; EM, effector memory; ISG, interferon-stimulated gene. (E) Overview of paired comparison in converters at their last ‘at-risk’ pre-symptomatic visit vs. time of their clinical RA diagnosis. (F) Normalized RNA expression of TNF in Core CD16 monocytes and CXCL10 in ISG+ CD16 monocytes. (G) Mean RNA expression of select inflammatory genes amongst monocyte subtypes. IL1B+ CD14 monocytes (red) have highest expression of pro-inflammatory genes. (H) Gene scores calculated by comparing marker genes from FOLR2+ICAM+ RA synovial tissue macrophages (30) among all monocyte cell types. STM, synovial tissue macrophages. (I) Frequency of IL1B+ CD14 monocytes within total CD14 monocytes. Violin plot shape is proportional to the density of data points (dots). Larger width represents a higher density of data points. Effect sizes and P values were determined by linear mixed effect models (C, D), paired Wald test (F), ANOVA (H) or paired Wilcoxon test (I). Nominal P values are indicated for (H, I). FDR values are indicated for remaining panels. Values below 0.1 were considered significant.

Fig. 3:. The B cell compartment exhibits a pro-inflammatory skewing during progression to clinical RA.

(A) Uniform Manifold Approximation and Projection (UMAP) plots of MBCs from ARI, HC1 and ERA showing B cell population labels (left) and Leiden clusters (right). (B) DEGs for Beff-C8 compared with Beff-C9 with selected genes labeled (left). Dot size in heatmap (right) indicates the fraction of cells with positive expression for selected effector population marker genes. (C) Specified IgH isotype identity, as frequency within each population, for Beff-C8 and Beff-C9. UND, undetermined. (D) IGHG3 gene expression by core naïve B cells of ARI and HC1. (E) Normalized expression of IgH constant gene germline transcription of IGHMD⁺ naïve B cells that differ between ARI and HC1. (F) CLR-transformed cytometry frequencies of naïve B cell cluster Bnve-S5 as CONV progress to clinical RA. Each participant’s longitudinal series is connected by a gray line, with a group trendline and 95% confidence interval in purple. (G) GSEA enrichment analysis with the top Reactome pathways among naïve B cells of ARI compared with HC1. (H and I) B cells were stimulated ex vivo and analyzed by intracellular flow cytometry. Experimental workflow (H) and RANKL+, IL-6+ and TNF+ cell frequencies within the stimulated naïve B cell populations of ARI and HC2 (I). Boxplots show median (centerline) and first and third quartiles (lower and upper bound of the box); whiskers show the 1.5x interquartile range of data. Violin plot shape is proportional to the density of data points (dots). Larger width represents a higher density of data points. P values were determined by a linear mixed model (B, G), Wald test (E, F), or Wilcoxon rank-sum test (J). FDR values are indicated for all panels. Values below 0.1 were considered significant.

Fig. 4:. Effector and memory T cells with pathogenic signatures expand during progression to clinical RA.

(A to C) Longitudinal analysis of CONV who progress to clinical RA. (A) RNA expression differences in central memory (CM) CD4 T cells. Genes associated with T cell activation are noted. (B) T cell RNA activation metric in CM CD4 T cells. See Supplemental Methods for details. Each participant’s longitudinal series is connected by a gray line, with a group trendline and 95% confidence interval in purple. (C) Frequency of CD4mem-C3 cells over time, as in (B). (D) Cells expressing Tfh gene program are distinguished based on the non-negative matrix factorization projection using a pre-computed weight matrix of CD4 T cells from (33) (left). A UMAP density plot of cluster CD4mem-C3 is shown (right). (E) Gene scores calculated by comparing the Tph/Tfh gene signature from Zhang et al. (28) among all CD4mem clusters. Cluster C3 (red) was expanded during progression to clinical RA. (F) Mean RNA expression of select genes across CD4mem clusters. Cluster C3 (red) was expanded during progression to clinical RA. (G) Normalized RNA expression of select genes that promote differentiation to B cell helper and Th17 cells in CD4mem-C3 (red) vs. remaining CD4mem clusters (blue). (H) Differentially expressed genes between CD4mem-C3 (red) and remaining CD4mem clusters (blue). Select genes associated with T cell activation are labeled. Violin plot shape is proportional to the density of data points (dots). Larger width represents a higher density of data points. P values were determined by linear mixed models (A, C), the Kruskal-Wallis test with pairwise Dunn’s posthoc test (E) or the Wilcoxon rank-sum test (G-H). Nominal P values are indicated for (B, E). FDR values are indicated for all other plots. Values below 0.1 were considered significant.

Fig. 5:. Naïve CD4 T cells in ARI have a multimodal signature of activation.

(A) PBMC TEA-seq was performed on a subset of ARI and HC2 samples. (B) Percentage of variance in each modality (surface protein, plasma protein, RNA, ATAC) explained by Multi-Omics Factor Analysis (MOFA) factors. (C) Factor 1 scores comparing ARI and HC2. (D) Scaled normalized expression of select genes in Calcium–Calcineurin–NFAT pathway in ARI and HC2. (E) Inferred accessibility for the top 15 transcription factors (TFs) positively or negatively associated with factor 1, ranked by weight. (F) RNA expression differences in core naïve CD4 T cells over time in CONV (orange) who progress to clinical RA (purple). Genes associated with T cell activation are annotated. (G) T cell RNA activation metric in core naïve CD4 T cells over time as CONV progress to clinical RA. Each participant’s longitudinal series is connected by a gray line, with a group trendline and 95% confidence interval in purple. (H) Enriched pathways associated with T cell activation in core naïve CD4 T cells over time as CONV progress to clinical RA. P values were determined by Wilcoxon rank-sum test and FDR were adjusted for all MOFA factors tested (C). P value was determined by linear mixed models (F, G). Normalized enrichment scores (NES) and adjusted P values, by GSEA, are shown (H). Nominal P value is shown in (G). FDR values are indicated for remaining plots. Values below 0.1 were considered significant.

Fig. 6:. Naïve CD4 T cells exhibit a bias toward B cell help in ARI.

(A) Louvain clusters in CD4 T cells by ATAC modality in TEA-seq. (B) Centered log-ratio (CLR)-transformed frequencies of ATAC clusters CD4nve-T2 and CD4nve-T5 in CD4 T cells. (C) Mean surface protein expression of select markers differentiating CD4 naïve, memory, regulatory T (T_reg), and cytotoxic CD4 T cells (CTL) across ATAC clusters. (D) ATAC UMAP overlaid with inferred gene activity scores for CXCR5 and IL21. (E) ChromVAR TF activity Z scores of BCL6 and STAT3 in CD4 T cells. (F) ATAC signal in ARI (orange) versus HC2 (green), and the delta between the two (red) at the IL21 locus. The gray box highlights a 500bp region containing differentially accessible peaks between ARI and HC2 (chr4: 122,617,500–122,617,999). Black arrows indicate the motif locations of BCL6 and STAT3 binding sites. Gene bodies are displayed on the bottom. Boxplots show median (centerline) and first and third quartiles (lower and upper bound of the box); whiskers show the 1.5x interquartile range of data. P values were determined by linear models (B) or by the zero-inflated Wilcoxon test (F). FDR values are indicated and values below 0.1 were considered significant.

Fig 7:. Gene signatures in CD4 T cells of converters reflect Abatacept treatment response.

Longitudinal DEGs as CONV progress to clinical RA were assessed within the context of RA patients with efficacious (responders) or non-efficacious (non-responders) clinical response to abatacept (ABT) or TNF inhibitor (TNFi) treatment (from (47, 48)). (A) Overview of the analysis strategy. mo., month. (B and C) Over-representation of CONV cell type-specific longitudinal DEGs amongst ABT (B) or TNFi (C) response DEGs. (D and E) Significant DEGs in ABT (D) or TNFi (E) responders compared to DEG changes over time as CONV progress to clinical RA in core naïve CD4 T cells (left) and CM CD4 T cells (right). Genes (dots) previously implicated in RA-like disease are labeled. (F and G) Normalized RNA expression of NABP1 over time as CONV progress to clinical RA (F) and pre- vs. post-ABT therapy in patients with RA (G). In (F), each participant’s longitudinal series is connected by a gray line, with a group trendline and 95% confidence interval in purple. In (G), each participant’s two samples are connected by a gray line. (H) Odds ratios of the number of longitudinal DEGs in core naïve and CM CD4 T cells from CONV that were reversed by ABT or TNFi treatment. Error bars indicate confidence intervals. P values were determined by hypergeometric enrichment tests (B, C), McNemar’s Chi-squared test (D, E), linear mixed models (F), or Wald test (G). Nominal P values are indicated for (C, D). FDR values are indicated for other panels. Odds ratios were determined by unconditional maximum likelihood estimation method and Z-test was used to compare odds ratios between Abatacept and TNFi treatment (H). Values below 0.1 were considered significant.