Dynamic spectrum of ectopic lymphoid B cell activation and hypermutation in the RA synovium characterized by NR4A nuclear receptor expression.

Ectopic lymphoid structures (ELS) can develop in rheumatoid arthritis (RA) synovial tissue, but the precise pathways of B cell activation and selection are not well understood. Here, we identify a synovial B cell population characterized by co-expression of a family of orphan nuclear receptors (NR4A1-3), which is highly enriched in RA synovial tissue. A transcriptomic profile of NR4A synovial B cells significantly overlaps with germinal center light zone B cells and an accrual of somatic hypermutation that correlates with loss of naive B cell state. NR4A B cells co-express lymphotoxins α and β and IL-6, supporting functions in ELS promotion. Expanded and shared clones between synovial NR4A B cells and plasma cells and the rapid upregulation with BCR stimulation point to in situ differentiation. Together, we identify a dynamic progression of B cell activation in RA synovial ELS, with NR4A transcription factors having an important role in local adaptive immune responses.

INTRODUCTION

Rheumatoid arthritis (RA) is a chronic autoimmune disease characterized by inflammation of the synovial tissue leading to joint damage and disability. Although dramatic advances in treatment options for RA have been achieved over the last two decades, a significant number of RA patients do not achieve remission or low disease activity (Drosos et al., 2019; Kaltsonoudis et al., 2019), highlighting the need for new therapies and biomarkers of response. B cells play a key role in RA disease pathogenesis both through autoantibody-mediated and antibody-independent functions and are key treatment targets in the disease (Derksen et al., 2017; Harre et al., 2012; Meednu et al., 2016; Sun et al., 2018). Although B cell clonal expansion is detected in RA blood (Lu et al., 2018) and synovium (Scheel et al., 2011), the mechanisms and location of antigen-specific B cell priming remains unclear.

B cell aggregates are present in the RA synovium (Humby et al., 2019; Lewis et al., 2019) at both early and late stages of disease and have been linked to important clinical outcomes including treatment response, bone erosion, and radiographic progression (Humby et al., 2019; Meednu et al., 2016). Ectopic lymphoid structures (ELS) are found in at least 40% of RA patient synovia with 10%–25% of those structures displaying features of a functional germinal center (GC) (Manzo et al., 2010; Rooney et al., 1988; Takemura et al., 2001; Yanni et al., 1993). ELS organization can recapitulate the architecture of secondary lymphoid organs (SLOs), with areas that resemble dark zones (DZs) (proliferating and AID-expressing B cells) and areas that resemble light zones (LZs) (rich in follicular dendritic cells and T cells in addition to B cells) (Nerviani and Pitzalis, 2018). ELS in RA synovial tissue has been shown to support antigen-driven selection and differentiation of autoreactive B cells, facilitating diversification, somatic hypermutation (SHM) (Manzo et al., 2010), and class switching (Humby et al., 2009) for in situ production of autoantibodies (Corsiero et al., 2016), suggesting that these structures are functional. Many critical factors driving ectopic lymphoid neogenesis, including lymphoid chemokines, are shared with factors driving SLO development, but the cellular components that produce these key factors may differ and are the subject of ongoing investigation. In mouse colitis, lymphotoxin (LT)-expressing B cells can support ELS formation and lead to severe inflammatory disease (Lochner et al., 2011). The role of B cells during ELS in RA synovial tissue, however, is unclear and which B cell subsets may be crucial for this process is still unknown.

Single-cell analyses of RA synovial tissue recently identified four B cell states including naive, memory, age-associated B cells, and plasma cells (PCs) (Zhang et al., 2019). To improve our understanding of the cellular, transcriptional, and antibody repertoire dynamics during human B cell activation in ELS, here we performed unbiased single-cell transcriptomic and repertoire profiling of RA synovial B cells. These single-cell antibody repertoires paired with single-cell transcriptomics allowed us to define transcriptional B cell states, resolving functional B cell heterogeneity and gene expression dynamics. We identified an abundant population of synovial B cells characterized by upregulation of NR4A1-3 along with other immediate-early response genes (EGR1/3, FOS, and JUN) and activation markers (CD69 and CD83). These cells show evidence of SHM, class switching, and clonal expansion. Gene expression profiling of synovial NR4A B cells demonstrated enrichment for genes associated with GC LZ B cells and expression of elevated levels of chemokine receptors and chemotactic factors involved in ELS, including LT-α, LT-β, and IL-6. This subset of B cells is enriched in synovial tissue but found at very low levels in peripheral blood and another autoimmune target tissue, the kidney in lupus nephritis. Notably, RA synovial B cells and PCs spontaneously express NR4A1 protein by flow cytometric and histologic analysis, and NR4A cluster genes are enriched in synovial samples with a lymphoid histologic pathotype (Lewis et al., 2019). In an RA flare cohort (Orange et al., 2020), an NR4A B cell transcriptomic signature peaked in the blood during flare. The presence of expanded and shared clones in this activated B cell state and NR4A-expressing synovial PCs and the rapid upregulation with BCR stimulation points to antigen-driven activation and differentiation in the synovium. Together our data support in situ activation of B cells with accrual of SHM in synovium ELS and NR4A expression as a potent readout of antigen stimulation, local adaptive immunity, and pathological B cell responses.

RESULTS

scRNA-seq analysis of RA B cells reveals multiple distinct subsets, including a NR4A B cell population in the synovium

To define the diversity of B cells in the synovium, we isolated B cells from lymphoid-rich RA synovial tissues and paired blood and performed single-cell transcriptomic and BCR sequencing on the 10× Chromium platform (Figure 1; Table S1). The mean frequency of B cells by flow cytometry was 12.65% (2.7%–21.6%, n = 4) of CD45+ synovial cells. The majority of synovial B cells were IgD−, whereas 80% of blood B cells were naive (IgD+CD27−). We also observed synovial enrichment of CD24−CD27hi plasmablasts/PCs (1.9%–7.6%) (Figure S1B).

Figure 1.

Overview of work flow

Synovial tissues were disaggregated using a validated protocol (Donlin et al., 2018) and B cells were sorted by flow cytometry cell sorting. We then performed scRNA-seq on sorted tissue and blood B cells using the 10× genomic platform with poly-A-selected, 5′ initiated expression and V(D)J libraries generated from each single cell. See also Table S1 and Figure S1.

We sequenced 4,115 B cells (2,589 synovial and 1,526 peripheral blood) with an average median genes per cell of 1,427 (755–1,868). After exclusion of cells with high mitochondrial RNA content, dimensionality reduction using t-distributed stochastic neighbor embedding (t-SNE) (van der Maaten and Hinton, 2008) and unsupervised clustering using Seurat package (Butler et al., 2018) on the combined dataset of 3,786 cells defined 8 conserved B cell subsets or clusters. Two subpopulations of PCs (PC(i) and PC(ii)) were characterized by high expression of transcription factors required for PC differentiation and maintenance (XPB1, IRF4, and PRDM1) and immunoglobulin genes (Figures 2A, 2B, and 2D). Expression of IGHD and TCL1A identified three naive B cell subpopulations (naive(i), naive(ii), and naive(iii)) (Figures 2A and 2B). Two B cell subpopulations expressed high levels of CD27, HOPX, S100A10, and S100A4, markers associated with memory B cells (Bautista et al., 2020; Descatoire et al., 2014; Klein et al., 1998; Riedel et al., 2020). The LMNA+ designated subpopulation also uniquely expressed the LMNA gene, not previously reported in memory B cells (Figures 2A, 2B, and 2D). The final subpopulation was labeled NR4A+ as it expressed high levels of the nuclear family receptor 4 (NR4A1, NR4A2, and NR4A3). This subpopulation also showed high expression of markers associated with activation and GC formation, such as CD69, CD83, and GPR183, the latter a molecule mediating B cell migration in lymphoid follicles (Gatto et al., 2011; Pereira et al., 2009) (Figures 2A, 2B, and 2D). NR4A was also expressed by the LMNA cluster and PC(i) (Figure 2B).

Figure 2.

Single-cell RNA sequencing identifies an enrichment of B cells expressing nuclear orphan family receptors (NR4As) in RA synovial tissue

(A) t-SNE visualization of 8 B cell clusters from 3,786 B cells from 4 RA synovial tissues (SYN) and one paired blood sample (BLD).

(B) Markers identified two clusters of plasma cells, three clusters of naive, a memory, an LMNA+, and an NR4A+ cluster.

(C) The same t-SNE map as in (A) with cells labeled by sample ID. See also Figure S2C.

(D) Heatmap displays top differentially expressed genes (DEGs) in each cell cluster. Top 10 DEGs and NR4A2 and NR4A3 of the NR4A+ cluster are magnified.

(E) Boxplots display frequency of NR4A+ B cells calculated using SingleR. Single-cell RNA sequencing of synovium (AMP) (n = 14) from (Zhang et al., 2019) (ImmPort SDY998), SLE (AMP) (n = 13) from (Arazi et al., 2019) (ImmPort SDY997), peripheral B cells from (Zheng et al., 2017), and RA134_BLD from this study (synovium, 10×). The horizontal blue lines are the mean and the blue boxes represent 90% CI.

(F) Bar graphs plotting log2 fold change of NR4A+ cluster genes in leukocyte-poor RA, leukocyte-rich RA, or RA biopsy versus OA using bulk RNA sequencing from (Zhang et al., 2019) (ImmPort SDY1299). *Significant at 5% and **significant at 1%.

The NR4A+ and LMNA+ subpopulations were restricted to synovial samples. Naive B cells were the main population in the blood, whereas the memory and PC subpopulations were present in both blood and synovial tissue (Figures 2C and S2C).

Unique enrichment of NR4A+ B cells in RA synovium compared with blood and other tissues

Using a supervised classification method SingleR (Aran et al., 2019), we next examined the presence of NR4A+ B cells in an independent RA synovial tissue cohort and other tissues where single-cell data were available. In SLE kidney (Arazi et al., 2019) and peripheral blood B cells (RA from our data and a published healthy donor [Zheng et al., 2017]), NR4A+ B cells were very rare at 0.7%–1.5% abundance. In contrast NR4A+ B cells were highly abundant (>40% abundance) in synovial tissue from AMP phase I (Zhang et al., 2019), similar to samples in our current study (Figure 2E). In the AMP phase I bulk RNA sequencing cohort (Zhang et al., 2019), we found that several markers, including CD83 and NR4A1, which defined the NR4A+ cluster, were upregulated in B cells from RA biopsy samples compared with B cells from RA and OA arthroplasty samples (Figure 2F). This NR4A cluster enrichment appeared to be driven by the inflammatory state of the tissue (1.7- and 2.1-fold increase in CD83 and NR4A1 gene expression leukocyte-rich RA versus OA, p < 0.05; NS for leukocyte-poor RA versus OA) (Zhang et al., 2019).

NR4A+ synovial B cells accrue SHM

To define the Ig repertoire in single B cells, we performed V(D)J sequencing on the same sorted CD19+ cells used for single-cell RNA sequencing (scRNA-seq), allowing a joint analysis of transcriptional profile and immune repertoire in the same cell (Figures 1 and S2). Similar to previously published data, all three naive clusters exhibited very low SHM frequencies, while memory and PC clusters showed the highest mutation rates (King et al., 2021; Klein et al., 1997; Tangye and Good, 2007) (Figures 3A and S2F). The LMNA+ cluster has a mutation rate comparable with that of memory B cells, supporting that LMNA+ are memory B cells that may have developed in the synovium. It is curious that BCR recovery in the LMNA+ cluster is low overall, but it is not clear if this is a technical or biologic issue. Although mutated cells are apparent in the NR4A+ cluster, the average SHM rate was not significantly higher than naive B cells (1.5%, 95% CI, 1%–4%) while memory and PCs had clear differences compared with naive, suggesting that NR4A+ cells may be in an earlier stage of the affinity selection and maturation process (Figure 3A). Furthermore, we detected IGHA expression in a small number of NR4A+ B cells, providing evidence of class-switch recombination (Figure S2E).

Figure 3.

B cells in NR4A+ clusters display somatic hypermutation (SHM) and clonal expansion

(A) Violin plots display SHM rate, averaged across detected chains, in each B cell cluster obtained from single-cell BCR sequencing. Each cluster’s SHM rate was tested against naive(i) using linear mixed models, with p < 0.05 displayed. See also Figures S2E and S2F.

(B) Scatterplots display the associations between percent SHM and expression level of CD27, NR4A1, NR4A2, and IGHD in NR4A+, memory and plasma cell (combined PC(i) and PC(ii)) clusters. Blue lines display a linear mixed model fit. p values are from a linear mixed effect model that adjusts for sample and number of genes detected. Only values of p < 0.05 are displayed.

(C) Number of cells with >97% DNA identity of the heavy-chain CDR3. The amino acid sequence of each putative clone is shown on the left.

(D) Heavy-chain BCR data were superimposed on cluster t-SNE plot. “None” (gray) indicates cells from which BCR was not recovered, “Recovered” (black) indicates cells from which BCR was recovered, “Expanded” (blue) shows cells with clonality in ≥2 cells, “CARHWRGKKPFDSW” (red) indicates cells with a shared clone recovered from both NR4A+ and plasma cells. See also Figure S2D.

Next, we looked at the relationship between SHM rate and the expression levels of CD27, NR4A genes, and IGHD. As expected, all subpopulations showed increased SHM with increased expression of CD27, while NR4A and memory subpopulations showed inverse associations of SHM with IGHD. However, particular to only the NR4A+ cluster, an accumulation of SHM positively associated with the expression levels of NR4A1 and NR4A2 (Figure 3B), suggesting a potential role for NR4A genes in the RA synovial immune reaction. There was no significant association found between the expression of either the NR4A genes or IGHD in naive and LMNA+ clusters. Importantly, we detected expansion of 43 families of putative clones in the synovium and five in the peripheral blood, most prominently featured in the PC and NR4A+ clusters. Clonal families were almost entirely confined to a single transcriptomic subpopulation (naive, memory, LMNA, NR4A, or PC), a finding not expected to occur by random chance (p < 0.001 by permutation test). The only exception was the clonal family CARHWRGKKPFDSW, which was detected in both the NR4A+ and PC cluster, consistent with a developmental relationship between these subsets (Figures 3C, 3D, and S2D).

NR4A+ synovial B cells are activated and resemble GC B cells

To further clarify the eight synovial B cell clusters, we performed gene set enrichment analysis. All three naive clusters were significantly enriched with naive B cells from blood and tonsil (Bm1+Bm2) (Figure 4A). In contrast, memory and LMNA+ clusters were enriched for gene sets expressed by tonsil Bm5 (a memory subset) (Bohnhorst et al., 2001) and, to a lesser extent, memory B cells from blood, supporting their identity as memory B cells (Figure 4A). The analysis also confirmed the identity of PC clusters (Li et al., 2014) (Figure 4A). The NR4A+ cluster was strongly enriched for genes expressed by GC LZ B cells (Figures 4A and S3B), including EGR1/2/3, CD83, BCL2A1, and GPR183 (Victora et al., 2012) (Figure S3A).

Figure 4.

NR4A+ cluster displays a spectrum of activation from naive state to an activated germinal center light zone (GC LZ) phenotype

(A) Heatmap displays signed-log10 (FDR q value) of gene set enrichment analysis (GSEA) between our B cell clusters and published blood and tonsil B cell subsets. See also Figure S3B. Principal component analysis projecting gene expression data of naive(i), naive(ii), naive(iii), NR4A+, LMNA+ and memory B cell subsets onto two-dimensional space. See also Figures S3E, S4A and S4B.

(C) Heatmap displays expression of the top 40 loading genes for dimension 1 with increasing score from naive(i), naive(ii), naive(iii), and N4A+ subsets.

(D) Plots show example of genes increase (NR4A1 and DUSP1) or decrease (TXNIP and CD79B) as cells transition from naive to NR4A+ state. Colors on the heatmap from blue to red represent normalized expression from low to high. The black graph color shows gene loading score along dimension 1 (PC1).

Recently, scRNA-seq of GC B cells from conventional SLOs has revealed a dynamic spectrum of DZ, LZ, and intermediate zone B cells (Holmes et al., 2020; Milpied et al., 2018). In particular, between the DZs and LZs there was a gradient of downregulation of classic DZ genes, such as CXCR4, and upregulation of classic LZ genes, such as CD83. Along this gradient, there was an initial upregulation of BCR signaling followed by NFκB signaling. The later stages of the intermediate zone GC B cells were characterized by upregulation of CD40 and MYC signaling, indicative of their potential to interact with T cells (Holmes et al., 2020). In our data, the NR4A+ B cells in RA synovial tissue displayed characteristics compatible with these intermediate cell stages. First, in contrast to classic LZs, NR4A+ cluster cells expressed CXCR4, BCL6, and FOXP1, genes associated with GC DZs (Caron et al., 2009; Holmes et al., 2020; Victora et al., 2012) (Figures S3A–S3C). Second, NR4A+ B cells showed evidence of BCR stimulation (Figure S3D) and high levels of immediate-early genes associated with recent BCR stimulation, including NR4A1, NR4A3, EGR1, and c-Fos (Monroe, 1988; Ollila and Vihinen, 2003; Ye et al., 2010) and NFκB signaling (Figure S3C). In contrast to classic LZ cells, the NR4A+ cluster does not show consistent upregulation of CD40 and MYC signaling, with the exception of BCL2A1 (CD40 signaling) and GPR183 (MYC signaling) (Figure S3C). Our NR4A+ B cell cluster shares remarkable transcriptomic similarity with a recently described activated B cell state in the human tonsil that appears to be on a trajectory to GC formation (King et al., 2021).

A gene expression activation continuum from naive to NR4A B cells

The variability in NR4A1, NR4A2, IGHD, and CD27 expression in the NR4A subset (Figures 2B and 2D) led us to hypothesize that B cells in synovium might exist along a continuum from a naive state. Projecting the gene expression data of the B cells onto two dimensions using principal component analysis, we sought to further define the relationship between B cell subsets. PCs were excluded from this analysis. Dimension 1 (4% variance) separated NR4A+ cells from naive. In contrast, dimension 2 explained the difference between the NR4A+ and memory subsets (2% variance) (Figure 4B). LMNA+ cells were positive in both dimension 1 and 2. In dimension 1, we found that genes with large gene loadings are the ones involved in lymphocyte activation and GC molecular events, including CDKN1A, NFKBID, GPR183, CD83, LY9, KLF6, and IGHG1, and are expressed in increasing levels from naive to NR4A+ and LMNA+ clusters. Consistently, the NR4A+ cluster is significantly enriched for gene sets associated with cell activation (Figure S3E). In addition, a number of DNA-binding transcription factors (NR4A2, NFKBID, NR4A3, JUND, DNAJA1, CREM, and FOS) gradually increase from naive to NR4A state, while genes characteristic of naive B cells (TCL1A, IGHM, and FCER2) are downregulated (Figure 4C). Along dimension 1, we observed cells in transition states between naive and NR4A+. Thus, expression of NR4A+ markers, including NR4A1 and DUSP1, are gradually increased, while, concurrently, markers associated with naive status, such as TXNIP and CD79B, showed a gradual decrease in expression (Figure 4D). In dimension 2, genes with increased expression levels are those associated with memory B cells, including HOPX, IGHA1, S100A4, and CD99. We also observed evidence of transitioning between NR4A+ and memory state with increased expression of memory genes (S100A6 and HOPX) and decreased expression of NR4A+ genes (JUN and DUSP1) (Figures S4A and S4B). Our data are consistent with a continuum of states from naive to NR4A+ with both loss of naive and acquisition of activation status. Whether LMNA+ cells differentiate from memory, or NR4A+, or both remains to be defined.

Highly expressed chemokines and cytokines in NR4A B cells

Chemokines and their interaction with corresponding chemokine receptors are vital for the organization of lymphoid organs and B cell functions (Nerviani and Pitzalis, 2018). A number of ELS promoting factors are upregulated in lymphoid-like RA synovial tissue (Manzo et al., 2007, 2008). To investigate the relevant chemokine and cytokine pathways in synovial tissue, we next examined differential gene expression of chemokines, cytokines, and their receptors in each B cell cluster. Both naive and NR4A+ B cells showed significant enrichment for CXCR5, CCR7, and CXCR4 (Figures 5A and 5B). These chemokine receptors are instrumental to orchestrate homing to B cell follicles (CXCR5), migration across high endothelial venules (CCR7), and dynamic mobilization inside GC (CXCR4) (Miyasaka and Tanaka, 2004; Okada et al., 2002), suggesting that these B cells can be recruited to the synovium and localize in response to specific chemokine signals. Receptors for B cell survival factors are also expressed by multiple B cell subsets. Naive and NR4A+ subsets showed significant upregulation of TNFRSF13C (BAFFR), while LMNA+ and memory shared upregulation of TNFRSF13B (TACI), a receptor that can bind to both B cell survival factors BAFF and APRIL. Consistent with previous reports, PC uniquely showed upregulation of TNFRSF17 (BCMA), a BAFF receptor important for PC survival (Figure 5A). Furthermore, NR4A+, LMNA+, and memory showed significant enrichment of TNFRSF1 (TNFR2), a receptor for TNF (Figure 5A). We also observed that LMNA+ and memory subsets uniquely upregulate SIGIRR, a negative regulator of IL-1 receptor signaling pathway, implying a response to IL1 family cytokines (Figure 5A).

Figure 5.

NR4A+ cluster expresses chemokine receptors and cytokines important for ELS formation

(A) Heatmap displays log2 fold change of cytokine and chemokine receptors in each B cell subset (n = 1,250 naive, 746 NR4A+, 633 LMNA+, 330 memory, and 779 PC). Volcano plot displays log2 fold change and −log10 p value of genes in the heatmap.

(B) Percentage of cells of each subset that express CXCR5, CCR7, CXCR4, or CCR6.

(C) Heatmap displays log2 fold change of cytokines and chemokines in each B cell subset. Volcano plot displays log2 fold change and −log10 p value of genes in the heatmap in NR4A+ and LMNA+ subsets. For (A) and (C), shown are genes that are expressed in at least 20% of the cells of at least one subset. p values are adjusted using the Benjamini-Hochberg correction for multiple tests.

(D) Percentage of cells of each subset that express IL6, TNFSF9, CD70, TGFB1, CCL5, or LTB.

(E) GSEA analysis revealed significant enrichment of PreM signature in NR4A+ and LMNA+ subsets, *p < 0.05. Heatmaps show log2 fold change of PreM genes (Holmes et al., 2020) in CCR6+ and CCR6− cells within NR4A+, LMNA+, or naive subsets. Violin plots show log2 expression of CCR6 in each cluster. Bottom right bar graph shows percentage of cells of each subset that co-express CCR6, CELF2, and BANK1. See also Figure S5.

B cells can play a direct role in the organization of SLOs by producing lymphotoxin β (Tumanov et al., 2002). More recently, B cell-derived IL-6 has been implicated in spontaneous GC formation in systemic lupus (Arkatkar et al., 2017). We found that all non-PC B cell subsets express lymphotoxin β (LTB) (Figures 5C and 5D). However, the NR4A+ subset uniquely expressed IL6 and co-expressed LTB and LTA (Figure 5C). As previously noted, GPR183 (EBI2) is highly expressed in the NR4A+ subset. GPR183 has pleiotropic roles in B cell positioning depending on the balance of other chemokines. GPR183, together with CXCR4, also expressed at high levels by the NR4A+ subset, are critical factors for the segregation of GC B cells in the DZs and LZs (Figure 5A). NR4A+ B cells also have elevated expression of a number of co-stimulatory and signaling molecules important for putative cell-cell interactions with APCs and T cells. In addition to IL6 (interaction with IL-6R on APCs and Tfh), this includes high levels of TNFSF9 (4-1BBL) (Figures 5C and 5D), a ligand that binds to 4-1BB on activated CD4 and CD8 T cells to provide co-stimulatory signals (Cannons et al., 2001; Vinay and Kwon, 2011). Overall, these data suggest that NR4A+ B cells may be a critical initiator of ELS formation in synovial tissue. Independent analysis of the RA cohort from AMP also supports these findings (Figure S5). Of note, we observed the presence of histologic ELS in our samples by immunofluorescent staining with classic GC markers (Figure S6).

One of the differentially expressed chemokine receptors in the NR4A+ subset is CCR6 (Figures 5A and 5B). Single-cell analysis of GC B cells from human tonsil has identified a distinct cell cluster of memory B cell precursors (PreM) originating within the LZ that shares a gene signature with previously reported murine CCR6+ PreM B cells (Holmes et al., 2020; Suan et al., 2017). The range of CCR6 expression within the NR4A+ cluster (Figure 5E) suggests that memory precursor cells are originating here and is consistent with memory precursor commitment to post-GC differentiation in the early stage of the LZ (Holmes et al., 2020). Consistent with CCR6 expression, we found that NR4A+ and LMNA+ cells are significantly enriched in a PreM signature (Figure 5E). In contrast to naive subsets, CCR6+ cells in the NR4A+ subset showed higher expression of most of the genes described as a signature of PreM, including CELF2, BANK1, and CD44, compared with CCR6− counterparts, while LMNA+ CCR6+ displayed a smaller subset of PreM signature genes (Figure 5E). In addition, the NR4A+ subset contained high fractions of cells co-expressing CCR6, CELF2, and BANK1 (Figure 5E).

Prominent expression of NR4A1 protein in synovial B cells and PCs in ELS

To further define the phenotype of synovial NR4A1 B cells and their localization in ELS, we next performed immunofluorescent staining on synovial tissue sections. RA synovial tissue histology identified clear GC structures, with CD21 staining for follicular dendritic cells (FDCs) in the LZ, proliferating nuclear antigen for B cells in the DZ, and IgD for the naive follicular mantle. NR4A1 expression was observed in B cells within and in close proximity to the GC and prominently in aggregates of CD138+ PCs (Figure 6A; see Figure S6 for additional representative images). Interestingly, NR4A1 was also expressed by select T cells localized within and near the GC (Figures 6A and S6). As endogenous NR4A1 (Nur77) expression is an indicator of antigen receptor signaling in B and T cells (Ashouri and Weiss, 2017; Mueller et al., 2015; Zikherman et al., 2012), synovial lymphocyte NR4A1 expression suggests in situ antigen activation, and persistent expression in synovial PCs suggests antigen selection and differentiation in the tissue. The majority of the PCs expressing NR4A1 are CD20+, indicative of a recently generated plasma cell or plasmablast (Figures 6A and S6, yellow overlay cells in CD20-NR4A1-CD138 panels). Histology of tonsil tissue confirmed the expression of NR4A1 in B cells in close proximity to CD21+ LZ FDC inside GC (Figure S7A) and in CD138+ PCs surrounding GCs (Figure S7B).

Figure 6.

NR4A1 is expressed by B cells and plasma cells in synovium ELS

(A) Synovial serial sections were stained with four antibody combinations, as indicated on the top of each image. The dashed yellow line outlines a GC containing CD20+IgD−PCNA+ centroblasts. Green and white boxes mark the area magnified in the right images to show details of cellular morphology and subcellular location of NR4A1. Corresponding single-color images for CD20, CD3, CD138, and NR4A1 are shown on the right. # and * indicate location of cropped areas. Green arrows point to NR4A1+CD20+ B cells, yellow arrows depict NR4A1+CD3+ T cells, and white arrows show areas with significant accumulation of NR4A1+CD20+CD138+ plasmablasts on the periphery of secondary B cell follicles in the synovium. Scale bars, 100 μm in (100× magnification pictures).

(B) Bar graph represents mean of percent NR4A1+ B cells in ex vivo sample from NC blood (n = 5), RA blood (n = 5), and RA synovial fluid (n = 6).

(C) Example dot plots showing the distribution of NR4A1+ B cells from synovial fluid. See also Figure S8A.

(D) Bar graphs display distribution (%) of NR4A1+ or NR4A1− B cells across B cells subsets: naive, double-negative (DN), unswitched memory (USM), switched memory (SM) (this gate may include PC), and plasma cell (PC) from synovial fluid (n = 7). PC is defined as CD19+IgD−CD27hiCD38hi.

(E) Bar graphs display percent of naive, DN, USM, SM, and PC that are NR4A1+.

(F) Histogram of NR4A1 expression by naive, USM, SM, DN, and PC, from synovial fluid (SynF) compared with blood naive and NR4A1 fluorescent-minus-one control (FMO). Right line graphs show NR4A1 geometric mean fluorescent intensity (gMFI) in naive, USM, SM, DN, and PC from synovial fluid (n = 6). Each line represents an individual sample.

(G) Dot plots display NR4A1+ and NR4A1− B cell distribution across B cell subsets from synovial tissue. Histogram represents NR4A1 expression by B cell subsets gated using CD38 and IgD. Line plot shows NR4A1 gMFI in different B cell subsets from synovial tissue (n = 1). In (B), (D), and (E), error bars are SEM. Statistical test was non-parametric one-way ANOVA and Tukey’s multiple comparisons test. *p < 0.05, **p < 0.01.

(H) NR4A1 expression in B cells from PBMC treated with anti-Ig(A + M + G) for 4 h. Scatterplots represent percent NR4A1+ B cells after treatment (S) compared with untreated control (U). **p < 0.01 by paired t test. See also Figures S8B and S8C.

Flow cytometry of B cells from disaggregated RA tissue and synovial fluid revealed that a significant fraction spontaneously express NR4A1. In marked contrast, only a small percentage of peripheral blood B cells (from healthy controls or RA) express NR4A1 (Figure 6B). NR4A1+ B cells from synovial fluid are predominantly IgD−, including IgD−CD27+ switched memory (SM), double-negative (DN) B cells, and PCs (Figures 6C and 6D). However, NR4A1 expression is seen in a fraction of all synovial B cells, regardless of phenotype. Memory B cells and PCs express significantly more NR4A1 than naive cells (Figures 6E and 6F). In synovial tissue, NR4A1 B cells are also predominantly IgD− (Figure 6G). Here, we characterized the B cells based on CD38 and IgD expression, as used to classify tonsil B cells (Bohnhorst et al., 2001). NR4A1 B cells include CD38+IgD− (Bm3/4 compartment containing GC cells and PCs) and CD38−IgD− (Bm5 compartment containing CD27+ memory and DN cells) (Figures 6G and S8A).

To establish the pattern of NR4A1 expression during normal B cell activation in SLOs, we performed additional flow cytometry analysis of tonsil B cells. Similar to the synovium, multiple subsets (Bm1+Bm2, Bm5, Bm3/4: GC and PC) of tonsil B cells express NR4A1, but with significantly higher expression in PC and GC B cells (Figure S7C). Furthermore, when we subset GC B cells based on CD83 and CXCR4, significantly more NR4A1+ cells were present in CD83+CXCR4+ intermediate GC B cells than in CD83−CXCR4+ (DZ), CD83+CXCR4− (LZ), or CD83− CXCR4− GC B cells, consistent with our RNA sequencing data that synovial NR4A+ B cells express high levels of both CXCR4 and CD83 (Figure S7D).

When we stimulated peripheral blood and synovial B cells with a cocktail of anti-IgG, anti-IgA, and anti-IgM to engage BCR of multiple isotypes, NR4A1 expression was significantly upregulated (Figure 6H). The NR4A1 mean fluorescent intensity (MFI) was significantly increased in stimulated B cells compared with untreated cells (Figure 6H). Notably, in synovial B cells, the NR4A1 MFI in untreated cells is already higher than in untreated blood B cells and further increased after stimulation (Figure 6H). In addition to NR4A1 upregulation, we also found that BCR stimulation upregulated NR4A2 and NR4A3 expression at the mRNA level (Figures S8B and S8C).

Lymphoid-specific NR4A subset genes correlate with RA synovial tissue pathotype and increase in the blood during RA flare

Next, we sought to relate NR4A subset signature genes to synovial histology and clinical characteristics in an independent cohort of synovial biopsies from treatment-naive patients: the Pathobiology of Early Arthritis Cohort (PEAC) (Lewis et al., 2019). Based on histology scores, synovial samples were classified as lympho-myeloid (CD20 B cell aggregate rich), diffuse-myeloid (CD68 rich in the lining or sub-lining layer but poor in B cells), or fibroid (paucity of immune-inflammatory cell infiltrate), as described previously (Lewis et al., 2019). We focused on the expression of CD83, GPR183, and Ly9 as more lymphoid-specific genes representative of the NR4A cluster as only bulk RNA sequencing of total tissue is available in this study. These genes were strongly enriched in the lymphoid pathotype tissues (Figure 7A) and also showed significant correlation with synovial histology scores for CD20, CD3, CD138, and sub-lining CD68 (Figure 7B). These data suggest that infiltration of multiple immune cell types associated with ectopic lymphoid responses in the synovial tissue may be linked to the B cell NR4A subset. In addition, CD83 and LY9 expression showed a relationship with ultrasonographic joint synovial thickness confirming that ELS gene expression strongly matches imaging signs of active joint inflammation in the particular joint undergoing biopsy (Figure 7B). Surprisingly, NR4A1, NR4A2, and NR4A3 whole synovial gene expression was lower in the lymphoid versus fibroid pathotype and inversely correlated with histologic immune cell infiltration (Figure 7B) possibly due to B cell independent roles for the NR4A family in the synovium in other cell populations that may be anti-inflammatory (Alivernini et al., 2020).

Figure 7.

Lymphoid-specific NR4A subset genes correlate with RA synovial tissue pathotype and increase in blood during RA flare

(A) Boxplots of histology by tertile demonstrating correlation with CD83 (Padj lymphoid versus myeloid = 3.0e–02; Padj myeloid versus fibroid = 3.4e–02; Padj lymphoid versus fibroid = 3.3e–07), GPR183 (Padj lymphoid versus myeloid = 5.4e–03; Padj myeloid versus fibroid = 4.6e–02; Padj lymphoid versus fibroid = 4.1e–08) and LY9 (Padj lymphoid versus myeloid = 1.4e–07; Padj myeloid versus fibroid = 8.0e–01; Padj lymphoid versus fibroid = 5.0e–07). p values were calculated by linear regression models.

(B) Correlation heatmap showing spearman correlations of CD3, CD20, CD138, CD68L, and CD68SL histology scores and ultrasound biopsy joint parameters (ST, synovial thickness; PD, power Doppler) at the biopsy joint (ultrasound ST/PD BJ) against CD83, GPR183, LY9, and NR4A1-3 expression. The stars represent the significance of the correlation coefficient: *p < 0.05, **p < 0.01, ***p < 0.001. The plots and p values were analyzed from Lewis et al. (2019).

(C) Relative expression of marker genes (see STAR Methods for details of calculation) from the synovial B cell subpopulations in the weeks before and after flare (Orange et al., 2020). See Table S2 for list of overlapping genes. Top panel shows the mean RAPID3 disease activity for the four donors studied. The dashed vertical line marks the time at flare. See also Figure S9.

Although our initial data suggested that NR4A B cells are rare in RA blood (Figure 6B), we wondered whether circulating numbers of this activated B cell population may vary over time depending on disease activity. To begin to understand this dynamic, we compared our synovial B cell cluster signatures with blood transcriptional profiles recently defined relative to flare (Orange et al., 2020). In this dataset, two transcriptomic clusters appeared in the blood preceding flare: AC2 was increased 2 weeks before flare and AC3 was increased 1 week before flare. AC2 was originally described as having characteristics of naive B cells, while AC3 is a unique CD45−CD31−PDPN+ mesenchymal population, termed PRIME (preinflammatory mesenchymal) cells (Orange et al., 2020). An additional cluster 1 was increased at the onset of flare. When we examined the overlap between these blood clusters and our single-cell B cell clusters (Figures 7C and S9A; Table S2) using hypergeometric tests for enrichment, we found that AC2 was significantly enriched with transcripts characteristic of peripheral blood dominant naive B cell clusters, although the synovial NR4A+ and LMNA+ clusters were represented (Figure S9A). In contrast, the flare cluster 1 was significantly enriched with naive(iii), NR4A+, and LMNA+ synovial populations (Figure S9B). This is well visualized in the line graph (Figure 7C), with resting naive (naive(i), naive(ii)) and memory B cells increasing pre-flare and the NR4A cluster genes increasing at the same time as the increase in disease activity (RAPID3). This is consistent with a trajectory of B cell activation whereby resting blood B cells are recruited from the blood to the synovium pre-flare and then activated within the synovial microenvironment to express NR4A during flare.

DISCUSSION

Our integrated single-cell analysis of RA synovial B cells has revealed a surprising level of B cell activation, with a large fraction of synovial B cells demonstrating evidence of recent antigen stimulation, as revealed by upregulation of the NR4A gene family. NR4A B cells appear primed for cell-cell interactions with APCs and T cells and, further, have evidence of molecular events associated with active GC reactions, including class-switch recombination, SHM, and clonal expansion with selection into the PC compartment. Our findings suggest a spectrum of activation of B cells within the synovial microenvironment with a pluripotent fate that includes differentiation to memory or PC lineage in situ.

Previous studies have supported local synovial activation of RA B cells with evidence of SHM and clonal expansion (Corsiero et al., 2016; Scheel et al., 2011) using a number of approaches, including enzymatic digestion or focal microdissection with bulk or single-cell BCR sequencing. An important advance in this study is the high-throughput combination of single-cell RNA sequencing and Ig gene repertoire. This allowed the unprecedented identification of a spectrum of B cell activation in synovial B cells downstream of BCR signaling as revealed by NR4A1 expression. A key finding from our study was the relationship between the degree of NR4A family gene expression and SHM. Notably, this dynamic process was unique to the synovium and not observed in B cells from RA peripheral blood or lupus kidney. NR4A1-3 encode a small family of orphan nuclear hormone receptors (NUR77, NURR1, and NOR1, respectively), which are rapidly induced by acute and chronic antigen stimulation in B and T cells (Ashouri and Weiss, 2017; Mittelstadt and DeFranco, 1993). Our results raise key questions regarding the precise synovial microenvironmental signals that promote NR4A family expression and second signals that may direct subsequent differentiation to memory B or PC. Another study did combine transcriptomic and repertoire analysis of peripheral blood RA B cells and interestingly suggested that ACPA-enriched B cell responses are imprinted with a T cell-dependent transcriptional network (Lu et al., 2018). Given the strong association of NR4A1 expression in B cells and BCR engagement (Zikherman et al., 2012), we speculate that the enrichment of B cell NR4A1 expression in the RA synovium represents chronic autoantigen stimulation and an ACPA response. This would be in keeping with other data in the literature indicating that B cells differentiating within ELS frequently target citrullinated peptides (Corsiero et al., 2016).

In terms of B cell fate, it is notable that NR4A family members may trigger B cell apoptosis after encounter with self-antigen in the absence of a second signal provided by T cell (or other) costimulation and act as a negative regulator of B cell activation (Tan et al., 2020). Thus, in murine B cells, the NR4A family, particularly NR4A1/NUR77, actively participates in B cell tolerance. NUR77 is upregulated in self-reactive B cells in response to chronic antigen stimulation and selectively restricts survival by limiting B cell survival factors such as BAFF (Tan et al., 2019). In addition, a small subset of GC LZ B cells exhibiting a high SHM rate upregulate NR4A1/NUR77, suggesting active BCR signaling and a role for NR4A1 in selection and regulation of GC differentiation pathways into the PC compartment (Mueller et al., 2015). Recently, NR4A1 and NR4A3 have been identified as partially redundant mediators of an antigen-dependent negative feedback loop in B cells to reinforce tolerance by increasing B cell dependence on T cell help and regulating clonal competition when T cell help is limited (Tan et al., 2020). These data point to an expanding role for the NR4A family in adaptive immune responses occurring in SLOs. However, the role of NR4A in ELS and autoimmunity had not been studied previously. Whether NR4A attempts to restrain B cell activation within the synovium, similar to its role in a recently reported anti-inflammatory synovial tissue-resident macrophage population (Alivernini et al., 2020) and exhausted T cells in the tumor microenvironment (Chen et al., 2019; Liu et al., 2019), remains to be defined. It is also not clear whether the synovial B cell activation reflected by NR4A expression is occurring by follicular and/or extra-follicular pathways (Jenks et al., 2018), although the high CXCR5 expression by NR4A B cells is at least consistent with follicular localization.

Our finding of shared clonality between the NR4A B cell and PC clusters and prominent expression of NR4A in synovial PCs suggests that NR4A B cells do indeed receive the signals required to drive their differentiation to memory and PC in situ. Similar findings have been reported in a mouse model of autoimmune arthritis where NR4A expression identified auto-antigen-specific CD4 T cells with Th17 differentiation potential (Ashouri et al., 2019). The signals that control differentiation of activated B cells into memory or PCs in ELS require further study, but likely depend on antigen receptor signaling and the strength and stability of interaction with T cells (Ise et al., 2018; King et al., 2021; Suan et al., 2017). We identified several genes within the NR4A cluster that may facilitate T-B cell interactions, including CD69 (which prevents activated lymphocyte egress from SLOs), CD86 (a ligand for CD28 on T cells), IL-6 (which can promote both APC and T cell interactions), and ICAM1 (which can support B cell interaction with T follicular helper cells) (Zaretsky et al., 2017). This is also in accord with a recent report showing that tetramer-identified ACPA-positive B cells in RA blood and synovial fluid express T cell-stimulating ligands and a persistently activated phenotype, suggestive of continuous antigenic triggering (Kristyanto et al., 2020). A key role for B cells in RA flare was also suggested by the recent discovery of a circulating expanded B cell subset (AC2) before the emergence of a circulating CD45−/CD31−/PDPN+, PRIME cell in RA patient blood before clinical flare (Orange et al., 2020). We find that this pre-PRIME circulating B cell cluster (AC2) was enriched with resting naive (naive(i), naive(ii)) and memory B cells, which were dominant in the blood sample that we included in our single-cell analysis. In contrast, the blood transcriptomic cluster appearing during flare (cluster 1) was enriched in the NR4A+, LMNA+, and naive(iii) clusters. This is consistent with a trajectory of B cell activation whereby blood B cells are recruited to the synovium pre-flare, activated within the synovial microenvironment to express NR4A, and then possibly re-circulate during flare. Re-circulation of B cells between the synovium and the periphery is also supported by a recent publication demonstrating shared clonal families between these compartments in RA (Elliott et al., 2020). The origin of the B cells increasing in the blood pre-flare remains to be determined. It also remains possible that the NR4A cluster genes enriched in the blood at the time of flare represent B cells activated elsewhere and recruited to the synovium.

B cells may contribute to RA pathogenesis in multiple ways, not only as the precursors of PCs and production of autoantibodies, but also by antibody-independent functions, including antigen presentation and cytokine production. B cells can produce cytokines and chemokines that promote ELS. Thus, lymphotoxin β expression by B cells has been shown to support ELS formation in a mouse colitis model (Lochner et al., 2011). Lymphotoxin β and α are elevated in RA synovium and relate to the level of inflammation in the tissue (O’Rourke et al., 2008), and LT-β is expressed by synovial B cells but also other cell populations (Takemura et al., 2001). The precise role of B cells in ELS formation in RA synovial tissue and which subpopulations may be critical are unclear. We found a significant relationship between NR4A B cell expression and the degree of B cell infiltration and lymphoid organization in tissue, implicating NR4A B cells as a potential ELS-initiating and -promoting B cell subset. In accord with this, NR4A B cells expressed LT-β and LT-α and were the dominant source of IL-6, the latter important for Tfh interactions and demonstrated to promote spontaneous GC formation in SLE (Arkatkar et al., 2017).

Overall, our results highlight the unique enrichment of NR4A-expressing B cells in the RA synovium and their identity as an antigen-experienced B cell in the synovial ELS with the capability of differentiating to PC or memory B cell in situ depending on additional signals. NR4A may serve as a biomarker of chronic autoantigen stimulation in the synovial microenvironment and a novel therapeutic target.

Limitations of the study

Although our study suggests that B cells are activated in situ in ELS, we acknowledge that the evidence is indirect. In addition, particularly the BCR sequencing data will need to be validated in a larger number of patients with more paired blood samples. The RA patients in our study were selected based on the presence of ELS. Although we found that NR4A B cell expression is enriched in leukocyte-rich RA in a larger more histologically and clinically diverse cohort, this reference dataset used bulk sequencing and did not examine the BCR. Despite these limitations, the enrichment of B cell NR4A1 expression in the RA synovium suggests chronic antigen stimulation. However, the antigen specificity of these cells remains to be determined. It is also unclear how the strength and duration of antigen stimulation regulates NR4A1 expression in synovial B cells and the functional outcome. Indeed, the different NR4A receptors influence T cell developmental fate via differential sensitivity to antigen receptor stimulation (Jennings et al., 2020). Furthermore, how other co-stimulatory signals control B cell activation and differentiation to PCs in the synovium and the role of NR4A in these processes requires further study. We speculate that the persistent high expression of NR4A in synovial PCs is indicative of antigen selection and differentiation in the tissue, but this remains to be defined. Finally, it will be important in future work to define the B cell migration dynamics between synovium and blood and the potential role of NR4A B cells in flare by comparing blood and tissue longitudinally in a larger number of patients.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Jennifer Anolik (Jennifer_anolik@urmc.rochester.edu).

Materials availability

This study did not generate new unique reagents.

Data and code availability

Single-cell RNA-seq data have been deposited at GEO depository. The records have been assigned GEO accession numbers and are publicly available as of the date of publication. This paper analyzes existing publicly available data. These accession number and references for the datasets are listed in the key resources table.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
rat anti-human podoplanin monoclonal antibody, PE	eBioscience	Cat#12-9381-42;RRID:AB_1582262
PerCP-cy5.5 mouse anti-human CD90	BD Biosciences	Cat#561557;RRID:AB_10712762
BV786 mouse anti-human CD45	BD Biosciences	Cat#563716; RRID: AB_2716864
Alexa Fluor 700 mouse anti-human CD14	BD Biosciences	Cat#557923; RRID: AB_396944
FITC mouse anti-human IgD	BD Biosciences	Cat#555778; RRID: AB_396113
PE/Dazzle 594 anti-human CD19	Biolegend	Cat#302251; RRID:AB_2,563,559
PE-cy5 mouse anti-human CD3	BD Bioscience	Cat#555334; RRID: AB_395741
Brilliant Violet 605 anti-human CD27	Biolegend	Cat#302830; RRID:AB_2561450
APC/Cyanine7 anti-human CD4	Biolegend	Cat#317418; RRID:AB_571947
PE/Cyanine7 anti-human CD11c	Biolegend	Cat#301608; RRID:AB_389351
Alexa Fluor 647 anti-human CD24	Biolegend	Cat#311109; RRID:AB_528783
APC/Cy7 mouse anti-human CD19	BD Bioscience	Cat#557791; RRID: AB_396873
Alexa Fluor 700 anti-human CD3	Biolegend	Cat#300424; RRID:AB_493741
BUV737 mouse anti-human IgD	BD Biosciences	Cat#612798;RRID: AB_2738894
PerCP-cy5.5 mouse anti-human CD38	BD Biosciences	Cat#551400;RRID: AB_394184
mouse anti-human Nur77 antibody PE	eBioscience	Cat#12-5965-82;RRID:AB_1257209
BV421 mouse anti-human CD183	BD Biosciences	Cat#562558;RRID: AB_2737653
PE/Dazzle 594 anti-human/mouse BCL6	Biolegend	Cat#358510;RRID:AB_2566194
PE-cy5 mouse anti-human CD69	BD Biosciences	Cat#555532;RRID: AB_395917
Alexa Fluor 647 anti-human CD83	Biolegend	Cat#305316;RRID:AB_2076531
BV786 mouse anti-human IgG	BD Biosciences	Cat#564230;RRID:AB_2738684
BUV395 mouse anti-human Ki-67	BD Biosciences	Cat#564710;RRID: AB_2738577
FITC mouse anti-human CD77	BD Biosciences	Cat#551353;RRID: AB_394165
PE-cy7 mouse anti-human CD184	BD Biosciences	Cat#560669;RRID: AB_1727435
IHC-plus polyclonal goat anti-human CD20	LSBio	Cat#LS-B11144
rabbit anti-human NR4A1 polyclonal antibody	Sigmaaldrich	Cat#HPA070142;RRID:AB_2732149
mouse anti-human CD21 monoclonal antibody	Invitrogen	Cat#MA5-11417;RRID:AB_10982851
monoclonal mouse anti-human Ki-67	Agilent Dako	Cat#M7240;RRID:AB_2142367
Mouse anti-CD138	LSBio	Cat# LS-B9360;RRID: AB_2877650
Goat anti-PCNA	Santa Cruz Biotechnology	Cat# Sc-9857; RRID:AB_2160372
Donkey anti-Goat IgG (H + L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 568	Invitrogen	Cat#A-11057;RRID:AB_2534104
Alexa Fluor® 488 AffiniPure F(ab’)2 Fragment Donkey Anti-Rabbit IgG (H + L)	Jackson ImmunoResearch Lab	Cat#711-546-152;RRID: AB_2340619
Alexa Fluor® 647 AffiniPure F(ab’)2 Fragment Donkey Anti-Mouse IgG (H + L)	Jackson ImmunoResearch Lab	Cat#715-606-150;RRID: AB_2340865
AffiniPure F(ab’)2 Fragment Goat Anti-Human IgM, Fc5μ fragment specific	Jackson ImmunoResearch Lab	Cat#109-006-129;RRID: AB_2337553
AffiniPure F(ab’)2 Fragment Goat Anti-Human IgA + IgG + IgM (H + L)	Jackson ImmunoReasearch Lab	Cat#109-006-064;RRID: AB_2337548
Biological Samples
Human synovial tissue and blood	Hospital for Special Surgerry, AMP Network (ImmPort SDY998)	N/A
Human tonsil	Surgical Pathology, University of Rochester Medical Center	N/A
Critical commercial assays
LIVE/DEAD™ Fixable Aqua Dead Cell Stain Kit, for 405 nm excitation	Invitrogen	Cat#L34966
Foxp3/Transcription Factor Staining Buffer Set	eBioscience	Cat#00-5523-00
Chromium Single-Cell5′ Library & Gel Bead Kit	10xGenomics	Cat#PN-1000006; Cat# PN-1000014,
the Chromium Single Cell V(D)J Enrichment Kit, Human B cell	10xGenomics	Cat#PN-1000016
Deposited data
PEAC RNA-seq data	Lewis et al. (2019)	ArrayExpress E-MTAB-6141; https://peac.hpc.qmul.ac.uk/
RNA-seq of RA flares	Orange et al. (2020)	https://www.nejm.org/doi/full/10.1056/NEJMoa2004114
SLE Kidney single cell RNA seq	Arazi et al. (2019)	ImmPort SDY997
RA and OA RNA-seq	Zhang et al. (2019)	ImmPort SDY1299
single cell RNA-seq of peripheral B cells	Zheng et al. (2017)	SRP073767; https://support.10xgenomics.com/single-cellgene-expression/datasets/1.1.0/b_cells.
Raw and analyzed data single cell RNA and BCR sequencing of RA synovial tissue and blood	This paper	GEO196150
Blood transcription modules	Li et al. (2014)	GSE52245, GSE13485, GSE29617, GSE29615
Tonsil DZ and LZ gene sets	Victora et al., 2012	GSE38696 and GSE38697
RA synovial tissue single cell RNA-seq	Zhang et al. (2019)	ImmPort SDY998
PreM gene set	Holmes et al. (2020)	GSE139891
BCR sequences	Human immunoglobin reference	The international ImMunoGeneTics information system; http://www.imgt.org/
Code for data analyses	This paper	DOI: 10.5281/zenodo.6000031; https://zenodo.org/record/6000031
Oligonucleotides
NR4A1(Hs00374226_m1)	Applied Biosystems	Cat#4331182
NR4A2 (Hs01117527_g1)	Applied Biosystems	Cat#4331182
NR4A3(Hs00545009_g1)	Applied Biosystems	Cat#4331182
PPIA(Hs04194521_s1)	Applied Biosystems	Cat#4448489
Software and algorithms
Cellranger	N/A	https://support.10xgenomics.com/single-cell-gene-expression/software/overview/welcome
SeqGeq v1.4	N/A	https://www.flowjo.com/
Seurat 2.3.4	Butler et al. (2018)	https://satijalab.org/seurat/
scran 1.8.4	N/A	http://packages.renjin.org/package/org.renjin.bioconductor/scran
CellaRepertorium version 0.8.1	N/A	http://bioconductor.org/packages/release/Bioc/html/CellaRepertorium.html
DESeq2 version 1.20.0	N/A	https://bioconductor.org/packages/release/bioc/html/DESeq2.html
clusterProfiler 3.8.1	N/A	https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html
Graphpad Prism	N/A	https://www.graphpad.com/
SingleR	N/A	https://bioconductor.org/packages/release/bioc/html/SingleR.html
FlowJo v.10	N/A	https://www.flowjo.com/
HighV-QUEST	Alamyar et al.,2012	http://www.imgt.org/IMGTindex/IMGTHighV-QUEST.php

All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Patient cohorts

Arthroplasty samples were acquired from RA patients after surgical removal as part of standard of care at Hospital for Special Surgery, NY (n = 3, 1 male, 2 female) or collected as part of the AMP network (ImmPort SDY998) (n = 1, female). Patient disease characteristics can be found in Table S1. Peripheral blood was collected prior to the procedure. Written informed consent for research use of patient samples was obtained prior to study inclusion at the time of sample collection. This study received approval from the institutional review board (IRB) at Hospital for Special Surgery and University of Rochester. Synovial tissue fragments were cryopreserved in CryoStor CS10 (BioLife Solution) for subsequent disaggregation at University of Rochester (Donlin et al., 2018). PBMC were isolated from whole blood by Ficoll-hypaque density gradient centrifugation and cryopreserved in CryoStor CS10 (BioLife Solution). Samples from the AMP network were collected as part of multi-center cross-sectional study of patients undergoing elective surgical procedures and a prospective observational study of synovial biopsy specimen from patients with RA aged ≥18 years, with at least one inflamed joint (Zhang et al., 2019). This cohort included synovial tissues from 36 RA biopsies and arthroplasties (29 female, 7 male) and 15 OA arthroplasties (10 female, 5 male). 14 RA biopsies have B cell single cell RNA sequencing data available for analysis. The Pathobiology of Early Arthritis cohort (PEAC) contains data from 87 RA patients (64 female, 23 male) fulfilling 2010 ACR/EULAR RA Classification Criteria that were enrolled as part of the Medical Research Council funded multi-center PEAC project (Lewis et al., 2019) (https://peac.hpc.qmul.ac.uk/). Longitudinal genomic analysis of 67 blood samples over 13 weeks from one female RA patient sampled before and during flare was analyzed from Orange et al., 2020.

Human subjects for PBMCs, and synovial cells for in vitro study

Normal control blood (n = 8, all female), RA peripheral blood (n = 6, 4 female, 2 male) and RA synovial cells (n = 9, 7 female, 2 male) were obtained from patients at University of Rochester under IRB study protocols approval from University of Rochester. PBMCs and synovial cells were isolated by Ficoll-hypaque density gradient centrifugation and cryopreserved in freezing medium (90%FBS/10% DMSO) until used.

Tonsil tissue samples

Tonsil tissue samples were obtained as discarded tissue from tonsillectomy procedures. Patient information including gender is not available because tissue was obtained as de-identified sample. Portions of tissue was fixed in formalin and paraffin-embedded for histologic staining. The remaining tissue was mechanically disaggregated to obtain cell suspension. The cell suspension was further purified by Ficoll-hypaque density gradient centrifugation. Isolated mononuclear tonsil cells were cryopreserved in freezing medium (90%FBS/10%DMSO) until used.

METHOD DETAILS

Preparation of synovial tissue and PBMC and flow cytometry cell sorting

Cryopreserved synovial fragments were thawed per SOP (Donlin et al., 2018). Synovial tissue fragments were then disaggregated by combination of an enzymatic digestion and mechanical disruption as described. PBMCs were thawed by rapidly warming the cryovial in a 37°C water bath. To stain cells for sorting, cells were first incubated with Fc-Block (Biolegend) for 5 min at room temperature (RT). Cells were then incubated with a mixture of fluorochrome-conjugated antibodies (Table S3) for 30 min on ice, washed with PBS and then resuspended in PBS solution containing DAPI (Biolegend). Cells were sorted on a FACSAria II (BD biosciences) into 5%FBS/RPMI media using a 100-μm nozzle. B cells were gated as CD45+CD3−CD14−CD19+ (Figure S1A).

Single cell RNA and BCR repertoire sequencing

Cellular suspensions were loaded on a Chromium Single-Cell Instrument (10× Genomics, Pleasanton, CA, USA). Single-cell RNA-Seq libraries were prepared using Chromium Single-Cell 5′ Library & Gel Bead Kit (PN-1000006, PN-1000014, 10× Genomics). The barcoded, full-length V(D)J segments were enriched by PCR with primers specific to the BCR constant regions using the Chromium Single Cell V(D)J Enrichment Kit, Human B cell (PN-1000016, 10× Genomics). For both B-cell enriched library and 5′ gene expression library, enzymatic fragmentation and size selection was used to optimize the cDNA amplicon size and indexed sequencing libraries were constructed by End Repair, A-tailing, Adaptor Ligation, and PCR. Final libraries contain the P5 and P7 priming sites used in Illumina bridge amplification. Paired end reads of 150nt were generated for each B-cell enriched library using the MiSeq sequencer (Illumina). The 5′ gene expression libraries were sequenced following 10× Genomics read length guidelines on the NovaSeq6000 sequencer (Illumina). The samples were disaggregated, sorted and loaded on 10× in multiple batches. Library preparations and sequencing for all samples were done in one batch.

Intracellular staining and flow cytometry

To stain for flow cytometry, PBMC or tonsil cells were first incubated with Fc-Block (Biolegend) for 5 min at RT. Then 100ul of surface antibody cocktail (Table S3) was added to each tube and incubated on ice for 30 min. Next, live/dead cell assay (Invitrogen) was applied for 15 min on ice. Intracellular staining was performed using Foxp3/Transcription factor staining buffer kit (eBioscience) following manufacture protocol. After intracellular staining, cells were washed and fixed in 100μL PBS/1% paraformaldehyde for 20 min at RT and then 200μL of PBS/5%BSA was added to dilute the fixative. Samples were analyzed on an 18-color LSR II flow cytometer (eBioscience) on the same day or after kept at 4°C overnight. Analyses were performed using FlowJo v.10 (Tree Star) and doublet exclusion was performed on all samples. Gating strategy is depicted in Figure S10.

In vitro stimulation

Freshly isolated or thawed PBMC were washed with complete medium and counted. For qPCR, B cells were positively selected using anti-CD19 magnetic beads (Miltenyi Biotech). 0.5–1.0 × 106 cells per well were seeded in round-bottom 96-well plates in 200μL 10% FBS/RPMI 1640. Cells were stimulated with anti-human IgM F(ab′2) or anti-human Ig(A + G + M) F(ab′2) (Jackson Immuno-Research Laboratories) for indicated times. After stimulation, cells were collected for flow cytometry or qPCR.

QUANTITATIVE REAL-TIME PCR (qPCR)

B cells were collected in RLTplus lysis buffer (Qiagen) supplemented with 1% 2-mercaptoethanol, and RNA were extracted using Quick-RNA kit (Zymo Research) and concentration was measured by NanoDrop (Thermo Scientific). Complementary DNA (cDNA) synthesis was done using a qScript cDNA synthesis kit (Quantabio). TaqMan expression assay system for NR4A1 (Hs00374226_m1), NR4A2 (Hs01117527_g1), NR4A3 (Hs00545009_g1) from Applied Biosystems were used to detect NR4A1, NR4A2 and NR4A3 mRNA expression, respectively. PPIA (Hs04194521_s1, Applied Biosystems) was used as a control housekeeping gene. The fold change in mRNA expression relative to control was calculated using the comparative threshold cycle (ΔΔCt) method.

Immunofluorescent staining

Primary antibodies: List of primary and secondary antibodies used for immunofluorescent staining can be found in Table S3. 5 μm formalin-fixed paraffin sections were incubated at 60°C overnight (ON). Tissue sections were quickly transferred to xylenes and gradually hydrated by transferring slides to absolute alcohol, 96% alcohol, 70% alcohol, and water. Slides were immersed in an antigen retrieval solution, boiled for 30 min, and cooled down for 10 min at room temperature (RT). Slides were rinsed several times in water and transferred to PBS. Non-specific binding was blocked with 5% normal donkey serum in PBS containing 0.1% Tween 20, 0.1% Triton X-100 for 30 min, at RT in a humid chamber. Primary antibodies were added to slides and incubated in a humid chamber at room RT, ON. Slides were quickly washed in PBS, and fluorescently labeled, secondary antibodies were incubated for 2 h at RT in a humid chamber. Finally, slides were rinsed for 1 h in PBS and mounted with Vectashield antifade mounting media with DAPI (H-1200, Vector Laboratories). Images were taken with a Zeiss Axioplan 2 microscope and recorded with a Hamamatsu camera.

QUANTIFICATION AND STATISTICAL ANALYSIS

Statistical analysis

The statistical tests used are as described in each figure legend. p values ≤ 0.05 are considered significant.

Single cell RNA and repertoire sequencing analysis

Raw reads generated from the Illumina basecalls were demultiplexed using Cell Ranger mkfastq version 2.1.1 utilizing bcl2fastq version 2.19.1.403. Cells were aligned and counted against the Cell Ranger reference refdata-cellranger-GRCh38–1.2.0 using Cell Ranger count version 2.1.1, utilizing only the library index, project, and sample ID for a given sample. We performed quality control to remove outlier events and doublets based on the cell’s total UMI size versus the number of distinct genes expressed in SeqGeq™ v1.4 (FlowJo, LLC) by gating on bivariate scatter plots. We also filtered putatively degraded cells based on the proportion of UMIs in the cell derived from the mitochondrial genome. These steps removed 252 cells. After filtering, the data were exported for further analysis in scran 1.8.4 and Seurat 2.3.4 (Butler et al., 2018). The top 1000 variable genes, excluding IGHV, IGKV and IGLV genes were used to derive 6 principal components. These were used to define clusters using the Louvain method and resolution parameter 0.5. Based on expression of non-B cell markers (COL1A2+, FN1+, MS4A1−), we identified 77 putative fibroblasts, and accordingly excluded these cells from further analyses. Differentially expressed genes (DEGs) were identified by calculating the difference between the average expression by cells in the subset versus cells not in the subset. BCR clonotypes were assigned by requiring 97% identity of the DNA sequence of the CDR3 region of the heavy or light chain using CellaRepertorium version 0.8.1 (http://bioconductor.org/packages/release/Bioc/html/CellaRepertorium.html). For somatic hypermutation (SHM) rate calculation, the consensus FASTQ sequences were aligned to the human immunoglobulin reference from the international ImMunoGeneTics information system (IMGT) using HighV-Quest (Alamyar et al., 2012) to calculate the percent identity in V and J in heavy and light chains which then averaged, weighting by the length of each segment.

Projection of clustering using SingleR

Gene-cell count matrices were downloaded for the AMP Phase I RA and SLE cohorts (ImmPort accessions SDY998 and SDY997). Filtered gene-cell count matrices for CD19+ B cells from Zheng et al., (2017) (Zheng et al., 2017) were downloaded from https://support.10xgenomics.com/single-cell-gene-expression/datasets/1.1.0/b_cells. 12848 gene symbols expressed in common in these datasets and our data were retained and the data merged. QC was evaluated using the quickPerCellQC function in scran and 874 cells were filtered. A SingleR model was trained on our data and its 8 cluster labels, and the inferred cluster identities in the other datasets were queried. Binomial mixed models were used to estimate the frequency of the NR4A+ cluster among these other data sets. For bulk RNAseq from AMP RA Cohort, data was downloaded (ImmPort accession SDY1299). DESeq2 version 1.20.0 was used to evaluate differential expression of NR4A+ marker genes in sorted B cells from this study.

Gene set enrichment analysis

For GSEA analysis, we used gene set from Blood transcription modules from Li et al. (2014) (Li et al., 2014), bulk RNA sequencing of sorted tonsil B cells (Bm1+Bm2 and Bm5) and DZ and LZ GC gene lists from Victora et al. (2012) (Victora et al., 2012). We used clusterProfiler 3.8.1 to conduct GSEA analysis. Heatmaps displaying the signed −log10(FDR q value) of the result was created using GraphPad Prism (GraphPad Software, LLC).

Comparison to longitudinal RNA sequencing

Orange et al. (2020) identified 2461 genes that were differentially expressed across the longitudinal time course. We took the intersection between these differentially expressed genes and marker genes for the synovial subpopulations. We calculated the relative expression of each gene across the time longitudinal time course as , where Y was batch-corrected expression in log2 reads per kilobase mapped of gene i at week t, Y. is the median expression of gene i over all weeks and subjects, and m is the median absolute deviation, and 0.6745 is a factor that scales the statistic to be asymptotically equal to a Z-score for normally distributed data. Figure 7C shows Z for the intersection between marker genes of our subpopulations and the 2,461 differentially genes from Orange et al. (2020). Figure S9B shows Z for these same sets of selected genes as a heatmap. Orange et al. (2020) defined 5 sets of genes: Cluster1, AC2, AC3, Cluster4, Cluster5. We tested for enrichment between these 5 sets and markers from our 9 subpopulations using hypergeometric tests. The −log10 transformed FDR q-values from these enrichment tests and the number of overlapping genes vs the number genes in the Orange et al. (2020) sets are shown in Figure S9A.