Molecular Portraits of Early Rheumatoid Arthritis Identify Clinical and Treatment Response Phenotypes

Abstract

There is a current imperative to unravel the hierarchy of molecular pathways that drive the transition of early to established disease in rheumatoid arthritis (RA). Herein, we report a comprehensive RNA sequencing analysis of the molecular pathways that drive early RA progression in the disease tissue (synovium), comparing matched peripheral blood RNA-seq in a large cohort of early treatment-naive patients, namely, the Pathobiology of Early Arthritis Cohort (PEAC). We developed a data exploration website (https://peac.hpc.qmul.ac.uk/) to dissect gene signatures across synovial and blood compartments, integrated with deep phenotypic profiling. We identified transcriptional subgroups in synovium linked to three distinct pathotypes: fibroblastic pauci-immune pathotype, macrophage-rich diffuse-myeloid pathotype, and a lympho-myeloid pathotype characterized by infiltration of lymphocytes and myeloid cells. This is suggestive of divergent pathogenic pathways or activation disease states. Pro-myeloid inflammatory synovial gene signatures correlated with clinical response to initial drug therapy, whereas plasma cell genes identified a poor prognosis subgroup with progressive structural damage.

Keywords: PEAC; Pathobiology of Early Arthritis Cohort study; RNA sequencing; ectopic lymphoid structures; lymphoid neogenesis; personalized medicine; rheumatoid arthritis; synovial biopsy; transcriptomics.

Graphical abstract

Figure 1

Synovium RNA Sequencing Correlates with Histological Pathotype in Early Rheumatoid Arthritis (A) Immunohistochemistry of synovial biopsies for CD20⁺ B cells, CD3⁺ T cells, and CD68⁺ macrophages in synovial lining or sublining layers and CD138⁺ plasma cells from treatment-naive individuals with early rheumatoid arthritis. Synovial biopsies were categorized as lympho-myeloid (B cell aggregates present), diffuse-myeloid (sublining macrophage infiltration), or pauci-immune fibroid (lack of or low inflammatory cell infiltrate). (B) Comparison of cell-specific RNA-seq gene module scores with histology scores. (C) Cell-specific gene scores compared across histology pathotypes. Statistical analysis by one-way ANOVA with Bonferroni post-test. (D) Clustering of lympho-myeloid, diffuse-myeloid and pauci-immune fibroid samples according to B cell, monocyte, and synoviocyte RNA-seq modules. (E) Heatmap showing hierarchical clustering of cell-specific gene module scores and collapsed module space (right) highlighting cellular composition of synovial biopsies from each pathotype.

Figure 2

Clinico-radiographic Correlates of Cell-Specific Gene Modules in Rheumatoid Arthritis Synovium (A) Correlation heatmap showing Spearman correlation of cell-specific gene modules against baseline clinical (ESR, erythrocyte sedimentation rate; CRP, C-reactive protein; CCP, anti-cyclic citrullinated peptide antibody titer; RF, rheumatoid factor titer; VAS, visual analog score; HAQ, health assessment questionnaire), ultrasonographic scores (ST, synovial thickness; PD, power doppler) at the biopsy joint (Ultrasound ST/PD BJ) or across 12 representative joints (Ultrasound ST/PD 12) and radiographic parameters (Total Sharp van der Heijde score). (B) Boxplots of clinical parameters by tertile demonstrating correlation with cell-specific gene modules. (C) Linear regression of ultrasound biopsy joint parameters against cell-specific gene modules. (D) Boxplots of total Sharp van der Heijde radiographic score by tertile correlated with cell-specific gene modules. p values were calculated by linear regression models.

Figure 3

Synovium and Blood RNA-Seq Comparison Reveals Differential Axes of Gene Expression (A and B) 3D cylindrical volcano plots of differentially expressed genes comparing RNA sequencing of (A) synovial tissue and (B) whole blood. Vectors for pathotype mean Z score per gene were projected onto a polar coordinate space analogous to RGB (red-green-blue) color space mapping to HSV (hue-saturation-value) as described in the STAR Methods. Lympho-myeloid, diffuse-myeloid, and pauci-immune fibroid vectors are mapped to 3 axes lympho-myeloid (L), diffuse-myeloid (M), and pauci-immune fibroid (F) using polar coordinates in the horizontal plane. The z axis shows –log₁₀ p value for likelihood ratio test. Genes with adjusted p value for likelihood ratio test < 0.05 (z axis) were considered significant (non-significant genes colored gray). Colors demonstrate pairwise comparisons (FDR < 0.05) between the 3 histological pathotypes: primary colors denote upregulation in one group only for lympho-myeloid (blue), diffuse-myeloid (red), and pauci-immune fibroid (green) compared to reference group with minimum gene expression; composite colors show genes significantly upregulated in two groups (myeloid+lymphoid, purple; fibroid+myeloid, yellow; lymphoid+fibroid, cyan). Lateral view and 2D polar plots are shown below. (C and D) Principal-component analysis of whole transcriptome RNA-seq data from untreated rheumatoid arthritis (C) synovium and (D) whole blood, showing separation of lympho-myeloid (blue) and pauci-immune fibroid (green) histological pathotypes on principal component 1 (PC1) for synovial RNA-seq, with separation of diffuse-myeloid (red) and pauci-immune fibroid samples on PC1 in whole blood.

Figure 4

Clustering and Pathway Analysis of Differentially Expressed Genes in Rheumatoid Arthritis Synovium (A) Heatmap of 2,964 RNA-seq genes differentially expressed between three histological pathotypes (lympho-myeloid, diffuse-myeloid, and pauci-immune fibroid) (FDR < 0.05, n = 87). Upper tracks show histological scores for CD3, CD20, CD68L, CD68SL, and CD138 and overall pathotype. Unsupervised hierarchical clustering demonstrated clustering of genes into four clusters, demonstrating some overlap between the three histologically determined pathotypes. (B) Ingenuity Pathway Analysis performed on synovial gene clusters produced by hierarchical clustering identified pathways by gene enrichment, using whole genome as background. Clusters S1 and S2 represent pauci-immune fibroid and diffuse-myeloid samples, and clusters S3 and S4 represent lympho-myeloid and diffuse-myeloid samples. Color scale and numbers depict –log₁₀ FDR-adjusted p values.

Figure 5

Modular Analysis of Synovial and Blood RNA Sequencing (A–C) Three axis polar plots of synovium (A) and blood (B) gene modules based on blood microarray modules (Li et al., 2014) and (C) synovium modules derived by weighted correlation network analysis (WGCNA), analyzed using QuSAGE. Modules are color-coded for statistical significance (FDR < 0.05) for upregulation in different pathotypes. WGCNA synovium modules were annotated against single-cell RNA-seq cell types (Stephenson et al., 2018). (D) Heatmap showing gene expression in selected gene modules in synovium and blood, grouped by pathotype. (E) Boxplots of summarized module scores in synovium and blood. Statistical analysis by QuSAGE with FDR correction: ^∗FDR < 0.05, ^∗∗FDR < 0.01, ^∗∗∗FDR < 0.001. (F) Comparison of number of significant (FDR < 0.05) synovium and blood gene modules (Li et al., 2014), which correlate with clinical and radiographic markers of disease activity and response to 6 months treatment with DMARDs, in either synovium only (blue) or blood only (red) or are concordantly correlated in both compartments (purple). Statistical analysis by Spearman correlation. (G) Bubble plot of –log₁₀ p values comparing correlation of gene modules in synovium and blood with disease activity measured by DAS28-CRP, showing significantly (FDR < 0.05) correlated modules found in synovium alone (blue), blood alone (red), or concordantly correlated with DAS28-CRP in both synovium and blood (purple).

Figure 6

Baseline Synovium Plasma Cell Gene Expression Is Associated with CCP Antibody Positivity and Worse Prognosis at 12 Months (A and B) Differential gene expression in synovium RNA-seq comparing (A) anti-CCP antibody (ACPA)-positive and ACPA-negative RA individuals and (B) individuals with progression of bone erosions on X-rays at 12 months compared to baseline versus non-progressors. (C and D) Single-cell RNA-seq-annotated WGCNA modular analysis shows that increased plasma cell module expression is associated with ACPA positivity (C) and is predictive of bone erosion radiographic progression at 12 months (D). (E) Upstream regulator analysis using Ingenuity Pathway Analysis showing upstream regulator effects of cytokines and chemokines associated with ectopic lymphoid structure development in synovium. (F and G) Sankey diagrams showing change in histological pathotype following 6-month treatment with disease-modifying anti-rheumatic drugs (DMARDs) for (F) whole cohort or (G) grouped by ACPA status. Statistical analysis by Fisher’s test. (H) Shift in pathotype between baseline and 6 months correlated against change in DAS28-ESR. Statistical analysis by Pearson correlation.

Figure 7

Association of RNA-Seq Modules and Response to 6 Months of DMARD Treatment (A and B) Correlation of gene modules in synovium (x axis) versus blood (y axis) with 6-month response to DMARD treatment measured by (A) delta DAS28-CRP and (B) delta ESR. Significantly correlated synovial modules (at FDR < 0.05) are shown in blue, significant blood modules in red, and modules, which concordantly correlate with each clinical parameter in both synovium and blood, are shown in purple. (C) Correlation of synovium gene modules with delta DAS28-CRP from baseline to 6 months. (D and E) Correlation of (D) synovium and (E) blood gene modules against change in ESR from baseline to 6 months following DMARD treatment. Statistical analysis by Spearman correlation with FDR adjustment (A–E). (F) Differential expression of synovial single-cell-annotated WGCNA modules between EULAR DAS28-CRP responders (good and moderate) and non-responders. Statistical analysis by QuSAGE with FDR adjustment.