Introduction

Synovitis is the clinical hallmark of rheumatoid arthritis (RA), and the analysis of synovial tissue from affected joints provides unique insights into mechanisms driving synovial inflammation and disease activity. Classically described as the seat of immune cell infiltration and resident cell hyperplasia, the RA synovium is a highly heterogeneous tissue. Histological as well as bulk transcriptomic studies also show wide variability of features, including the extent and composition of immune cell infiltration, between patients.1 Extra-synovial factors such as disease duration and treatment contribute to variability in synovial signals.1 Recent studies have however confirmed wide variations in immune alterations even among patients with untreated, early rheumatoid arthritis (ERA).2 3 This has led to a keen interest in exploring the potential of ERA synovial tissue analysis in defining clinically relevant patient subgroups; the much sought-after goal is effective patient stratification for truly tailored (type and strength of) therapy.4 5 Different patterns of synovial inflammation associated with baseline clinical features have therefore been explored for their ability to earmark different disease mechanisms and/or treatment responses. A current paradigm distinguishes between lympho-myeloid, diffuse-myeloid and pauci-immune subgroups (pathotypes), defined based mainly on the relative abundance in synovia of T cells, B cells and macrophages (as assessed by histology and molecular signature), and suggested to predict disease evolution and treatment response.2

Somewhat in contrast with these observations, we previously reported that different disease-modifying antirheumatic drugs (DMARDs) downregulate a robust set of genes in common, notwithstanding their direct targets or mechanisms of action.6 This would suggest a high degree of interdependence between the different cell populations and processes in the synovium that therefore exert knock-on effects on each other. The degree of downregulation of these ‘shared’ genes seemed to depend primarily on their pretreatment expression levels, and we could not discern distinct signatures significantly associated with clinical response to any particular DMARD.6 The patients in these studies had longstanding disease and undergone multiple, heterogeneous rounds of therapy leading up to the specific DMARD being tested, begging the question of whether biopsies performed early in disease would be more informative. Here, we apply unsupervised clustering to RNASeq data from a large series of ERA synovial biopsies, for a data-driven assessment of (1) the most variable transcriptomic features among patients with untreated RA at early disease stage, (2) their baseline clinical correlates and (3) their ability to define clinically relevant patient subgroups according to response to methotrexate.

Materials and methodsPatients and samples

Patients with active, early (<12 months disease duration) rheumatoid arthritis, naïve to DMARD and steroids, were recruited for this study. All patients met the American College of Rheumatology/European Alliance of Associations for Rheumatology (ACR/EULAR) 2010 RA classification criteria. At inclusion, patients were assessed for disease activity and presence of bone erosions on radiographs. Processing and analysis of biopsies was performed in accordance with EULAR minimal reporting requirements for synovial tissue research.7 Synovial tissue from a clinically affected joint was obtained using needle arthroscopy or ultrasound (US)-guided biopsy prior to starting treatment. Three to five synovial fragments were obtained from each joint by knee needle arthroscopy or US-guided small joint biopsy and stored at −80°C after overnight incubation in RNA-Later solution (Invitrogen). Where possible, an additional three to five fragments from the same biopsy were fixed overnight in 10% formalin buffer at pH 7.0 and embedded in paraffin for histology and immunofluorescence (n=36 from knees, 2 from wrists, 1 from metatarsophalangeal joint (MTP), 1 from metacarpophalangeal joint (MCP)). Patients were started on conventional treatment (methotrexate±steroids as bridging therapy, n=51) or on other, heterogeneous treatments (n=23) (table 1) (online supplemental table 1). Effects of therapy on disease activity in methotrexate-treated patients was assessed after 3 months.

Table 1

Baseline clinical characteristics of n=74 patients with RA included in the study

Patient and public involvement

Patients and the public were not involved in the design, conduct or reporting of this research.

Transcriptomic profiling

Total RNA was extracted from synovial biopsies using the NucleoSpin RNA II extraction kit (Macherey-Nagel) including DNAse treatment. RNA quality control was performed on a 2100 Bioanalyser (Agilent). RNASeq was performed on 350 ng of RNA, for all (n=74) samples that yielded a RNA Integrity Number (RIN) of 5 or more. Sequencing was performed by Macrogen Europe using a standardised protocol. Briefly, cDNA libraries were constructed with the TruSeq Stranded mRNA LT Sample Prep Kit and sequenced on an Illumina sequencer. Data were processed using a standard RNAseq analysis pipeline including FastQC-0.11.8 for quality control of raw data, and Trimmomatic-0.38 to remove low quality reads and trim adapters from sequences. Hisat2-2.1.0 was used to align reads to the genome (GRCh38).8 Gene expression levels were evaluated using featureCounts from subread-2.0.09 and Homo_sapiens.GRCh38.94.gtf. Sequencing data were deposited in the Gene Expression Omnibus of the National Centre for Biotechnology Information and are accessible through Gene Expression Omnibus accession number GSE225731.

Histological evaluation

Degree of lining layer hyperplasia, inflammatory infiltrate, vascular proliferation and fibrinoid necrosis was graded on a 0 (absent) to 3 (maximum) scale on all H&E-stained sections from (n=39 of) the same biopsies used for RNASeq, by an expert pathologist (CG) blinded to clinical and transcriptomic data.

Multiplex immunofluorescence (IF)

Multiplex IF was performed on 5 µm serial sections from the same (n=39) biopsies as above. After deparaffinisation, tissue sections were processed according to the protocol described in the literature.10 Briefly, endogenous peroxidases were inhibited with 3% hydrogen peroxide in methanol for 20 min. Antigen retrieval was performed in 10 mM citrate buffer pH 5.7 and blocking in Tris buffered saline (TBS)−5% bovine serum albumin (BSA)−0.1% Tween20. Sections were incubated in the primary antibody overnight at 4°C (in TBS−1% BSA−0.1% Tween20), followed by the corresponding horseradish peroxidase-conjugated secondary antibody (Envision, Dako) for 1 hour at room temperature, and visualised by tyramide signal amplification using AlexaFluor (AF)-conjugated tyramides (Thermo Fisher Scientific). Antibodies were stripped by incubating in citrate buffer. This protocol was applied sequentially with each of the primary and secondary antibodies, with different AF-tyramides, for a total of three markers: CD45 (1:200, Dako M0701) with AF555 (cat. no. B40955), CD68 (1:500, Dako M0876) with AF488 (cat. no. B40953), CD206 (1:2500, Cell Signaling 91992) with AF647 (cat. no. B40958). After a final wash in phosphate buffered saline (PBS), nuclei were stained with Hoechst 33 342 (Thermo Fisher Scientific; in TBS−10% BSA−0.1% Tween 20) and washed in TBS−0.1% Tween20. Sections were mounted with fluorescence mounting medium (Dako). Multispectral image acquisition was performed using an Axioscan.z1 (Zeiss, ×20 objective).

Quantitative and semiquantitative evaluation of multiplex IF

Entire multiplex-stained paraffin sections were quantified using the image analysis tool Author V.2022.01.4 (Visiopharm). Tissue sections were first automatically delineated at low digital magnification (×5). Following cell detection and segmentation at high resolution (×20, based on the Hoechst staining), postprocessing steps were applied to separate unstained and stained cells for each marker (CD45, CD68 and CD206). Results were expressed in terms of proportion of marker-positive cells/total cells (nuclei), or as ratios of different marker-positive cells. Classification of biopsies based on the macrophage composition of the synovial lining, as well as evaluation of synovial lining thickness (visualised using CD68, on a scale of 1 to 3), was performed on Zen lite image analysis software (Zeiss) by two scorers blinded to transcriptomic and clinical data. Samples were classified as having a (CD45+CD68+CD206-) M1-predominant lining (M1L) or a (CD45+CD68+CD206+) M2-predominant lining (M2L) (see online supplemental figure 1A–D).

Statistical analyses

Raw RNASeq read counts were normalised using DESeq2 Bioconductor Package.11 Normalised expression data were log2-transformed, scaled to a mean of 0 and variance of 1, and imported into Qlucore Omics Explorer for further analysis. We filtered for genes with variance (σ/σmax) ≥0.4. Principal component analysis (PCA) and linear regression analysis were performed and samples and genes were classified into clusters based on the gene expression data using unsupervised agglomerative hierarchical clustering. Differentially expressed genes between patient clusters were calculated using t-tests with correction for multiple testing using false discovery rate (FDR)-adjusted q-values. Pathway enrichment analyses (top 3 with FDR-adjusted q-values (figure 1), top 20 with p values (online supplemental figures)) were performed on Metascape.12 CibersortX13 was used to extrapolate relative abundance of immune cell subsets from gene expression data (using the LM22 reference dataset). Other statistical analyses (Spearman’s correlation coefficient, Mann-Whitney, t-tests with Benjamini-Yekutieli correction for multiple testing) were performed on GraphPad Prism V.9.

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.

ResultsHeterogeneity of synovial gene expression profile in ERA

RNASeq was performed on synovial biopsies from 74 patients with ERA. Patient and sample characteristics are shown in table 1 and online supplemental figure 2A–G. We applied variance filtering to identify the genes whose expression showed the greatest variation across the cohort (σ/σmax≥0.4, n=894 most variable genes) (figure 1A,B). PCA based on the expression of these genes yielded a continuous distribution (rather than distinct clusters) of samples, broadly associated with baseline DAS28CRP (figure 1A). Unsupervised clustering was used to demark two patient clusters (PtC): PtC1 (n=52) and PtC2 (n=22) (figure 1B), based on their expression of these 894 genes. Pathway analysis of these genes revealed significant enrichment of immune system genes, in the inflammatory response and rheumatoid arthritis pathways (gene cluster 1: GC1) and B cell and plasma cell-related pathways (GC2) (figure 1C, online supplemental figure 3A,B, online supplemental table 2), upregulated in PtC1. Almost 50% of genes in GC2 are immunoglobulin (Ig) genes, potentially overwhelming signals from additional pathways. We therefore performed pathway analysis after removing all Ig genes, and found that B cell activation remains the most significantly enriched (q-value=10-2.3, not shown). The B cell transcriptomic signature in PtC1 was also confirmed using a validated synovial B cell gene set (online supplemental figure 3D).14 The third cluster of genes (GC3, figure 1C, online supplemental figure 3C, online supplemental table 2), related to fatty acid metabolism and PPAR signalling, is instead upregulated in PtC2. These pathways have been associated with M2-like (in vitro) or non-conventional tissue resident macrophages (in vivo), described to have anti-inflammatory, pro-homeostatic functions.15 16

Unbiased ERA synovial transcriptome-driven clustering defines two clinically meaningful groups

Having used the most variably expressed genes for unbiased clustering of samples, we set out to explore the differences between the two patient groups thus defined. Interestingly, PtC1 and PtC2 were characterised by distinct clinical features (figures 1B and 2A, and online supplemental figure 4A–G): Patients in PtC1 displayed higher baseline indices of all disease activity score components (resulting in a median DAS28CRP of 5.56 in PtC1 vs 4.09 in PtC2; p value=0.0003) except Patient Global Assessment. They also had higher rates of baseline radiological erosions (erosive disease in 34.6% vs 10%; p value=0.0252) (figure 2A). Notably, ACPA/RF positivity was not significantly different between the groups.

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.

Baseline SJC28 was the most significantly distinct clinical value between PtC1 and PtC2 (online supplemental figure 4A). Genes whose expression correlated (n=572; q-value <0.005 using regression analysis) with SJC28 were highly enriched in lymphocyte activation and adaptive immune response pathways (online supplemental figure 5 and online supplemental table 3). We previously showed that these pathways are downregulated by multiple DMARDs, including methotrexate.6 We therefore assessed for whether PtC1 and PtC2 differ in terms of response to first-line methotrexate treatment (used in n=51 patients) at 3 months. A significant effect of methotrexate was indeed observed in PtC1 but not PtC2 (figure 2B–D). As previously observed by others,17 however, treatment response is correlated with baseline disease activity, whether in terms of absolute (figure 2E) or relative drop (not shown) in DAS28CRP (r=−0.67, p<0.0001 and r=−0.55, p<0.0001, respectively). Notably, the relationship between baseline disease activity and drop in DAS28CRP was similar in PtC1 and PtC2 (Spearman’s r=−0.621 and −0.593, respectively). This confirms that the synovial transcriptome in early, untreated RA, like that in later-stage, heterogeneously treated disease, is strongly tied to disease activity, a major determinant of treatment response amplitude.

Distinct pathways, histological features and immune cell signatures distinguish the two synovial transcriptome-based ERA groups

We found that the PtC1 and PtC2 groups have significant, clinically relevant characteristics. We therefore set out to perform a detailed supervised analysis of the differentially expressed genes between these two groups (figure 3A). Genes upregulated in PtC1 synovia genes (top 500, ranked by FDR-adjusted q-value) show significant enrichment of pathways related to B cell biology, as well as other immune pathways with critical roles in RA such as lymphocyte activation, TCR signalling and cytokine signalling (figure 3B, online supplemental figure 6A and online supplemental table 4). By contrast, PtC2 is characterised by the upregulation of genes (top 500, ranked by FDR-adjusted q-value) involved in lipid metabolism-related pathways, although with lower statistical support for pathway enrichment (figure 3C, online supplemental figure 6B and online supplemental table 4). In keeping with the upregulation of immunity/inflammation-related genes in PtC1 synovia, multiple histological features of inflammation were found to be higher in PtC1 synovia than in PtC2 (figure 3D–G), with fibrinoid necrosis being the most pronounced.

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.

M2-like:M1-like macrophage ratio is associated with lower disease activity and synovial inflammation

We next used CibersortX to extrapolate information on the relative proportions of immune cell subpopulations in synovia, from gene expression data.13 In line with pathway analyses, PtC1 had a higher transcriptomic signature (inferred relative frequency) of B and plasma cells, T cells and M1-like macrophages than PtC2. PtC2 instead had relatively higher M2-like macrophage and resting mast cell signatures (figure 4A). Accordingly, the M2-like macrophage signature was inversely correlated with the B cell, T cell and M1-like macrophage signatures (online supplemental figure 7), as well as with clinical indices of baseline disease activity (figure 4B–E and online supplemental figure 7). M2-polarised macrophages have been shown to have a metabolic profile characterised by predominant fatty acid oxidation,18–20 and we indeed found that the M2-like macrophage signature correlated with expression of GC3 genes across samples (online supplemental figure 8A,B).

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/M1RNA (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.

We used multiplex IF (for CD45, a pan-immune cell marker; CD68, a pan-macrophage marker; CD206, an M2-like macrophage marker) to confirm the inverted M2/M1 ratios suggested to distinguish PtC1 and PtC2 samples by CibersortX. Digitally quantified IF data correlated well with transcriptomic data (M2/M1 ratios, figure 5A), and confirmed the inverse correlation between immune infiltration of the synovium (CD45+ fraction by IF) and M2/M1 ratio (figure 5B). Semiquantitative evaluation of multiplex IF added a spatial dimension to our analysis, revealing marked variations in M2 versus M1 content of synovial lining among samples. Samples were classified as M1L (M1-Lining: predominantly CD206- lining) or M2L (M2-Lining: predominantly CD206+ lining) (see the Materials and methods section, figure 5C, online supplemental figure 1A–D). M2L samples had higher M2/M1IF ratios than M1L samples (figure 5D); they also showed significantly lower (CD68+) synovial lining thickness on IF (figure 5E).

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/M1RNA), and digital quantification of multiplex immunofluorescence (number of CD45+CD68+CD206+ cells/CD45+CD68+CD206- cells, across entire biopsies) (M2/M1IF). Spearman’s r and p value are shown. (B) Correlation between M2/M1RNA ratio and immune cell infiltration based on digital quantification of multiplex immunofluorescence (number of CD45+ cells/total cells (nuclei) across entire biopsies) (CD45+ fractionIF). 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/M1IF 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 thicknessIF) in M1L and M2L classified samples. Bars represent mean and SD. P value of Mann-Whitney is shown.

Discussion

Histological and transcriptomic analyses of synovial biopsies have provided invaluable insights into tissue inflammatory processes in RA. Most studies have focused on established disease6 21–23 or have used a priori classification approaches.2 4 24 Here, we instead exploited the transcriptional heterogeneity of 74 DMARD-naïve ERA synovial biopsies for unbiased, data-driven clustering of samples. The two groups thus defined (PtC1 and PtC2) showed distinct clinical and histological features, suggesting biologically meaningful clustering. To the question of whether early, treatment-naïve transcriptomes may allow for patient stratification or prediction of treatment response: PtC1 showed a greater decrease in (absolute and relative, ie, baseline-normalised) disease activity 3 months into treatment with methotrexate. Given that the synovial transcriptome correlates strongly with baseline disease activity, we cannot yet conclude the former is an independent (or better) predictor of treatment response.

That the synovial transcriptome captures local (histological) inflammation is not surprising. The degree to which it echoes disease extent (as measured in clinical scores), while not unexpected, is striking. It is also consistent across studies.2 3 Increased numbers of tender and swollen joints no doubt tend to be linked to greater local inflammation; nevertheless, the association between synovial transcriptomic profile and global disease activity indices speaks to their robustness and continuing utility. While not designed to specifically address the question, we assessed for previously described pathotype-signatures in our series.2 We did not find evidence for distinct (lymphoid, myeloid and fibroid) subgroups.2 4 24 This must be interpreted with some caution: The two studies differ significantly in design and analysis; in particular, unbiased versus hypothesis-driven, reflected in the use of whole transcriptome data (RNASeq) versus a selected gene-panel (NanoString). The former likely entails a far larger cohort for pathotype-defining genes to reach genome-wide significance. Hierarchical clustering based on the expression of published pathotype-gene sets however also resulted in a similar continuum of samples as in our analysis: from fibroid with lower disease activity, to high inflammatory (lymphoid+myeloid) with higher disease activity, without evident lymphoid vs myeloid clustering (online supplemental figure 9A–C).

Based on the genes differentially expressed between the two groups, synovial biopsies from patients in PtC1 have strong signatures of inflammation and adaptive (particularly B cell lineage) immunity. By contrast, synovia from patients in PtC2, who tend to have less active and erosive disease, have a low inflammatory profile and intriguingly, evidence an M2-like macrophage signature.22 The M2 marker CD206 is expressed by tissue-resident macrophages, thought to play a critical role in maintaining synovial tissue homeostasis.25 26 The presence of CD206+ macrophages in the synovial lining has indeed been shown to be inversely correlated with clinical parameters of disease activity and histological parameters of synovial hyperplasia and inflammation in long-standing RA.27 These synovial tissue macrophages (STMs) are phenotypically, transcriptionally and metabolically reminiscent of IL4/IL13 polarised, M2-like macrophages, described in vitro.16 Here, we show, in a large series of samples, that (CD206- M1-like vs CD206+ M2-like) macrophage composition of the synovial lining is associated with synovial inflammation even in early, untreated RA. Although further mechanistic studies will be required to separate cause from consequence, it is possible that CD206+ STMs (or particular subsets thereof27) play an important protective role in dampening synovial inflammation in RA. In this regard, their assessment in unaffected joints of patients with RA, and longitudinally in follow-up biopsies, could be informative. An open question is whether CD206+ STMs in healthy synovia are replaced by CD206- macrophages in RA, or whether they are reprogrammed in the inflammatory context—and if so, what factors may allow for more efficient (or lasting) reversal of this plasticity.

Overall, our data suggest that the synovial immune profile in RA reflects a continuum from high to low inflammatory burden and a corresponding gradient of disease activity; whether one presages the other remains to be determined. It is also possible they reflect different stages of the same, dynamic pathogenic process. While an important strength of our study is the inclusion of only untreated patients with recent onset of RA symptoms (<12 months), it is known that synovial alterations may occur long before the development of clinically overt synovitis.28–30 The proportion of knee joints was significantly higher in the PtC2 than in the PtC1 group (table 1), yet it is unlikely that the groups cluster simply based on joint type, as we and others have observed similar patterns of synovial inflammation and gene expression when comparing large versus small joints from the same patients.31 32 Moreover, previous studies have not found an enrichment of knee joints in low inflammatory synovitis as assessed by histology.2 We applied variance filtering on gene expression exclusively to the (n=55) synovia from knees, and obtained very similar results in terms of both overlap between gene-lists (and therefore, pathway enrichment), and expression-based unsupervised clustering of samples (data not shown).

We acknowledge several limitations in this study. First, this is an observational cohort in which steroid treatment was not controlled. Longer term follow-up and repeat biopsies, especially in poor responders to methotrexate therapy, would in addition be of great value in determining the clinical relevance of our findings. Consequently, confirmation in a large, controlled, prospective setting is required. While (a minimum of 3) fragments sampled from different locations within joints were used for transcriptional profiling, our methodology may not fully account for intra-joint heterogeneity.33 The robust correlations observed between RNAseq signatures and histological data (from additional fragments from the same joints) nevertheless support the representativeness of our sampling. In particular, the association of M2-like/tissue resident macrophages with lower disease activity and synovial inflammation is striking and provides an avenue for further mechanistic insights. Spatial single cell transcriptomic profiling will be a powerful starting point, particularly in examining the subphenotypes of STMs,27 their proximal synoviocytes and the interaction between these compartments.