Article Text
PDF
Abstract
Objectives Immune and stromal cell communication is central in the pathogenesis of rheumatoid arthritis (RA) and psoriatic arthritis (PsA), however, the nature of these interactions in the synovial pathology of the two pathotypes can differ. Identifying immune-stromal cell crosstalk at the site of inflammation in RA and PsA is challenging. This study creates the first global transcriptomic analysis of the RA and PsA inflamed joint and investigates immune-stromal cell interactions in the pathogenesis of synovial inflammation.
Methods Single cell transcriptomic profiling of 178 000 synovial tissue cells from five patients with PsA and four patients with RA, importantly, without prior sorting of immune and stromal cells. This approach enabled the transcriptomic analysis of the intact synovial tissue and identification of immune and stromal cell interactions. State of the art data integration and annotation techniques identified and characterised 18 stromal and 14 immune cell clusters.
Results Global transcriptomic analysis of synovial cell subsets identifies actively proliferating synovial T cells and indicates that due to differential λ and κ immunoglobulin light chain usage, synovial plasma cells are potentially not derived from the local memory B cell pool. Importantly, we report distinct fibroblast and endothelial cell transcriptomes indicating abundant subpopulations in RA and PsA characterised by differential transcription factor usage. Using receptor–ligand interactions and downstream target characterisation, we identify RA-specific synovial T cell-derived transforming growth factor (TGF)-β and macrophage interleukin (IL)-1β synergy in driving the transcriptional profile of FAPα+THY1+ invasive synovial fibroblasts, expanded in RA compared with PsA. In vitro characterisation of patient with RA synovial fibroblasts showed metabolic switch to glycolysis, increased adhesion intercellular adhesion molecules 1 expression and IL-6 secretion in response to combined TGF-β and IL-1β treatment. Disrupting specific immune and stromal cell interactions offers novel opportunities for targeted therapeutic intervention in RA and PsA.
- Autoimmunity
- Arthritis, Rheumatoid
- Arthritis, Psoriatic
Data availability statement
Data are available in a public, open access repository. Data are available upon reasonable request. All raw and processed files as well clinical information for each sample are deposited on national center for biotechnology information (NCBI), ascension number GSE200815 and are publically available without any restrictions of their subsequent use. Additionally, while detailed vignettes and base code is available on the bioinformatics platform GitHubgithub (as indicated in methods) for all packages used in the analysis, if specific parts of the code are needed, they will become available upon reasonable request.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See: http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ABOUT THIS SUBJECT?
Previous single cell RNA sequencing and flow cytometric analysis of sorted immune cells revealed the presence of peripheral helper T and follicular helper T cells and pathogenic B cells in the inflamed joint of patients with rheumatoid arthritis (RA).
THY1+ sublining synovial fibroblasts are expanded in RA.
Notch signalling driven synovial fibroblast and EC crosstalk contributes to synovial inflammation.
WHAT DOES THIS STUDY ADD?
First study to perform transcriptomic analysis of unsorted synovial tissue single cell suspensions of the RA and psoriatic arthritis (PsA) inflamed joint.
First time characterisation of immune-stromal cell interactions via the utilisation of receptor–ligand interaction networks on a global scale in the inflamed joint for RA and PsA.
Identification of differential fibroblast subpopulation involvement in PsA compared with RA.
Evidence regarding the origin of the T cell and B cell populations in the joint, potentially impacting the current paradigms.
HOW MIGHT THIS IMPACT ON CLINICAL PRACTICE OR FUTURE DEVELOPMENTS?
The data presented in this study will impact our understanding of RA and PsA synovial inflammation pathogenesis and will reveal new opportunities for targeted therapeutic intervention based on inhibition of specific cell–cell interactions.
Introduction
Rheumatoid arthritis (RA) and psoriatic arthritis (PsA) are common autoimmune and autoinflammatory diseases of unknown aetiology characterised by complex synovial pathology with a detrimental effect on the patient’s quality of life.1 2 RA and PsA are characterised by a spectrum of clinical manifestations that can be similar in both conditions, however, there are significant differences at a number of levels including clinical, anatomical, genetic, cellular and molecular.1–3 The most defined differences focus on the presence/absence of autoantibodies, synovial vascular morphology, the pattern of periarticular inflammation, bone erosion and new bone formation at the entheseal complex of peripheral and spinal joints.1–8 These differences may explain certain distinct clinical manifestations of the two diseases, and more importantly, may account for different responses to specific therapies impacting on disease outcome and prognosis.1–8 The complexity of synovial inflammation, associated with different pathotypes, is further increased by immune and stromal cell involvement.3–8 Recent implementation of single cell transcriptomic analysis of sorted synovial cells has revealed the diverse cellular landscape of the RA synovial stromal and immune cell compartments.9 While these studies have identified unique synovial cell clusters, increasing our understanding of potential pathogenetic mechanisms involved in RA, no studies to date have examined the synovial landscape of PsA, in addition to characterising the differential and complex immune-stromal cell crosstalk that may define the distinct synovial pathotypes observed in RA and PsA.
T cells have been implicated in RA and PsA synovial pathogenesis. Synovial polyfunctional CD4 and CD8 T cells expressing multiple pro-inflammatory cytokines simultaneously, associate with disease progression in PsA, with polyfunctional CD4 T cell responses recently reported in the synovial tissue of patients with RA.8 10 Synovial T cell functional plasticity is also highlighted by PD-1highCXCR5− peripheral helper T (Tph) cells sharing features with follicular helper T (Tfh) cells in promoting B cell antibody responses in RA.7 Clonally expanded CXCR3-expressing memory CD8 T cells with diverse phenotypes have been identified in the synovial fluid of patients with PsA and CD8 T cell clonal convergence between patients provides evidence for common MHC-I–antigen complex involvement and potential for T cell–stromal cell crosstalk.11 Along with T cells, macrophages are the predominant immune cells in synovial tissue. Macrophages form distinct subsets in the joint of patients with RA and exhibit immune regulatory and pro-inflammatory features, with MerTK+ macrophages associated with disease remission.12 Tissue-resident synovial macrophages characterised by expression of CX3CR1 form a physical barrier at the lining layer of the joint.13 These self-renewing macrophages have characteristics akin to epithelial cells and contribute to the homeostasis of the joint.13
Similarly to the diverse profile and roles of immune cells involved in synovial inflammation, emerging evidence suggests specific fibroblast cell subsets contribute to RA disease pathogenesis.6 FAPα and THY1 define two functionally distinct synovial fibroblasts with FAPα+ THY1+ fibroblasts mediating bone erosion whereas FAPα+ THY1− fibroblasts contribute to inflammation via the production of chemokines that promote immune cell trafficking to the inflamed joint.14
Despite recent advances in the resolution of the RA synovial tissue composition, several key questions remain unanswered, while additionally, the cellular landscape of PsA has not been explored at this level. To achieve precision medicine in RA and PsA, minimise lost time with exploratory treatments and reduce potential adverse effects, a better understanding of specific cell–cell interactions in RA and PsA is required. In this study transcriptomic analysis of intact RA and PsA synovial tissue cell suspensions was performed allowing for characterisation of immune-stromal cell interactions and the identification of overlapping and differential pathways of inflammation.
The cellular landscape of RA and PsA reveals points of convergence and distinct underlying mechanisms of synovial inflammation with utilisation of receptor–ligand interaction networks providing evidence of T cell and macrophage synergy in shaping the transcriptome of proinflammatory fibroblasts in RA.
Results
Single cell RNA sequencing reveals distinct synovial tissue immune and stromal cell clusters in patients with RA and PsA
Following the implementation of novel, strict, quality control measures (as described in the Methods section), we analysed a total of 178.196 cells derived from four patient with RA and five patient with PsA synovial tissue samples. Clinical characteristics of the patients at time of arthroscopic surgery are summarised in online supplemental table 1. Following data integration with Harmony to minimise sample to sample variation, cells were divided into nine megaclusters akin to distinct cell types. These megaclusters include: fibroblasts (88 953 cells), endothelial cells (24 207 cells), pericytes (4182 cells), macrophages (25 315 cells), dentritic cells (DC) (4103 cells), B cells (5902 cells), plasma cells (3098 cells), T cells (18 420 cells) and natural killer T cells (NKT) (1517 cells) (figure 1A). The annotation of the distinct cell type clusters was based on manual (prior knowledge-defined) and automated (scCATCH) identification of cell type-specific markers. Cell type-specific markers were identified following comparison of expression values of a specific cluster to all other synovial cells, leading to the generation of a list of cluster-specific markers (figure 1B).15 The identified cell type-specific megaclusters were divided further resulting in a total of 18 stromal cell clusters (11 fibroblast clusters, 6 endothelial cell clusters and 1 pericyte cluster) and 13 immune cell clusters (3 macrophage clusters, 5 T cell clusters and 2 B cell and 2 plasma cell clusters) (figure 1C). Distribution of expression on a single cell level using non-linear, stochastic Uniform Manifold Approximation and Projection for Dimensionality Reduction (UMAP) and expression level per cell type-specific megacluster of key markers is shown in figure 1D,E.
Supplemental material
High dimensionality single cell RNA sequencing analysis identifies specific cell clusters in patient with RA and PsA synovial tissue biopsies. (A) UMAP representation of 9 ‘mega’-clusters based on 178 196 cells across all cell types and synovial tissue biopsies (n=4 and 5 for patient with RA and PsA biopsies, respectively). (B) Differential gene expression analysis identifies upregulated or downregulated marker genes of the observed mega clusters. (C) Division of the nine identified mega clusters into a total of 33 subclusters. (D) Feature plots for the expression and distribution of the indicated genes in all cells. (E) Violin plots depicting log normalised expression per cluster of key markers used in cluster annotation. DC, dentritic cells NKT, natural killer T cells, IL, interleukin; PsA, psoriatic arthritis; RA, rheumatoid arthritis; UMAP, Uniform Manifold Approximation and Projection.
Differential fibroblast cluster distribution between PsA and RA synovial tissue samples
Recent studies have highlighted that patient with RA synovial fibroblasts are highly heterogeneous with a newly described synovial subset, characterised by expression of FAP and THY1, exhibiting effector function via the secretion of pro-inflammatory cytokines.14 Analysis of the fibroblast megacluster of 88.953 cells, resulted in the identification of 11 distinct fibroblast clusters in RA and PsA synovial biopsies (figure 2 and online supplemental figure S2). Importantly, there is differential abundance of the synovial fibroblast clusters (expressed as frequency of each cluster as part of all synovial fibroblasts per sample), separating PsA from patient with RA samples with a significantly higher abundance of F1 fibroblasts in PsA (***p<0.001) and F8, F9 and F11 fibroblasts in RA (**p=0.006, ***p<0.001, ***p<0.001, respectively) (figure 2A,D). We then examined the expression of FAP and THY1 by synovial fibroblast clusters, with the enriched RA F11 synovial fibroblast cluster harbouring the highest number of FAPα-expressing and/or THY1-expressing cells and the enriched in PsA F1 synovial fibroblast cluster showing almost no FAP or THY1 expressing cells (figure 2B). Expression of THY1 and FAP by cells of a specific fibroblast cluster is an under-representation of co-expressing cells due to dropouts in the sampling of RNA and sequencing. Therefore, to examine the degree of FAP and THY1 co-expression, we performed data imputation.16 Data imputation algorithms use the transcriptional profile of neighbouring cells to infer the expression of genes that may be affected by increased sparsity.17 Fibroblasts co-expressing FAP and THY1 are found in the enriched RA F11 cluster but not in the enriched PsA F1 cluster (figure 2C). Due to the differential abundance of clusters F11 and F1 between RA and PsA synovial tissue samples, we performed pathway enrichment analysis using differentially expressed genes between fibroblast cluster F11 and cluster F1. Pathway enrichment identified several pathways, previously implicated in RA pathogenesis, which are enriched in F11 compared with F1 fibroblasts. These pathways include the extracellular matrix (ECM) receptor, focal adhesion and RA pathways (figure 2E). Interestingly principal component analysis (PCA) of enriched pathways in synovial fibroblasts, revealed a separation of the fibroblast clusters, with F1 and F11 synovial fibroblasts on opposite ends of the spectrum (online supplemental figure S2C). Common genes and the upregulation or downregulation of specific members of these pathways in the comparison between F11 and F1 fibroblasts are shown in figure 2F.
Supplemental material
Distinct fibroblast cluster distribution in RA compared with PsA synovial biopsies. (A) Abundance of fibroblast clusters in patient with RA and PsA synovial biopsies. (B) Expression and percentage of positive cells per fibroblast cluster for FAP and THY1. (C) Scatterplots showing the relation between THY1 and FAP expressing cells before and after data imputation for RA and PsA fibroblast cluster F1 and fibroblast cluster F11. Fibroblast clusters with significantly different abundances between RA and PsA are indicated by green (higher in RA) and blue (higher in PsA) boxes. (D) Frequency of fibroblast clusters (calculated as a percentage of all fibroblasts/sample) in patient with PsA and RA synovial biopsies (n=4–5), data are presented as box and whiskers plots (min to max), symbols represent individual samples, statistical significance was determined by two-way analysis of variance with Sidak’s multiple comparisons test (***p<0.001, **p=0.0062). (E) Analysis of pathways enriched in fibroblast cluster 11 compared with fibroblast cluster 1, colour intensity represents significance and dot size the number of genes within each pathway that are differentially expressed. (F) Term plot of the indicated pathways with significant enrichment in fibroblast cluster 11 compared with fibroblast cluster 1. Colour indicates up or downregulation of specific genes within the pathway and dot size represents statistical significance of change. ECM, extracellular matrix; IL, interleukin; PsA, psoriatic arthritis; RA, rheumatoid arthritis.
Differential abundance and distinct transcriptional profile of specific endothelial cell clusters between PsA and RA
A pivotal first step of synovial inflammation in RA and PsA is increased angiogenesis which facilitates immune cell infiltration into the synovial tissue. Of the six endothelial and one pericyte cell cluster, endothelial cell cluster E1 is significantly (*p=0.02) elevated in RA compared with PsA (figure 3A,B, online supplemental figure S3). Interestingly, E1 endothelial cell cluster shows the highest expression of the VEGF receptors VEGFR1 (FLT1) and VEGFR2 (KDR) and high expression of NOTCH family members, specifically NOTCH4, NOTCH1 and their ligand, DLL4 (figure 3C). VEGF and NOTCH signalling result in fate decisions of endothelial cell specialisation towards stalk, tip or intermediate cell phenotypes that impact angiogenesis.18 In order to identify regulators of the E1 transcriptional profile in PsA and RA, we performed transcription factor (TF) usage estimation by analysing the expression of known, TF-regulated genes that are differentially expressed between PsA and RA. Interestingly, PsA and RA E1 endothelial cells show stark differences in TF usage, with TEA domain 1 (TEAD1) and myocyte enhancer factor 2A being the highest scored TF in RA E1 cells (figure 3D).19 Contrary to RA, PsA E1 cells show potential involvement of FOXP1 (figure 3D). The differences in endothelial cell TF usage, are potentially a reflection of differential transcriptional regulation, indicative of the distinct synovial angiogenesis in RA and PsA. Angiogenesis is the result of a highly regulated, orchestrated process, characterised by cell–cell interactions that define the fate and specialisation of endothelial cells.20 In order to examine potential cell interactions of endothelial cells belonging to cluster E1, differentially expressed receptors of cluster E1 were identified. Based on prior knowledge of receptor–ligand interaction potential, the heatmap of figure 3E depicts the top ligands for receptors expressed by endothelial cell cluster E1 (figure 3E). We then assessed the expression on all synovial cells of the top ligands for receptors of cluster E1 (figure 3F). Interestingly, endothelial cell cluster E1 shows potential for interaction with other endothelial cell clusters due to the high expression of several ligands by endothelial cells. While, limited, specific interactions between E1 cells and synovial fibroblasts and immune cells can be inferred from the extent of potential receptor ligand interactions (figure 3F).
Distinct endothelial cell profiles between patient with PsA and RA synovial biopsies. (A) Abundance of endothelial and pericyte cell clusters in patient with RA and PsA synovial biopsies. (B) Frequency of endothelial cell clusters in patient with PsA and RA synovial biopsies (n=4–5), data are presented as box and whiskers plots (min to max), symbols represent individual samples, statistical significance was determined by two-way analysis of variance with Sidak’s multiple comparisons test, *p=0.019. (C) Dotplot for the average scaled expression levels of angiopoietin receptor (TIE1 and TEK), VGF receptor (KDR and FLT1) and notch signalling elements (DLL4, NOTCH4 and NOTCH1). Dot size represents the percentage of cells per cluster expressing the indicated genes. (D) DoRothEA analysis of transcription factor usage by endothelial cell cluster E1 cells, based on expression of known downstream ligands. VIPER score difference between RA and PsA is shown. (E) Top 20 ligands with known and predicted interactions with receptors expressed by E1: endothelial cell cluster. (F) Dotplot depicting the potential sources of top ligands for cells of the E1 endothelial cell cluster. (G) Violin plots for the expression of VEGF receptor (KDR and FLT1), neuropilin-1 (NRP1), podocalyxin (PODXL), NOTCH4 and CD31 (PECAM1) by RA and PsA E1: endothelial cell cluster. DC, dentritic cells; NKT, natural killer cells; PsA, psoriatic arthritis; RA, rheumatoid arthritis.
The VEGF receptors, FLT1 and KDR and the VEGF-binding neuropilin-1, which modulates KDR expression, are upregulated in RA compared with PsA E1 endothelial cells21(figure 3G). Angiogenic NOTCH4 is upregulated by endothelial cells in response to VEGF, and previous histological analysis has revealed high NOTCH4 expression in the synovial tissue of patients with RA and PsA.22 Consistent with the upregulated VEGF receptor expression by RA E1 endothelial cells, NOTCH4 shows higher expression in RA compared with PsA (figure 3G). Platelet And Endothelial Cell Adhesion Molecule 1 (PECAM1), involved in endothelial cell adhesion and motility during angiogenesis, and podocalyxin, a key modulator of apical-basal endothelial cell polarisation and lumen formation, are also upregulated in RA compared with PsA E1 endothelial cells23 24 (figure 3G). The potential capacity of endothelial cells to interact with stromal and immune cells of the joint, and the identified transcriptomic differences between PsA and RA, indicate that the altered synovial blood vessel morphology between the two disease pathotypes is potentially the result of complex alterations in endothelial cell–cell crosstalk.
Identification of IL-1B expressing synovial macrophage cell cluster in RA and PsA synovial tissue
Macrophages are the most abundant immune cells of the synovial tissue with known protective, as well as pro-inflammatory, roles in RA disease pathogenesis.12 13 The abundance, calculated as the frequency of each cell cluster as part of all macrophages/DC per sample, of the three identified synovial tissue macrophage and one synovial DC clusters are comparable in PsA and RA (figure 4A,C). Interestingly, the macrophage cell cluster MC1, shows high level of IL-1B expression (figure 4B). Pathway enrichment analysis of the MC1 cluster of macrophages, shows enrichment of pathogenic signalling pathways including the tumor necrosis factor (TNF), IL-17, chemokine and cytokine–cytokine receptor pathways (figure 4D). PsA MC1 cluster macrophages use reduced myelocytomatosis proto-oncogene (MYC) compared with their RA counterparts (figure 4E). MYC can dictate the activation threshold for macrophages and early metabolic reprogramming by suppressing their response to lipopolysaccharide (LPS)-dependent stimulation, instead MYC has been shown to induce genes associated with an M2 transcriptional profile.25 26 To investigate which synovial cells have the highest potential to respond to the MC1 cluster macrophage-derived IL-1β, the expression of the IL-1β receptor, IL-1R1 was assessed. Synovial fibroblast clusters showed increased expression of IL-1R1 compared with other synovial cells, however, not all fibroblasts exhibited the same level of IL-1R1 expression (figure 4F). Interestingly, the F11 cluster, enriched in FAP+THY1+ synovial fibroblasts has higher expression of IL-1R1 compared with fibroblast cluster F1 (figure 4F).
Distinct macrophage cell transcriptomic profiles between patient with PsA and RA synovial biopsies. (A) Abundance of macrophage and DC cell clusters in patient with RA and PsA synovial biopsies. (B) Dot plot for the indicated markers in PsA (blue) and RA (green) macrophage and DC clusters. (C) Frequency of macrophage and DC cell clusters (calculated as a percentage of all macrophage and DC cells per sample) in patient with PsA and RA synovial biopsies (n=4–5), data are presented as box and whiskers plots (min to max), symbols represent individual samples, statistical significance was determined by two-way analysis of variance with Sidak’s multiple comparisons test, *p<0.05 were considered significant. (D) Term plot of the indicated pathways with significant enrichment in PsA compared with RA macrophage cluster MC1 following pathway enrichment analysis. Colour indicates up or downregulation of specific genes within the pathway and dot size represents significance. (E) Estimation of transcription factor activity by macrophages of cluster MC1 in RA compared with PsA. Transcription factor usage is estimated based on the differentially expressed genes with prior knowledge of genes regulated by specific transcription utilising bioinformatics package DoRothEA. Transcription factor enrichment score difference between RA and PsA is shown. (F) Violin plots for the normalised expression of IL-1R1 (IL-1B receptor) by all identified synovial cell clusters. IL, interleukin; PsA, psoriatic arthritis; RA, rheumatoid arthritis.
Limited in situ synovial T cell proliferation and high expression of TGFB1 in RA
We identified five clusters of synovial T cells, T1: CD8 T cells, T2: CD4 T cells, T3: CD74 T cells, T4: proliferating T cells and T5: NKT cells, in the synovial tissue of patients with PsA and RA. T cell cluster abundances, calculated as the frequency of each cell cluster as part of all T cells per sample, did not differ between PsA and RA except for a significant (*p=0.016) enrichment of T cell cluster T1: CD4 T cells in RA compared with PsA (figure 5A). T cells are key mediators of synovial inflammation, however, whether the primary mechanism of T cell accumulation in the synovial tissue is migration, or in situ synovial tissue expansion, remains poorly understood.27 MKI67, a marker of cell proliferation, was primarily detected in synovial T cell cluster T4 (figure 5B). Computational analysis of transcriptional profiles associated with cell proliferation based on the relative expression of 54 G2/M phase associated genes and 43 s phase associated genes, revealed that the T4 cluster was the only actively proliferating T cell cluster (figure 5C). T cell cluster T4, represents only 1.4%±0.6% and 8.8%±5.17% of synovial T cells in RA and PsA, respectively. Pathway enrichment analysis of synovial tissue T cells revealed differential enrichment of genes of the cytokine–cytokine receptor interaction and chemokine signalling pathways. Additionally, Tph/Tfh associated gene expression and indicative of tissue residency and early activation expression of CD69 are higher in RA compared with PsA synovial T cells (online supplemental figure S5)9. Demarcation of synovial T cell subsets and chemokine expression of PsA and RA T cells may be indicative of differential T cell involvement, however extensive further analysis and subclustering of T cells is required. Interestingly, transforming growth factor (TGF)B1 expression was increased in RA compared with PsA synovial T cells (figure 5D). Inhibition of TGF-β can limit synovial fibroblast hyperplasia in murine models of RA.28 29 With T cells being a critical source of TGF-β, TGFB1 expression was examined further. Patient with RA synovial tissue T cells of clusters T1, T3 and T5 exhibited high expression of TGFB1 compared with their PsA counterparts (figure 5E). TGF-β is a pleiotropic cytokine and determining its role in RA disease pathogenesis has been challenging. However, previous studies have identified that signalling pathways associated with TGF-β are enriched in RA but not osteoarthritis (OA) synovial fibroblasts and TGF-β1 messenger RNA expression correlates with patient with RA C-reactive protein (CRP) levels.30
Patient with RA synovial T cells express TGFB1. (A) Frequency of T cell clusters (calculated as a percentage of all T cells per sample) in patient with PsA and RA synovial biopsies (n=4–5), data are presented as box and whiskers plots (min to max), symbols represent individual samples, statistical significance was determined by two-way analysis of variance with Sidak’s multiple comparisons test, *p=0.016. (B) Violin plot for the log normalised expression of MKI67 by synovial T cell clusters. (C) Computational cell cycle analysis of synovial cells, based on the relative expression of 54 G2/M phase associated genes and 43 s phase associated genes, depicting cells in different stages of the cell cycle. (D) Term plot of the indicated pathways with significant enrichment in RA compared with PsA T cell clusters following pathway enrichment analysis. Colour indicates up or down regulation of specific genes within the corresponding pathway and dot size represents significance. (E) Violin plot for the log normalised expression of TGFB1 by RA and PsA synovial T cell clusters. HLA, human leukocyte antigen; IL, interleukin; PsA, psoriatic arthritis; RA, rheumatoid arthritis; TGF, transforming growth factor; UMAP, Uniform Manifold Approximation and Projection.
The majority of synovial tissue plasma B cells are potentially not derived from synovial tissue memory B cells
We identified four clusters of B cells, clusters B1 and B2 consisting primarily of memory B cells and clusters B3 and B4 consisting of plasma cells. Relative abundances expressed as the frequency of each cluster as a percentage of the total B cells for each sample did not differ between patient with PsA and RA synovial tissue (figure 6A). PCA plot of all enriched pathways per cluster following pathway enrichment analysis shows separation of B cells and plasma cells with plasma cell clusters grouping together while B cell clusters appear more dissimilar in the pathways that are being used (figure 6B). Despite the absence of any noticeable difference in the abundance of synovial B cell clusters between PsA and RA, specific pathways including the B cell receptor (BCR) signalling pathway were enriched in RA, compared with PsA, B cells (figure 6C). Ectopic lymphoid structure formation is a characteristic of aberrant RA synovial inflammation. It has been hypothesised that synovial plasma B cells emerge in the aforementioned structures as a result of in situ memory B cell differentiation.31 To evaluate this hypothesis, we examined κ and λ light chain usage by synovial tissue B cells. Due to allelic exclusion, a process that ensures B cells express one monospecific BCR following rearrangement of the light chains in early stages of B cell development, B cells express either a κ or a λ light chain.32 Synovial tissue B cells of cluster B1 and B2 and plasma cells of cluster B4 showed high expression of the κ light chain constant region (IGKC). Contrary to the majority of synovial plasma cells, cluster B3 demonstrates a clear preference for the λ light chain constant region (IGLC2) (figure 6D). Due to reports of a small population of B cells with dual BCR expression, the relationship between IGKC and IGLC2 expression was examined. Indeed, dual κ-expressing and λ-expressing synovial B cells were identified without data imputation. These B cells were primarily confined within the B2 B cell cluster (figure 6E). To assess the potential progression of synovial memory B cells to plasma cells, trajectory analysis was performed. Trajectory analysis uses gene expression to reconstruct the progression of cells along a lineage.33 Pseudotime, a measure of the distance of the cells from the starting point of the trajectory is used to infer the progression of the cells from the basal condition.34 The starting point of the trajectory was decided based on maximum pseudotime from B cells to plasma cells (figure 6F). Analysis of groups of co-regulated genes (modules) on the trajectory shows separation of plasma and B cell clusters (figure 6F). Interestingly, different gene modules achieve high scores between B cell cluster 1 and 2 (figure 6G). Plotting the dynamics of IGKC and IGLC2 expression as a function of pseudotime accentuates the separation between κ light chain-expressing B cells and λ light chain-expressing plasma cells (figure 6H). The distribution of IGKC and IGLC2 expression in relation to pseudotime invites the question of whether synovial B cells revise their BCR from κ to λ light chain, a phenomenon previously only observed in very early stages of B cell development.35 Therefore, the pseudotime was divided into segments and expression of the differentially expressed genes of pseudotime segment B (pseudotime distance 1 to 2) was evaluated as a function of pseudotime. Interestingly, differentially expressed genes of segment B showed high expression only in segment B (figure 6I).
Synovial plasma cells show biased usage of antibody λ over κ light chains. (A) Abundance of B cell and plasma cell clusters in patient with RA and PsA synovial biopsies, the frequency of each cluster was calculated as a percentage of all B cells and plasma cells per sample in PsA and RA synovial biopsies (n=4–5). Data are presented as box and whiskers plots (min to max), symbols represent individual samples, statistical significance was determined by two-way analysis of variance with Sidak’s multiple comparisons test, *p<0.05 were considered significant. (B) Principal component analysis plot of enriched pathways following pathway enrichment analysis of the identified B cell and plasma cell clusters, shows separation of B cells and plasma cells. (C) Term plot of the indicated pathways with significant enrichment in RA compared with PsA B cell clusters following pathway enrichment analysis. (D) Violin plot for the log normalised expression of IGKC (κ chain) and IGLC2 (λ chain) by the identified clusters. (E) Scatter plot of the relation between IGKC and IGLC2 expression of al B cell clusters. (F) Trajectory analysis of B cell and plasma cell clusters, arrows indicate starting point of pseudotime analysis. Due to the branching of the trajectory, in order to identify starting point of pseudotime analysis, analysis was initially performed with the indicated starting point. The highest pseudotime difference was then identified and used as the new starting point so trajectory analysis progresses from B cells to plasma cells. (G) Heatmap of co-regulated genes expressed per cluster as a function of pseudotime. Co-regulated genes were found using the find_gene_modules function in Monocle3 which runs UMAP on the genes rather than cells to group genes into modules using Louvain community analysis. (H) Expression of IGKC and IGLC2 as a function of pseudotime. (I) Differentially expressed genes of pseudotime fragment 1–2 expressed over the length of pseudotime. PsA, psoriatic arthritis; RA, rheumatoid arthritis; UMAP, Uniform Manifold Approximation and Projection.
Synovial T cell-derived TGF-β and macrophage IL-1β drive the transcriptome of proinflammatory synovial fibroblasts
As this study included unsorted synovial tissue single cell suspensions, it had the advantage of being able to examine potential networks of immune-stromal cell interaction involved in RA and PsA, thus reflecting the joint microenvironment. Synovial fibroblast clusters F1 and FAP+THY1+ F11, enriched in PsA and RA, respectively, were assigned the role of receiver cells to generate receptor–ligand interaction networks (figure 7A,B). Importantly, examination of the top receptor–ligand interactions (receptors expressed by fibroblasts; ligands expressed by all other synovial cells), indicate that the transcriptional profile of the proinflammatory fibroblast cluster F11 is potentially driven by synovial T cell derived TGFB and macrophage derived IL1B (figure 7B). Top 20 ligands with high receptor–ligand interaction potential with F11 fibroblast-expressed receptors include IL1B, TGFB, migration inhibitory factor (MIF), vascular cell adhesion molecule 1 (VCAM1) and NOTCH ligand Jagged1 (JAG1) (figure 7C). To assess the influence of IL1B and TGFB on the transcriptome of the fibroblasts of cluster F11 in RA and PsA, we used machine learning with random forest generation to evaluate to what extent IL1B and TGFB can predict the top per cent of differentially expressed genes of cluster F11 positioned downstream of the IL1B and TGFB receptors (figure 7D). IL1B but not TGFB could significantly (*p=0.028) predict the expression of downstream genes of fibroblast cluster F11 in PsA (figure 7D). Conversely, neither TGFB nor IL1B, could predict the expression of downstream genes of fibroblast cluster F11 in RA, however, the combination of both TGFB and IL1B shows high significance (***p<0.001) in predicting the downstream expression of differentially expressed genes of F11 fibroblasts in RA (figure 7D).
Identification of ligand receptor interactions that promote proinflammatory fibroblasts. (A) Circo plot depicting the top ligand and downstream target interaction for enriched in PsA synovial fibroblast cluster F1. (B) Circo plot depicting the top ligand and downstream target interaction for enriched in RA, proinflammatory synovial fibroblast cluster F11. (C) Heatmap of ligand receptor interactions for synovial fibroblast cluster F11. (D) Percentage of gene targets downstream of IL1B, TGFB1 or IL1B+TGFB as part of the top targets regulated by F11 fibroblast cluster receptors. Asterisks indicate significance of ligand-target interactions. (E) DoRothEA analysis of transcription factor usage by RA compared with PsA fibroblast cell cluster F11, based on expression of known downstream ligands, VIPER score difference between RA and PsA is shown. IL, interleukin; PsA, psoriatic arthritis; RA, rheumatoid arthritis; TGF, transforming growth factor.
Transcription factor usage analysis based on expression of known transcription factor-regulated genes, support a potentially increased usage of MYC and HIF1A by in RA F11 fibroblasts compared with PsA (figure 7E).
IL-1β and TGF-β synergistically drive metabolic adaptation of patient with RA synovial fibroblast and pro-inflammatory markers
The transcriptomically identified synergy of IL-1β and TGF-β in RA was assessed by in vitro characterisation of patient with RA synovial fibroblasts treated with combination of IL-1β and TGF-β. Previous studies have shown intercellular adhesion molecules 1 (ICAM-1) is upregulated in lining layer fibroblasts and facilitates tissue invasion and immune cell adhesion.36 Flow cytometric analysis of patient with RA synovial fibroblast showed no increase in ICAM-1 expression in response to TGF-β, however a significant increase in ICAM-1 (***p=0.0004) following treatment with IL-1β was observed (figure 8A). Importantly, the combined treatment with TGF-β and IL-1β resulted in a significant (**p=0.0047) increase in ICAM-1 expression compared with IL-1β only treated synovial fibroblasts (figure 8A). Fibroblasts are the main source of IL-6 in RA with pro-inflammatory synovial fibroblast secreting high levels of IL-6 in response to TNF-α.6 37 Similarly to the expression of ICAM-1, treatment of synovial fibroblasts with TGF-β did not lead to increased IL-6 secretion compared with untreated synovial fibroblasts, however the combined treatment with TGF-β and IL-1β resulted in a significant (**p=0.0018) increase in IL-6 compared with IL-1β only treated synovial fibroblasts (figure 8B). Dysregulation of synovial fibroblast metabolism with increased reliance on glycolysis has previously been associated with fibroblast pathogenic behaviour, interestingly, T cell derived soluble mediators have also been shown to enact metabolic switch of fibroblast towards glycolysis.38–40 Bioinformatic characterisation of gene modules based on oxidative phosphorylation or glycolysis reference pathways, deposited on Kyoto Encyclopedia of Genes and Genomes (KEGG) (pathways hsa00190 and hsa00010, respectively) showed on a transcriptional level, increased glycolysis and oxidative phosphorylation involvement in RA synovial fibroblasts of cluster F11 compared with PsA (figure 8C,D). To evaluate on a functional level that the potential synergistic effect of TGF-β and IL-1β regulate a pro-inflammatory synovial fibroblast we next used real-time Seahorse metabolic profiling to examine the extracellular acidification rate (ECAR) which measures glycolysis and the comparable oxygen consumption rate (OCR) which measured oxidative phosphorylation (figure 8E). While both TGF-β and IL-1β alone increased ECAR, the combined TGF-β and IL-1β treatment leads to a significant (*p<0.05) increase in all ECAR measurements (baseline glycolysis, glycolytic capacity, glycolytic reserve) compared with IL-1β only treated synovial fibroblasts (figure 8E). While there was no significant difference in the OCR profile in response to TGF-β and IL-1β stimulation alone, the combination of TGF-β and IL-1β, resulted in a decrease in maximal spare respiratory capacity compared with either alone. This led to an increase in the ECAR/OCR ratio, indicative of the cell’s reliance on glycolysis over oxidative phosphorylation, which is significantly higher in TGF-β and IL-1β treated synovial fibroblast compared with IL-1β or TGF-β alone (*p=0.029 and **p=0.0014, respectively) (figure 8F). The metabolic energy map demonstrated a shift in the overall metabolic profile of synovial fibroblasts where the combination of TGF-β and IL-1β resulted in a highly glycolytic synovial fibroblast phenotype (figure 8G). This shift in metabolic profile was paralleled by changes in synovial fibroblast mitochondrial function and morphology in response to TGF-β and IL-1β. Combined treatment with TGF-β and IL-1β resulted in significantly reduced tetramethylrhodamine methyl ester (TMRM) staining intensity compared with TGF-β or IL-1β (**p=0.0056, **p=0.0024, respectively), indicative of reduced mitochondrial membrane potential (figure 8H1). Decreased mitochondrial aspect ratio is indicative of reduced mitochondrial fusion, TGF-β and IL-1β treated synovial fibroblasts had significantly (*p=0.028) reduced aspect ratio compared with IL-1β only treated fibroblasts (figure 8H1) and aligns with the increased reliance of TGF-β and IL-1β treated fibroblasts to glycolysis over oxidative phosphorylation, as mitochondrial fusion supports oxidative phosphorylation.41
Effect of IL-1β and TGF-β on patient with RA synovial fibroblasts. (A) Flow cytometric analysis of CD54 (ICAM-1) expression by patient with RA synovial fibroblasts following treatment with IL-1β, TGF-β or a combination of both. Statistical significance was determined by one-way analysis of variance (ANOVA), symbols indicate individual samples (n=8), **p=0.0047, ***p=0.0007. (B) IL-6 secretion by patient with RA synovial fibroblasts treated under the conditions indicated. Statistical significance was determined by one-way ANOVA, symbols indicate individual samples (n=8), **p=0.0018, ***p=0.0008. (C) Violin plots of gene module expression scores in F11 fibroblasts, generated using the AddModuleScore function in Seurat, for oxidative phosphorylation-related and glycolysis-related genes derived from the Kegg gene sets ‘hsa00010’ and ‘hsa00190’. (D) Scatter plot of the relationship between oxidative phosphorylation and glycolysis gene module expression scores in the F11 fibroblast cluster. Scatter plots are coloured by density and Pearson correlation scores estimate the relationship between oxidative phosphorylation and glycolysis in RA-derived and PsA-derived F11 fibroblasts. (E) Seahorse bioenergetic analysis of patient with RA synovial fibroblast, ECAR and OCR measurements are shown under basal conditions or following treatment with TGF-β, IL-1 β or a combination of both. Statistical significance was determined with two-way ANOVA, p<0.05 was considered significant, n=10, points and lines represent mean values. (F) ECAR to OCR ratio of patient with RA fibroblasts for the indicated conditions, statistical significance was determined by one-way ANOVA, n=10, **p=0.0014, *p=0.019, ***p<0.0001. (G) Bioenergetic profile graph of patient with RA synovial fibroblasts under the indicated treatments, n=10. (H) Representative multiphoton microscope images of TMRM stained synovial fibroblast mitochondria under the indicated conditions. (I) Fluorescent intensity of TMRM and mitochondrial aspect ratio. Statistical significance was determined with one-way ANOVA, p<0.05 were considered significant. ECAR, extracellular acidification rate; IL, interleukin; OCR, oxygen consumption rate; PsA, psoriatic arthritis; RA, rheumatoid arthritis; TGF, transforming growth factor; TMRM,tetramethylrhodamine methyl ester perchlorate.
Discussion
Synovial inflammation in RA and PsA has a complex aetiology and is defined as the outcome of several underlying immunological mechanisms. Despite recent advances and increased availability of therapeutic options due to the introduction of biologics, patients often undergo exploratory treatments until they show an adequate response.42 For patients to experience sustained remission, achieving remission early is fundamental, therefore, lost time at initial stages of disease can have serious, lasting effects for the patients’ quality of life.43 Even when successful therapeutic intervention is achieved, long-term toxicity can have an impact on the patient.44 To advance towards precision medicine, it is crucial that we achieve a better understanding of the complex immune environment of the inflamed joint. The complexity of RA and PsA pathogenesis is confounded by multifaceted synovial and stromal cell interactions. Identifying and therapeutically targeting specific immune-stromal cell interactions has the potential to greatly reduce toxicity and improve therapeutic outcomes for both patients with RA and PsA. While significant advances have been made with the introduction of single-cell RNA sequencing (scRNAseq) and other omic approaches in RA, to our knowledge, this is the first time that a transcriptomic analysis of intact synovial single cell suspensions of the inflamed joint in RA and PsA has been performed,9 allowing for in-depth comparative cellular analysis of these two pathotypes.
Using high numbers of cells from intact synovial biopsy single cell suspensions for scRNAseq analysis offers distinct advantages. Previous studies have used sorted immune or stromal cells on the basis of CD45 expression or sorted specific populations, while prior knowledge of the cells included can expedite cluster analysis and annotation, it makes the generation of cell–cell interaction networks challenging.11 By not sorting synovial cells prior to RNAseq analysis, we remove an important potential source of variation between experimental data and the in situ environment of the joint. Additionally, cell frequencies of the populations analysed more faithfully mirror their relative abundances in the inflamed joint and allow for the generation of cell–cell interaction networks between immune and stromal cells.
The resulting transcriptomic analysis of the inflamed joint in RA and PsA revealed several previously unappreciated aspects of synovial inflammation. Limited T cell proliferation indicates that infiltrating T cells may be have a more important role in maintaining synovial T cells than previously anticipated. Additionally, differential light chain expression by synovial memory and plasma B cells leads to the hypothesis that in part, synovial plasma cells are recruited to the inflamed joint. Importantly, we identified differential abundance of synovial fibroblasts and their transcriptome in RA and PsA, alluding to disease specific mechanisms of synovial inflammation.
Recent studies have identified the existence of synovial fibroblasts with distinct functional characteristics in RA.14 Mizoguchi et al, have used RNAseq analysis of sorted synovial fibroblasts from two patient with RA and two patient with OA samples to identify three populations of synovial fibroblasts based on the expression of CD34 and THY1.6 Synovial fibroblasts negative for CD34 but expressing THY1 are expanded in RA and potentially contribute to synovial inflammation via the production of pro-inflammatory cytokines.6 Dividing synovial fibroblasts into functionally distinct subsets is an emerging field of study, an additional categorisation of synovial fibroblasts into distinct populations has been proposed where synovial fibroblasts are divided into two populations based on expression of FAP and THY1.14 FAP+THY1+ RA synovial fibroblasts express elevated levels of pro-inflammatory cytokines and chemokines including IL-6, chemokine (C-C motif) ligand 5 (CCL5) and CCL2, and are necessary to maintain synovial inflammation in a murine model of RA.14 Available information on distinct functions of synovial fibroblasts in PsA is scarce, however, recent studies show that PsA synovial fibroblasts can promote angiogenesis through regulation of endothelial cells.45 Previous studies have identified FAP+THY1+ synovial fibroblasts as pro-inflammatory with increased expression of C-C chemokine receptor type 2 (CCR2) and reduced expression of matrix metalloproteinase-3 (MMP3) compared with THY1– synovial fibroblasts, supporting our analysis of F11 compared with F1 synovial fibroblasts. FAP expression has been detected on synovial fibroblasts at early stages of inflammation in RA, indicating possible contribution of FAP+THY1+ fibroblasts early in disease pathogenesis, however, little is known regarding functionally distinct synovial fibroblast clusters in PsA.46 Herein, we report increased abundance of THY1– synovial fibroblast cluster F1 in PsA compared with RA, indicative of differences in fibroblast involvement in synovial inflammation between the two pathotypes.
Extensive angiogenesis is a characteristic of both RA and PsA, required to support the egress of immune cells and O2 from the periphery to the otherwise hypoxic environment of the inflamed joint.47 Despite the central role of angiogenesis in RA and PsA, morphological differences are evident with PsA synovial blood vessels presenting a tortuous, elongated and dilated phenotype, similar to that observed in tumour vasculature.45 Endothelial cell contribution to synovial inflammation is more complex and extends beyond pathogenic angiogenesis. Recent studies show a stromal crosstalk between synovial fibroblasts and endothelial cells in RA, with the latter providing NOTCH3-activating ligands to promote THY1-expressing synovial fibroblasts.20 In our analysis, five transcriptionally distinct synovial endothelial cell clusters were identified, however, only one cluster showed evidence of differential abundance between RA and PsA. Interestingly, VEGF receptor expression and NOTCH expression is higher in RA compared with PsA E1 endothelial cells. VEGF and NOTCH signalling cascades decide the fate of endothelial cell specialisation towards stalk, tip or intermediate cell phenotypes that impact angiogenesis. Differences in the angiogenic process in RA and PsA are additionally reflected by the differential TF usage with FOXP1, a TF that is required for neoangiogenesis and endothelial cell sprouting, upregulated in PsA compared with RA, while master regulator of endothelial cell metabolic reprograming during sprouting, TEAD1, shows enhanced usage in RA compared with PsA.19 48 49 Importantly, endothelial cells harbour high potential for interaction not only with distinct synovial fibroblast clusters but also with immune cells. Endothelial cell clusters show a plethora of potential interactions with other endothelial cell clusters; interactions that could be pivotal in the organisation of new blood vessels contributing to the pathogenesis of RA and PsA.
There is a great body of evidence regarding autoantibody involvement in RA disease pathogenesis.50 In addition to autoantibodies, novel functions of synovial B cells and synovial B cell populations have recently been described.4 51 52 The presence of ectopic lymphoid structures in RA has led to the popular hypothesis that synovial plasma cells are generated in the synovial tissue from clonally expanded, peripheral blood B cell infiltrates.53 However, no direct connection leading from synovial B cells to plasma cells has previously been described. B cells are monospecific and express a BCR consisting of two identical heavy and two identical light chains, the monospecificity of the B cell is ensured by the process of allelic exclusion.54 Following successful functional rearrangement of the heavy chain-encoding allele immunoglobulin heavy chain (IGH), the light chain-encoding loci are rearranged. Rearrangement of the light chain starts at the κ chain locus and, if no functional κ light chain emerges, recombination proceeds with the λ chain locus.55 Due to this process, the ratio of κ/λ chain usage by antibodies is biased towards κ light chains (κ/λ, 2:1).55 Interestingly, in patients with RA, anti–citrullinated protein antibody (ACPA)-expressing B cells show increased bias towards λ light chains.56 In agreement with this study, synovial tissue plasma cells show a clear preference for the expression of IGLC2 compared with IGKC. Surprisingly, most synovial B cells express IGKC, which indicates that synovial plasma cells are not derived from synovial B cells. It has to be noted that several rounds of BCR editing can result in a transition from κ light chain usage to λ light chain, this process however, has been reported at the very early stages of B cell development and there is no direct evidence to suggest that it can occur after the onset of somatic hypermutation.35 57 Despite the effectiveness of allelic exclusion, dual BCR-expressing B cells can emerge.58 59 Expression of two BCRs with different specificities could help autoreactive B cells evade central tolerance mechanisms. The data presented in this study suggest that synovial plasma cells are not derived from synovial B cells, therefore the role of ectopic lymphoid organs may be secondary to the differentiation of plasma cells. However, the presence of a small population of λ light chain expressing B cells or BCR editing in the synovial tissue could be contributing in the emergence of λ light chain expressing synovial plasma cells. Further studies are required to evaluate the origin of synovial plasma cells and assess their connection to plasma cells recruited from the periphery, and dual κ+λ+ light chain-expressing B cells.
The most abundant immune cells of the inflamed joint are synovial macrophages with distinct protective as well as pro-inflammatory roles in RA disease pathogenesis.13 60 61 In this analysis we have identified three macrophage and one DC cluster. Importantly we describe a macrophage population with high pathogenic capacity, characterised by high expression of IL1B. While present in similar abundances in both RA and PsA, IL1B-expressing macrophages are differentially regulated between RA and PsA, with increased nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) activation potential in PsA compared with RA. IL-1β, could be responsible for inducing pro-inflammatory programming by synovial fibroblasts, therefore, we assessed the expression of the IL-1β receptor by all synovial cells. Interestingly, highest degree of IL-1β receptor expression is observed by synovial fibroblasts and endothelial cells, however, not all fibroblast and endothelial cell clusters express the IL-1β receptor to a similar degree, thus, certain fibroblast and endothelial cell populations, such as fibroblast cluster F11, are more susceptible to programming by IL-1β.
Polyfunctional T cell responses with a bias towards Th17-like and Th1 have been reported in PsA and RA, respectively.8 62 Even within RA, synovial T cell cytokine contributions are not uniform and are indicative of discrete endotypes of disease.27 Further characterisation of T cell subsets and their potential for crosstalk with stromal cells in RA and PsA could significantly increase our understanding of T cell involvement in synovial inflammation. Specific RA T cell populations show high expression of the immunomodulatory TGF-B1 inhibition of TGF-β can limit synovial fibroblast hyperplasia in murine models of RA and, due to synergistic effects with other cytokines, it could be an attractive target for future therapeutic intervention.28 29
Despite current efforts, there is no consensus on the proportion and function of synovial fibroblasts in RA or PsA. An additional level of complexity is added when attempting to decipher the interactions that dictate the pro-inflammatory attributes of synovial fibroblasts. Recent studies have suggested a stromal crosstalk between endothelial cells and synovial fibroblasts as the driving force of the transcriptome of potentially pathogenic synovial fibroblasts.20 Indeed, decoding cell–cell interactions at the site of inflammation in RA and PsA can lead to the identification of novel avenues of targeted therapeutic intervention; therefore, instead of targeting systemic immunological pathways or entire immune populations, specific context-dependent cell–cell interactions can be disrupted leading to resolution of inflammation with minimal side effects for the patient. Herein we describe one such potential immune-stromal cell interaction: synovial T cell TGF-β and macrophage-derived IL-1β synergistically drive the transcriptome, cell adhesion molecule expression, pro-inflammatory cytokine secretion and metabolic profile of potentially pathogenic fibroblasts that are enriched in RA but not in PsA.
Sample heterogeneity may impact some of the findings, and while both patients with RA and PsA had comparable active moderate to high disease activity based on disease activity score-28 (DAS28) and disease activity index for psoriatic arthritis (DAPSA), respectively, and all biopsies were obtained from the same joint type, analysis of additional samples would allow for the assessment of transcriptomic profiles, disease status and response to treatment in RA and PsA. Additionally, while we have performed functional characterisation of synovial tissue fibroblasts following treatment with IL-1β and TGF-β, to validate the bioinformatically identified cellular crosstalk, further targeted in vitro studies will be required for the assessment of endothelial cell transcriptomic differences in RA and PsA and their potential impact on the characteristic vascular morphology of the two disease pathotypes. The balance between synovial T cell proliferation potential and homing from the periphery will require further flow cytometric analysis and proliferation assays. The herein presented scRNAseq and complementary functional assays importantly, highlight the need for the implementation of novel antibody multiplexing techniques using DNA barcoded antibodies allowing for target co-detection by indexing and spatial transcriptomic analysis for further characterisation of the proposed receptor–ligand interactions.
The first analysis of patient with intact RA and PsA synovial tissue single-cell suspensions is a significant step towards precision medicine and reveals previously unappreciated aspects of synovial inflammation. The potential reliance of the synovial T cell and plasma cell pools on renewal from the periphery and the identification of immune-stromal cell interactions can become a paradigm shift in the development of novel therapeutic options for inflammatory arthritis.
Methods
Patient sample collection and study approval
Patients with RA and PsA (defined by the American College of Rheumatology (ACR) and CASPAR Criteria, respectively) were recruited from the Rheumatology Department, St. Vincent’s University Hospital, UCD and Tallaght University Hospital, TCD. Patient with RA and PsA synovial tissue samples from knee joints with active inflammation were obtained under local anaesthetic using Wolf 2.7 mm needle arthroscopy or ultrasound guided biopsy as previously described, please see online supplemental table S1 for patient clinical characteristics.63 Patients with RA and PsA had comparable moderate to high disease activity (DAS28 4.6±1.1 and DAPSA 24.2±4.9, respectively) and biopsies had lymphocyte infiltrates and lining layer hyperplasia as scored by a clinical pathologist. The research was performed in accordance with the Declaration of Helsinki.
Synovial tissue sample preparation
Synovial tissue single cell suspensions were generated following enzymatic and mechanical digestion of the synovial biopsies as described previously.27 Briefly, approximately 15 synovial biopsies per patient were digested using the GentleMACS Tumor Dissociation Kit, human (Miltenyi Biotech) as per manufacturers’ instructions. Immediately after isolation, biopsies are washed with RPMI (Merck) before being placed in 4.7 mL RPMI supplemented with 200 µl of enzyme H, 100 µl enzyme R and 25 µl enzyme A in a gentleMACS C Tube followed by initial mechanical disruption of the tissue using programme h_tumor_01 on a gentleMACS Dissociator. Samples are enzymatically digested for a total of 1 hour at 37°C under continuous rotation using the MACSmix Tube Rotator with further applications of the gentleMACS Dissociator at the halfway point and at the end of the 1 hour incubation. The cell suspension is then passed through a 70 µM cell strainer. Viability of the cells is assessed with trypan blue exclusion staining and immediately cryopreserved in sterile filtered 10% dimethyl sulfoxide (DMSO)/fetal bovine serum (FBS) at a concentration of 1×106 cells per mL scRNAseq.
Frozen synovial biopsy cell suspensions were thawed quickly in a 37°C waterbath and transferred to sterile tubes with warm RPMI media (10% FBS). After washing and counting, a dead cell removal kit (Miltenyi cat#130-090-101) was implemented to increase viability. Using the Chromium Next GEM Single Cell 3’ Reagent Kits V.3.1 (10X Genomics), cells were loaded onto the GEM Chips. The 10X Genomics Chromium Next GEM Single Cell 3’ user manual was followed for all steps to generate complementary DNA (cDNA) libraries for each sample. cDNA quantifications and quality control were determined using the Agilent TapeStation. Final libraries were normalised, quantified (Illumina/ROX low, Kappa Biosystems), pooled based on 40 k reads/sample. Pooled libraries were sequenced on the Illumina NovaSeq using and S2 NovaSeq 6000 Reagents V.1 kit and a 100-cycle sequencing run.
ScRNA-seq data analysis
Initial processing
The gene expression raw sequencing data for the synovial tissue single cell suspensions were processed using Cell Ranger V.3.1 (10X Genomics, California, USA), with the 10X human transcriptome GRCh38.3.0.0 serving as a reference. Single-cell reads for each sample were converted to Seurat objects using the R package Seurat (V.4.0.3) in R (V.4.1). For each object representing an individual patient synovial tissue sample, data were filtered with genes detected in less than 3 cells, excluded from downstream analysis. Empty droplets were removed with function EmptyDrops (DropletUtils, V.1.12.1, code available here: https://github.com/MarioniLab/DropletUtils), followed by removal of cell doublets. Cell doublets were removed by using a newly described computational approach—DoubletFinder. DoubletFinder intersects transcriptional data and a data set specific artificial population of doublets, generated by averaging gene expression of randomly selected pairs of cells, in order to identify cell doublets (code available here: https://github.com/chris-mcginnis-ucsf/DoubletFinder).64 Apoptotic cells were removed by eliminating cells with a mitochondria associated gene expression of over 25%.9 One patient sample with a high frequency of cells over the mitochondrial gene expression threshold was completely excluded from further analysis due to potentially compromised/stressed live cells. (Data scaling) and normalisation were performed with the newly described sctransform (V.0.3.2) package. Sctransform uses non-heuristic approaches in order to scale the data based on Pearson residuals of a negative binomial regression, as such, it is superior to the widely used unique molecular identifier method since it is less susceptible to technical variations associated with widely different sequencing depths between deeply and shallowly sequenced cells of the same data set (code available here: https://github.com/satijalab/sctransform).65 This approach resulted in 178 804 cells from four patient with RA and five patient with PsA synovial tissue samples available for downstream cell clustering and analysis.
Clustering of major cell populations
PCA using the sctransform scaled data identified variable genes as input. Prior to clustering, integration of the Seurat objects representing synovial tissue samples from nine patients was performed with Harmony (V.0.1.0) (code available here: https://github.com/harmony-one/harmony).66 Harmony reduces variation associated with technical differences between samples that may otherwise, ‘mask’ biological differences, this is achieved by cell specific correction of the cell’s PCA coordinates. Clusters were identified with FindCluster function of Seurat and visualised on a UMAP plot (code available here: https://github.com/satijalab). Clustering efficiency was independently assessed by calculating the observed to expected edge weight ratios for all pairs of clusters (online supplemental figure S6). These ratios were calculated with function pairwiseModularity of package bluster of all off-diagonal cluster pairings had less observed to expected edge weight ratios compared with cells belonging to the same cluster. Further analysis of cluster stability was performed with bootstrapping the data in order to calculate the probability of a cell being randomly co-assigned to two clusters (function bootstrapStability, package bluster), clustering was efficient and stable. Cell identity was calculated with automated and prior knowledge approaches. SignleR and scCATCH were unsuccessful in annotating the majority of the cell clusters primarily due to a lack of synovial tissue specific data in reference data sets.15 Therefore, several clusters were annotated based on differential gene expression (DEG) profiles. DEGs were derived using the FindMarkers function of Seurat with Wilcoxon test and p values adjusted by Bonferroni correction (code available here: https://github.com/satijalab/seurat). DEGs were filtered on the basis of a minimum 0.25 of the cell cluster expression, a minimum log2 fold change of 0.5 and a p value below 0.05. Cell numbers per cluster per patient are included in online supplemental file 3.
Supplemental material
Trajectory analysis
Trajectory analysis of synovial tissue B cells was performed with Monocle3 (V.1.0). Co-regulated genes over pseudotime were identified by the find gene modules function. Additional analysis was performed by clustering cells in pseudotime fragments and then identifying DEG per pseudotime clustered cells (code available here: https://github.com/cole-trapnell-lab/monocle-release). DEG of specific pseudotime fragments were used as modules and their expression assessed over the complete pseudotime.
Cell cycle analysis
Cell cycle analysis was performed by scoring the relative expression of 54 G2/M phase associated genes and 43 s phase associated genes as per function CellCycleScoring of package Seurat. The base code used can be found as part of the cell cycle analysis vignette here: https://github.com/satijalab/seurat/blob/master/vignettes.
Cell–cell interaction analysis
Cell–cell interactions were identified with the nichenetr (V.1.0) package following the nichenetr vignette and code available here (https://github.com/saeyslab/nichenetr/tree/master/vignettes).67 Potential cell–cell interactions were identified based on gene expression and predetermined, based on prior-knowledge, receptor–ligand interaction pathways. One cell cluster was assigned the role of the ‘receiver’ population with its expression data intersected with known receptors and all other cells were assigned the role of ‘sender’ cells with their expression data intersected with known ligands. Receptors, downstream target genes of interest and ligands were based on DEG between a defined condition of interest and a reference condition. To evaluate to what extent TGF-β and IL-1β may regulate the differences between RA and PsA F11 fibroblasts a multi-ligand random forest model was used. This model uses the regulatory potential scores of TGF-β and IL-1β to predict the transcriptional programme of RA F11 fibroblasts and PsA F11 fibroblasts. The per cent of RA-specific or PsA-specific genes which belong to the 5% most strongly predicted targets were visualised and a one-sided Fisher’s exact test was used to test the significance of the association between the RA-specific and PsA-specific genes and whether they are part of the 5% most strongly predicted targets.
Pathway analysis
Pathway enrichment analysis was performed with pathfindR (V.1.6.2) active subnetwork analysis.68 DEG identified by FindMarkers were filtered based on log2 fold change and adjusted p value followed by run_pathfinR based on the KEGG database. Gene modules of oxidative phosphorylation-related and glycolysis-related genes were generated using the gauge package (https://bioconductor.org/packages/release/bioc/html/gage.html) to access the Kegg gene sets ‘hsa00010 Glycolysis/Gluconeogenesis’ and ‘hsa00190 Oxidative Phosphorylation’. The AddModuleScore function in Seurat calculated the average expression of each gene module in each disease group on a single cell level and subtracted the aggregate expression of 100 control genes. Gene modules were visualised by violin plot and scatterplot and Pearson correlation scores were calculated.
TF usage analysis
TF usage was estimated with package dorothea (V.1.4.1) with human regulons A, B, C.69 For visualisation, viper score differences between RA and PsA cell clusters were calculated.
Data visualisation
Plots are generated via ggplot2 (V.3.3.5), pheatmap (V.1.0.12), Seurat (V.4.0.3) and pathfindR (V.1.6.2) functions in R (V.4.0.1). Cell cluster abundance box and whisker plots (min to max) were generated in Prism based on relative frequency data; each symbol represents an individual sample.
Data availability and public access
All raw and processed files as well clinical information for each sample are deposited on national center for biotechnology information (NCBI), ascension number GSE200815 and are publically available without any restrictions of their subsequent use. Additionally, while detailed vignettes and base code is available on the bioinformatics platform GitHub (as indicated in methods) for all packages used in the analysis, if specific parts of the code are needed, they will become available on reasonable request.
Cellular bioenergetic function analysis
To examine the metabolic profile of IL-1B/TGF-B stimulated RA-fibroblast like synoviocytes (RAFLS), OCR and ECAR, reflecting oxidative phosphorylation and glycolysis, respectively, were measured using the Seahorse-XFe96 analyser (Seahorse Biosciences). RAFLS were seeded at 15 000 cells per well in a 94-well cell culture XFe microplate (Seahorse Biosciences) and allowed to adhere overnight. Following this, cells were then treated with either IL-1B (1 ng/mL), TGF-B (10 ng/mL), or a combination of IL-1B (1 ng/mL) and TGF-B (10 ng/mL) for 24 hours. Cells were then washed with assay medium (unbuffered DMEM supplemented with 10 mM glucose, pH-7.4) before incubation with assay medium for 1 hour at 37°C in a non-CO2 incubator. Basal oxidative phosphorylation/glycolysis was calculated by the average of three baseline OCR/ECAR measurements, respectively, obtained before injection of specific metabolic inhibitors; oligomycin (ATP-synthase-inhibitor), (2 µg/mL; Seahorse Biosciences, UK) trifluorocarbonylcyanide phenylhydrazone (FCCP) (mitochondrial uncoupler) (5 µM; Seahorse Biosciences) and antimycin A (complex-III inhibitor) (2 µM; Seahorse Biosciences). Oligomycin was injected to evaluate both the maximal glycolytic rate and ATP synthesis, determined by subtracting the amount of respiration left after oligomycin injection from baseline OCR. FCCP was injected to evaluate the maximal respiratory capacity (average of three measurements following injection) and respiratory reserve. Maximal respiratory capacity was determined by subtracting baseline OCR from FCCP-induced OCR and the respiratory reserve (baseline OCR subtracted from maximal respiratory capacity).
Cytokine measurements
To assess the effects of treatment with IL-1B and TGF-B singly and in combination on the production of pro-inflammatory mediators by the RAFLS, RAFLS were seeded in 24-well plates at a density of 5×105 per well and allowed to attach overnight. Cells were then incubated in 1% RPMI-1640 for 24 hours and subsequently stimulated with IL-1B (1 ng/mL), TGF-B (10 ng/mL), or a combination of IL-1B (1 ng/mL) and TGF-B (10 ng/mL). Supernatants were then harvested and levels of IL-6 (IL-6, R&D systems, UK,) were determined according to manufacturer’s conditions.
Flow cytometric analysis
Surface marker expression of RAFLS following stimulation with IL-1B (1 ng/mL), TGF-B (10 ng/mL), and IL-1B (1 ng/mL)+TGF B (10 ng/mL) was analysed by multiparameter flow cytometry. For extracellular staining, cells were seeded at 5×105 cells/well in a 24-well plate and stimulated with 1% complete Roswell Park Memorial institute medium (cRPMI) supplemented with the specific cytokines IL-1B and TGF-B prior to staining. For the gating strategy, the cells were initially gated based on forward and side scatter and doublets were removed. LIVE/DEAD fixable NIR (Thermo Fisher) viability dye was used to eliminate dead cells. To eliminate non‐specific binding of monoclonal antibodies to the Fc‐γ receptor (FcγR), samples were blocked with a human FcγR‐binding inhibitor (TruStain FcX Receptor blocking solution (BioLegend)) prior to antibody staining. The following antibodies were used in combination to investigate surface markers expressed by the stimulate RAFLS: Podoplanin FITC (Clone NC-08) (BioLegend), Human FAP Alexa Fluor 700 (Clone 427819) (RnD), CD90 Brilliant Violet 421 (Clone 5E10) (BioLegend), CD34 Brilliant Violet 510 (Clone 581) (BioLegend), CD54 Brilliant Violet 605 (Clone HA58) (BD), CD45 Brilliant Violet 650 (Clone HI30) (BioLegend), CD146 Brilliant Violet 711 (Clone P1H12) (BioLegend), HLA-DR Brilliant Violet 785 (Clone L243) (BioLegend), FAS-L PE (Clone NOK-1) (BioLegend) and CD309 PE/Cy7 (Clone 7D4-6) (BioLegend). Samples were acquired using the LSR Fortessa Flow Cytometer (BD) and analysed using FlowJo (V.10) software. Fluorescent minus one gating controls used were appropriate.
Tetramethylrhodamine methyl ester staining and analysis
Mitochondria Imaging was performed using a custom upright (Olympus BX61WI) laser multiphoton microscopy system equipped with a pulsed (80 MHz) titanium: sapphire laser (Chameleon Ultra, Coherent, USA), water-immersion 25× objective (Olympus, 1.05NA) and temperature controlled stage at 37°C. Fibroblasts were seeded in 35 mm petri-dishes and stained at 37°C for 30 min with 250 nM of tetramethylrhodamine methyl ester and then washed with phosphate-buffered saline. Two-photon excitation was performed at 850 nm and fluorescence emission was collected at 580–638 nm. Fluorescence images were acquired and quantified for their intensity and mitochondria morphology using CellProfiler with a custom built project pipeline, script available on request.70 71
Patient and public involvement
In this study we analysed specific immune and stromal cells obtained from patients with RA and PsA. While no patients were involved in setting the research question, the outcome measures or recruitment plans for the study, the group have hosted a number of patient information evenings where we have described the research, current project and the importance of patient engagement. No patients were asked to advice on interpretation or writing up of results, however, in collaboration with patient partners we developed a series of patient partnership workshops where lay dissemination of the study research to relevant patient groups was performed, with patient feedback now incorporated in future studies.
Data availability statement
Data are available in a public, open access repository. Data are available upon reasonable request. All raw and processed files as well clinical information for each sample are deposited on national center for biotechnology information (NCBI), ascension number GSE200815 and are publically available without any restrictions of their subsequent use. Additionally, while detailed vignettes and base code is available on the bioinformatics platform GitHubgithub (as indicated in methods) for all packages used in the analysis, if specific parts of the code are needed, they will become available upon reasonable request.
Ethics statements
Patient consent for publication
Ethics approval
The study was approved by the Institutional Ethics Committees of St Vincent’s University Hospital UCD and Tallaght University Hospital, TCD. Participants gave informed consent to participate in the study before taking part.
Acknowledgments
We would like to thank all the patients who consented to be involved in this study.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
DJV and UF are joint senior authors.
Handling editor Josef S Smolen
Twitter @Sarah_M_Wade, @mary_canavan, @drmgmonaghan
Contributors Guarantor, AF,UF, Conceptualisation, AF, UF, DJV; methodology, AF, CMS, OT; software, AF, CMS, VK; validation, AF, CMS, OT; formal analysis, AF, UF, DJV, SN, CMS, OT, NN; investigation, AF, UF, DJV; resources, UF, DJV; data curation, AF, UF, DJV, VK, MC, SN, OT, NN, MGM, VM, JMF, SC, MH, CC, SMW, L-YH; writing—original draft preparation, AF, CMS, UF, DJV; writing—review and editing, AF, CMS, OT, UF, DJV; visualisation, AF, UF, DJV, CMS, NN; supervision, UF, DJV; project administration, UF, DJV; funding acquisition, UF, DJV.
Funding This research was funded by Health Research Board of Ireland, grant number ILP-POR-2017–047; Centre for Arthritis and Rheumatic Diseases, CARD-2019–01; and Arthritis Ireland.
Competing interests None declared.
Patient and public involvement Patients and/or the public were not involved in the design, or conduct, or reporting, or dissemination plans of this research.
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.