Defining inflammatory cell states in rheumatoid arthritis joint synovial tissues by integrating single-cell transcriptomics and mass cytometry

Abstract

To define the cell populations that drive joint inflammation in rheumatoid arthritis (RA), we applied single-cell RNA sequencing (scRNA-seq), mass cytometry, bulk RNA sequencing (RNA-seq) and flow cytometry to T cells, B cells, monocytes, and fibroblasts from 51 samples of synovial tissue from patients with RA or osteoarthritis (OA). Utilizing an integrated strategy based on canonical correlation analysis of 5,265 scRNA-seq profiles, we identified 18 unique cell populations. Combining mass cytometry and transcriptomics revealed cell states expanded in RA synovia: THY1(CD90)+HLA-DRAhi sublining fibroblasts, IL1B+ pro-inflammatory monocytes, ITGAX+TBX21+ autoimmune-associated B cells and PDCD1+ peripheral helper T (TPH) cells and follicular helper T (TFH) cells. We defined distinct subsets of CD8+ T cells characterized by GZMK+, GZMB+, and GNLY+ phenotypes. We mapped inflammatory mediators to their source cell populations; for example, we attributed IL6 expression to THY1+HLA-DRAhi fibroblasts and IL1B production to pro-inflammatory monocytes. These populations are potentially key mediators of RA pathogenesis.

Conflict of interest statement

COMPETING FINANCIAL INTERESTS

The authors declare no competing financial interests.

Figures

Figure 1.

Overview of synovial tissue workflow and pairwise analysis of high-dimensional data. a. We acquired synovial tissue, disaggregated the cells, sorted them into four gates representing fibroblasts (CD45−CD31−PDPN+), monocytes (CD45+CD14+), T cells (CD45+CD3+), and B cells (CD45+CD3−CD19+). We profiled these cells with mass cytometry, flow cytometry, sorted low-input bulk RNA-seq, and single-cell RNA-seq. Here, we use Servier Medical Art by Servier for the joint picture. b. Presence and absence of five different data types for each tissue sample. c. Schematic of each dataset and the shared dimensions used to analyze each of the three pairs of datasets with canonical correlation analysis (CCA). d. CCA finds a common mapping for two datasets. For bulk RNA-seq and single-cell RNA-seq, we first find a common set of g genes present in both datasets. Each bulk sample si gets a coefficient ai and each cell ci gets a coefficient bi. The linear combination of all samples s1…n arranges bulk genes along the canonical variate CVs1 and the linear combination of all cells c1…m arranges single-cell genes along CVc1. CCA finds the coefficients a1…n and b1…m that arrange the genes from the two datasets in such a way that the correlation between CVs1 and CVc1 is maximized. After CCA finds the first pair of canonical variates, the next pair is computed on the residuals, and so on.

Figure 2.

Distinct cellular composition in synovial tissue from OA, leukocyte-poor RA, and leukocyte-rich RA patients. a. Histological assessment of synovial tissue derived from OA (n = 15 independent tissue samples), leukocyte-poor RA (n = 17 independent tissue samples), and leukocyte-rich RA (n = 19 independent tissue samples). b. Cellular composition of major synovial cell types by flow cytometry. c. Synovial T cells, B cells, and monocytes by flow cytometry in samples from OA (n = 15), leukocyte-poor RA (n = 17), and leukocyte-rich RA (n = 19). Leukocyte-rich RA tissues were significantly higher infiltrated in synovial T cells (Student’s one-sided t-test P = 4×10−9, t-value = 8.92, df = 22.27) compared to leukocyte-poor RA and OA. Leukocyte-rich RA tissues were significantly higher infiltrated in synovial B cells (Student’s one-sided t-test P = 1×10−3, t-value = 3.50, df = 20.56) compared to leukocyte-poor RA and OA. Center value is mean. Statistical significance levels: ****P<1×10−4 and ***P<1×10−3. d. Quantitative histologic inflammatory scoring of both sublining cell layer and lining layer. Leukocyte-rich RA samples (n = 19) exhibited higher (Student’s one-sided t-test P = 1×10−3, t-value = 3.21, df = 30.66) Krenn inflammation scores than leukocyte-poor RA (n=15) and OA tissues (n = 10) samples. Center value is mean. e. Correlation between leukocyte infiltration assessed by cytometry with histologic inflammation score (n = 44 biologically independent samples). Student’s one-sided t-test P = 3×10−09, t-value = 7.15, df = 46.51. f. tSNE visualization of synovial cell types in OA, leukocyte-poor RA, and leukocyte-rich RA by mass cytometry density plot.

Figure 3.

High-dimensional transcriptomic scRNA-seq clustering reveals distinct cell type subpopulations. a. 18 clusters across 5,265 cells from all cell types on a tSNE visualization. b. Cluster abundances across donors. c. Fibroblasts: three types of THY1+ sublining fibroblasts (SC-F1, SC-F2, and SC-F3) and CD55+ lining fibroblasts (SC-F4). d. Monocytes: two activated cell states of IL1B+ pro-inflammatory (SC-M1) and IFN-activated (SC-M4) monocytes. e. T cells: CD4+ subsets: SC-T1, SC-T2, SC-T3, and CD8+ subsets: SC-T4, SC-T5, and SC-T6. f. B cells: HLA+ (SC-B1, SC-B2, and SC-B3) and plasmablasts (SC-B4). The cluster colors in c-f are consistent with (a).

Figure 4.

Distinct synovial fibroblast subsets defined by cytokine activation and MHC II expression. a. scRNA-seq analysis identified three sublining subsets, CD34+ (SC-F1), HLAhi (SC-F2), and DKK3+ (SC-F3) and one lining subset (SC-F4). Differential analysis between leukocyte-rich RA (n = 16) and OA (n = 12) bulk RNA-seq fibroblast samples shows marker genes upregulated or downregulated in leukocyte-rich RA. Fold changes with 95% confidence interval (CI). b. By querying the leukocyte-rich RA (n = 16) and OA (n = 12) fibroblast bulk RNA-seq samples, scRNA-seq cluster HLA-DRAhi (SC-F2) and CD34+ (SC-F1) fibroblasts are significantly overabundant (two-sided Student’s t-test P=2×10−6, t-value=6.2, df = 23.91 and P=2×10−3, t-value = 3.20, df = 25.41, respectively) in leukocyte-rich RA relative to OA. Lining fibroblasts (SC-F4) are overabundant (two-sided Student’s t-test P=5×10−7, t-value=−5.31, df =21.97) in OA samples. Fold changes with 95% CI. c. Pathway enrichment analysis for each cluster. Two-sided Kolmogorov-Smirnov test with 105 permutations; Benjamini-Hochberg FDR is shown. d-e. Identified subpopulations from fibroblasts (n = 25,161) and disease status from 6 leukocyte-rich RA, 9 leukocyte-poor RA, and 8 OA by mass cytometry on the same gating with scRNA-seq. f-g. Normalized intensity of distinct protein markers shown in tSNE visualization and averaged for each cluster heatmap. h. CCA projections of mass cytometry clusters and bulk RNA-seq genes. First two canonical variates (CVs) separated genes upregulated in leukocyte-rich RA from genes upregulated in OA. HLAhi genes are highly associated with THY1+CD34−HLA-DRhi by mass cytometry. i. Integration of mass cytometry clusters with scRNA-seq clusters based on the top markers (AUC > 0.7) for each scRNA-seq cluster using top 10 canonical variates in the low-dimensional CCA space. We computed the spearman correlation between each pair of scRNA-seq cluster and mass cytometry cluster in the CCA space and performed permutation test 104 times. Z-score is calculated based on permutation p-value. We observed HLAhigh sublining fibroblasts by scRNA-seq are strongly correlated with THY1+CD34−HLA-DRhi fibroblasts by mass cytometry.

Figure 5.

Unique activation states define synovial monocytes heterogeneity. a. scRNA-seq analysis identified four subsets: IL1B+ pro-inflammatory monocytes (SC-M1), NUPR1+ monocytes (SC-M2) with a mixture of leukocyte-poor RA and OA cells, C1QA+ (SC-M3), and IFN-activated monocytes (SC-M4). Differential analysis by bulk RNA-seq on leukocyte-rich RA samples (n = 17) and OA samples (n = 13) revealed upregulation/downregulation of cluster marker genes. Effect sizes with 95% CI are given. b. By querying the bulk RNA-seq, we found scRNA-seq cluster IL1B+ pro-inflammatory monocytes (two-sided Student’s t-test P=6×10−5, t-value=4.56, df =26.33) and IFN-activated monocytes (two-sided Student’s t-test P=6×10−3, t-value=3.28, df =23.68) are upregulated in leukocyte-rich RA (n = 17) compared to OA (n = 13), while SC-M2 is depleted (two-sided Student’s t-test P=2×10−5, t-value=−5.62, df=26.81) in leukocyte-rich RA. Error bars indicate mean and 95% CI. c. Pathway enrichment analysis indicates the potential pathways for each subset. Two-sided Kolmogorov-Smirnov test with 105 times permutation was performed; Benjamini-Hochberg was used to control the FDR of multiple tests. The standard names for the immunological gene sets from up to bottom are: Genes down-regulated in neutrophils versus monocytes (GSE22886); Genes down-regulated in healthy myeloid cells versus SLE myeloid cells (GSE10325); Genes down-regulated in control microglia cells versus those 24 h after stimulation with IFNG (GSE1432); Genes down-regulated in unstimulated macrophage cells versus macrophage cells stimulated with LPS (GSE14769); Genes up-regulated monocytes treated with LPS versus monocytes treated with control IgG (GSE9988); Genes up-regulated in monocytes versus myeloid dendritic cells (mDC) (GSE29618); Genes up-regulated in monocytes versus plasmacytoid dendritic cells (pDC) (GSE29618). d. Detection of pro-inflammatory IL-1β in inflamed synovium by multicolor immunofluorescent staining with antibodies CD14 (red), IL-1β (green), and counterstained with DAPI (blue) identified CD14+IL-1β+ cells (white arrow). The experiment was repeated > 5 times with staining of 6 independent leukocyte-rich RA samples with similar results. Image was acquired at 200 magnification. Scale bar is 50 μm. e–f. Identified subpopulations from monocytes (n = 15,298) and disease status from 6 leukocyte-rich RA, 9 leukocyte-poor RA, and 11 OA by mass cytometry on the same gating with scRNA-seq. g-h. Normalized intensity of distinct protein markers by tSNE visualization and averaged for each cluster in heatmap. i. Integration of identified mass cytometry clusters with bulk RNA-seq reveals genes that are associated with CD11c+CD38+ and CD11c+CCR2+, like IFITM3, CD38, HBEGF, ATF3, and HLA+ genes. j. Integration of mass cytometry clusters and scRNA-seq clusters revealed that CD11c+CD38+ by mass cytometry are significantly associated with IL1B+ pro-inflammatory (SC-M1) monocytes.

Figure 6.

Synovial T cells display heterogeneous CD4 and CD8 T cell subpopulations in RA synovium. a. scRNA-seq analysis identified three CD4+ subsets: CCR7+ (SC-T1), Treg cells (SC-T2), and Tph and Tfh (SC-T3); and three CD8+ subsets: GZMK+ (SC-T4), CTLs (SC-T5), and GZMK+GZMB+ (SC-T6). Differential expression analysis on leukocyte-rich RA (n = 18) comparing with OA (n = 13) on sorted T cell bulk RNA-seq samples revealed that CXCL13 is most significantly enriched in leukocyte-rich RA compared to OA. Effect sizes with 95% CI are given. b. Disease association of scRNA-seq clusters by aggregating top markers (AUC>0.7) by comparing leukocyte-rich RA (n = 18) with OA (n = 13) using bulk RNA-seq. Tph and Tfh cells (SC-T4) are upregulated (two-sided Student’s t-test p=0.01, t-value=2.73, df =29.00) in leukocyte-rich RA. Error bars indicate mean and 95% CI. c. Pathway analysis based on immunologic gene set enrichment indicates the potential enriched T cell states pathways. Two-sided Kolmogorov-Smirnov test with 105 times permutation was performed; Benjamini-Hochberg was used to control the FDR of multiple tests. The brief description of the standard names from up to bottom are: Genes up-regulated in CD4 high cells from thymus: Treg versus T conv (GSE42021); Genes up-regulated in comparison of effector CD8 T cells versus memory CD8 T cells (GOLDRATH); Genes down-regulated in comparison of effector memory T cells versus central memory T cells from peripheral blood mononuclear cells (PBMC) (GSE11057); Genes up-regulated in comparison of effective memory CD4 T cells versus Th1 cells (GSE3982); Genes up-regulated in comparison of T follicular helper (Tfh) cells versus Th17 cells (GSE11924). d. Detection of CD3+CD8+IFNγ+ (white arrow) in inflamed RA synovium by multicolor immunofluorescent staining with antibodies CD3 (green), CD8 (red), IFNγ (white), and counterstained with DAPI (blue). The experiment was repeated > 5 times with staining of 6 independent leukocyte-rich RA samples with similar results. Image was acquired at 200 magnification. Scale bar is 50 μm. e-f. Identified subpopulations from T cells (n = 19,985) and disease status from 6 leukocyte-rich RA, 9 leukocyte-poor RA, and 11 OA by mass cytometry. g-h. Distinct patterns of protein markers by tSNE and heatmap that define these clusters. i. Integration of identified mass cytometry clusters with bulk RNA-seq using CCA reveals bulk genes that are associated with CD4+PD-1+ICOS+ and CD8+PD-1−HLA-DR+ by mass cytometry. j. Integration of mass cytometry clusters with scRNA-seq clusters on the top markers (AUC>0.7) for each scRNA-seq cluster in the top 10 canonical variates. Z-score based on permutation test reveals that CD4+PD-1+ICOS+ and CD8+PD-1+HLA-DR+ by mass cytometry are highly associated with Tph and Tfh (SC-T3) by scRNA-seq; CD8+PD-1−HLA-DR+ T cells by mass cytometry are highly associated with CD8+ T cells (SC-T4, SC-T5, and SC-T6).

Figure 7.

Synovial B cells display heterogeneous subpopulations in RA synovium. a. scRNA-seq analysis identified naive B cells (SC-B1), memory B cells (SC-B2), autoimmune-associated B cells (ABCs) (SC-B3), and plasmablasts (SC-B4). Differential expression analysis is given by comparing leukocyte-rich RA (n = 16) with OA (n = 7) using bulk RNA-seq B cell samples. Effect size with 95% CI are given. b. Pathway enrichment analysis using immunologic gene sets indicates the distinct enriched pathways for each scRNA-seq cluster. Two-sided Kolmogorov-Smirnov test with 105 times permutation was performed; Benjamini-Hochberg was used to control the FDR of multiple tests. The standard names for the immunological gene sets from up to bottom are: Genes up-regulated in plasma cells versus memory B cells (GSE12366); Genes up-regulated in comparison of B cells versus plasmacytoid dendritic cells (pDC) (GSE29618); Genes up-regulated in B lymphocytes: naive versus plasmablasts (GSE42724); Genes up-regulated in B lymphocytes: human germinal center light zone versus dark zone (GSE38697); Genes up-regulated in comparison of memory IgM B cells versus plasma cells from bone marrow and blood (GSE22886); Genes up-regulated in comparison of memory IGG and IGA B cells versus plasma cells from bone marrow and blood (GSE22886). c. Detection of CD20+T-bet+CD11c+ (white arrow) in inflamed synovium by multicolor immunofluorescence. Immunofluorescent staining with antibodies CD20 (red), CD11c (white), T-bet (green), and counterstained with DAPI (blue). The experiment was repeated > 5 times with staining of 6 independent leukocyte-rich RA samples with similar results. Image was acquired at 200 magnification. Scale bar is 50 μm. d-e. Identified subpopulations of B cells (n = 8,179) and disease status from 6 leukocyte-rich RA, 9 leukocyte-poor RA, and 8 OA by mass cytometry. f-g. Distinct expression patterns of protein markers by tSNE and averaged for each cluster in heatmap. h. Integrating mass cytometry clusters with bulk RNA-seq data using CCA shows that CD38+CD20−Ig− (plasmablasts) population is highly associated with gene expression of plasma cells makers, like XBP1. i. Integration of mass cytometry clusters with scRNA-seq clusters suggested that CD38++CD20−IgM+HLA-DR+ and CD38++CD20−IgM−IgD− are significantly associated with plasmablast (SC-B4); IgM−IgD−HLA-DR++CD20+CD11c+ B cells are associated with ABCs (SC-B3).

Figure 8.

Transcriptomic profiling of synovial cells reveals upregulation of inflammatory pathways in RA synovium. a. Pathway enrichment using bulk RNA-seq identified shared and unique inflammatory response pathways for each cell type. Two-sided Kolmogorov-Smirnov test with 105 permutations was performed on 18 leukocyte-rich RA, 17 leukocyte-poor RA, and 14 OA. b. Bulk RNA-seq profiling of genes obtained from the significantly enriched pathways from (a) shows the averaged gene expression from each group (18 leukocyte-rich RA, 17 leukocyte-poor RA, and 14 OA) normalized across all cell type samples. c. scRNA-seq profiling resolved that inflammatory cytokines/chemokines, interferon responsive, and inflammatory responsive genes were driven by a global upregulation within a synovial cell type or discrete cell states.