Disease activity drives transcriptomic heterogeneity in early untreated rheumatoid synovitis

Abstract

Objectives: Transcriptomic profiling of synovial tissue from patients with early, untreated rheumatoid arthritis (RA) was used to explore the ability of unbiased, data-driven approaches to define clinically relevant subgroups.

Methods: RNASeq was performed on 74 samples, with disease activity data collected at inclusion. Principal components analysis (PCA) and unsupervised clustering were used to define patient clusters based on expression of the most variable genes, followed by pathway analysis and inference of relative abundance of immune cell subsets. Histological assessment and multiplex immunofluorescence (for CD45, CD68, CD206) were performed on paraffin sections.

Results: PCA on expression of the (n=894) most variable genes across this series did not divide samples into distinct groups, instead yielding a continuum correlated with baseline disease activity. Two patient clusters (PtC1, n=52; PtC2, n=22) were defined based on expression of these genes. PtC1, with significantly higher disease activity and probability of response to methotrexate therapy, showed upregulation of immune system genes; PtC2 showed upregulation of lipid metabolism genes, described to characterise tissue resident or M2-like macrophages. In keeping with these data, M2-like:M1-like macrophage ratios were inversely correlated with disease activity scores and were associated with lower synovial immune infiltration and the presence of thinner, M2-like macrophage-rich synovial lining layers.

Conclusion: In this large series of early, untreated RA, we show that the synovial transcriptome closely mirrors clinical disease activity and correlates with synovial inflammation. Intriguingly, lower inflammation and disease activity are associated with higher ratios of M2:M1 macrophages, particularly striking in the synovial lining layer. This may point to a protective role for tissue resident macrophages in RA.

Keywords: arthritis, rheumatoid; methotrexate; synovitis.

Figure 1

Transcriptomic heterogeneity of synovia from patients with untreated, early rheumatoid arthritis (ERA). (A) Principal components analysis (PCA) plot showing distribution of the 74 samples based on the expression of n=894 most variable genes. Samples are coloured by disease activity score (DAS28CRP colour-scale). Samples in patient cluster 2 (PtC2) from (B) are circled. (B) Heatmap of scaled, normalised expression of the 894 most variable genes across samples. Samples are classified into two patient clusters (PtC): PtC1 (blue) and PtC2 (yellow), using unsupervised clustering (average linking rule). Gender, DAS28CRP, presence of bone erosions at baseline and ACPA seropositivity are indicated above the heatmap. Genes fall into three groups (gene clusters: GC) using unsupervised clustering (average linking rule): GC1 (pink), GC2 (purple) and GC3 (orange); in white: genes located on the Y chromosome, expressed only in male patients across both PtC. These are excluded from pathway analyses. (C) Pathway enrichment analysis using Metascape (top three based on FDR-adjusted q-value) of GC 1–3, coloured as in (B). ACPA, anti-citrullinated peptide antibody; DAS28CRP, Disease Activity Score 28-joint count C reactive protein; FDR, false discovery rate.

Figure 2

Clinical characteristics and response to treatment of patients in PtC1 and PtC2. (A) Baseline clinical characteristics of patients in PtC1 and PtC2. Values shown are Q50 (Q25-Q75) or percentages. P value of Mann-Whitney or Fisher’s exact test are shown. *Data unavailable for one patient from PtC1; assigning this patient to one or the other group (steroid/no steroid) does not modify p value. (B, C) Effect of therapy on DAS28CRP in methotrexate-treated patients from PtC1 (n=35) (B) and PtC2 (n=15) (C); p values (paired Wilcoxon) shown. (D) Change in DAS28CRP from baseline to 3 months (∆DAS M3-BL) in PtC1 and PtC2; p value (Mann-Whitney) shown. (E) Correlation between baseline (BL) DAS28CRP and delta DAS28CRP from baseline to 3 months (∆DAS M3-BL) in n=50 patients treated with methotrexate (data unavailable at M3 for one patient in PtC1). Dots coloured according to PtC; p value and Spearman’s r are shown. ACPA, anti-citrullinated peptide antibodies; HAQ, Health Assessment Questionnaire; PhGA, Physician Global Assessment; PtGA, Patient Global Assessment; RIN, RNA integrity number; RF, rheumatoid factor; SJC28, Swollen Joint Count across 28 joints; TJC28, Tender Joint Count across 28 joints.

Figure 3

Transcriptomic and histological differences between PtC1 and PtC2. (A) Volcano plot of all differentially expressed genes (no restriction on variance) between PtC1 and PtC2 using t-test. Dashed line: p-value threshold=0.01. Genes belonging to GC 1 (pink), 2 (purple) and 3 (orange) coloured as in figure 1B. Gene symbols of the top 10 upregulated genes in each PtC (based on t-test p values) are shown. (B, C) Pathway enrichment analysis on top 500 upregulated genes based on t-test (with FDR-adjusted q-values), in PtC1 (B) and PtC2 (C). (D–G) Scoring for features of synovial inflammation on H&E sections: lining layer hyperplasia (D), inflammatory infiltrate (E), fibrinoid necrosis (F) and vascular proliferation (G). P values of Mann-Whitney are shown on the graph. FDR, false discovery rate; PtC, patient cluster.

Figure 4

M2:M1 ratios are inversely correlated with clinical disease activity scores. (A) CibersortX-inferred immune cell proportions in samples from PtC1 and PtC2 (mean proportion±SD of each immune cell subset). *****p<0.000001, ****p<0.00001, ***p<0.0001, **p<0.001, *p<0.01 (unpaired t-tests with 5% FDR (Benjamini-Yekutieli)). (B–E) Correlation between M2/M1_RNA (M2:M1 ratio by CibersortX on RNASeq) and baseline clinical disease activity scores (n=74 samples): DAS28CRP (B), Tender Joint Count (C), Swollen Joint Count (D) and CRP (E). P value and Spearman’s r are shown. DAS28CRP, Disease Activity Score 28-joint count C reactive protein; FDR, false discovery rate; PtC, patient cluster.

Figure 5

M2:M1 ratios are associated with differences in synovial inflammation, and synovial lining macrophage composition and thickness. (A) Correlation between M2/M1 ratio by CibersortX on RNASeq (M2/M1_RNA), and digital quantification of multiplex immunofluorescence (number of CD45⁺CD68⁺CD206⁺ cells/CD45⁺CD68⁺CD206^- cells, across entire biopsies) (M2/M1_IF). Spearman’s r and p value are shown. (B) Correlation between M2/M1_RNA ratio and immune cell infiltration based on digital quantification of multiplex immunofluorescence (number of CD45⁺ cells/total cells (nuclei) across entire biopsies) (CD45⁺ fraction_IF). Spearman’s r and p value are shown. (C) Representative examples of multiplex immunofluorescence images from samples with M1-predominant lining (M1L, left panels) and samples with M2-predominant lining (M2L, right panels). Nuclei in blue, CD68 in green (top panels) and CD206 in orange (bottom panels). (D) M2/M1_IF ratio (number of CD45⁺CD68⁺CD206⁺ cells/CD45⁺CD68⁺CD206^- cells, across entire biopsies) in M1L and M2L classified samples. Boxplots showing mean and IQR. P value of Mann-Whitney is shown. (E) Semiquantitative scores of CD68⁺ lining of biopsies by immunofluorescence (lining thickness_IF) in M1L and M2L classified samples. Bars represent mean and SD. P value of Mann-Whitney is shown.