Notch signalling drives synovial fibroblast identity and arthritis pathology.

The synovium is a mesenchymal tissue composed mainly of fibroblasts, with a lining and sublining that surround the joints. In rheumatoid arthritis the synovial tissue undergoes marked hyperplasia, becomes inflamed and invasive, and destroys the joint1,2. It has recently been shown that a subset of fibroblasts in the sublining undergoes a major expansion in rheumatoid arthritis that is linked to disease activity3-5; however, the molecular mechanism by which these fibroblasts differentiate and expand is unknown. Here we identify a critical role for NOTCH3 signalling in the differentiation of perivascular and sublining fibroblasts that express CD90 (encoded by THY1). Using single-cell RNA sequencing and synovial tissue organoids, we found that NOTCH3 signalling drives both transcriptional and spatial gradients-emanating from vascular endothelial cells outwards-in fibroblasts. In active rheumatoid arthritis, NOTCH3 and Notch target genes are markedly upregulated in synovial fibroblasts. In mice, the genetic deletion of Notch3 or the blockade of NOTCH3 signalling attenuates inflammation and prevents joint damage in inflammatory arthritis. Our results indicate that synovial fibroblasts exhibit a positional identity that is regulated by endothelium-derived Notch signalling, and that this stromal crosstalk pathway underlies inflammation and pathology in inflammatory arthritis.

We examined the transcriptional heterogeneity of synovial stromal cells by fluorescence-activated cell sorting for CD45- cells followed by droplet-based single cell RNA-sequencing (scRNAseq). After stringent quality control, we obtained data on 35,153 high quality stromal cells that we corrected for donor specific effects with Harmony[6] (Extended Data Fig. 1a). We used graph-based clustering[7], differential expression analysis, and canonical markers to identify 4 major stromal subpopulations (Fig. 1a): VWF+ vascular endothelial cells, PRG4(Lubricin)+ fibroblasts, THY1+(CD90) fibroblasts, and MCAM(CD146)+ mural cells (Extended Data Fig. 1b, Supplementary Table 1). Based on previous studies[4,8], we designated PRG4+ cells as lining and THY1+ cells as sublining fibroblasts (Extended Data Fig. 1c). Using flow cytometric analysis, we confirmed that lining and sublining fibroblasts, mural cells, and vascular endothelial cells are the predominant stromal cell types in synovia (Extended Data Fig. 1e). In flow cytometry data, lining fibroblasts represent the majority (48%) of stromal cells in OA, while the RA synovium is characterized by the expansion of sublining fibroblasts (52% vs 31%, p=0.002) (Extended Data Fig. 1f). However, in our scRNAseq dataset, lining and sublining fibroblasts did not separate into entirely disjoint clusters. Instead, the two populations were connected, suggesting the presence of intermediate subpopulations (Fig. 1a). This is exemplified by the continuous, non-mutually-exclusive distribution of the canonical sublining and lining fibroblast markers THY1[9] and PRG4[10] in UMAP space (Extended Data Fig. 1c).

Extended Data Figure 1.

Characterization of synovial stromal cells by single cell RNAseq and flow cytometry.

a, UMAP projection of single cell RNA-seq data of 35,153 stromal synovial cells from 12 donors before Harmony integration (left) and after Harmony integration (right). Each cell is coloured by donor source. b, Mean expression (colour) and percentage (size) of stromal markers among lining fibroblasts, sublining fibroblasts, mural cells, and endothelial cells. c, Expression of lining marker PRG4 and sublining marker THY1 in UMAP projection. d, UMAP embedding of fibroblasts and mural cells coloured by their positional identity (grey to blue). e, Representative flow cytometric gating scheme for analysis of synovial stromal cells. Intact cells are identified based on forward scatter (FSC-A) and side scatter (SSC-A) characteristics. Dead cells (PI+), red blood cells (CD235a+) and leukocytes (CD45+) are excluded from analysis. Endothelial cells (CD31+ CD146+) and mural cells (CD31- CD146+) can be distinguished from synovial fibroblasts based on CD146. Within CD31-CD146- gate, lining fibroblasts (PDPN+CD90-) and sublining fibroblasts (PDPN+ CD90+) can be identified. Red window indicates intermediates between cell types. f, Flow cytometric quantification of stromal cell populations in RA (n = 9) and OA (n = 11) synovia. Shown are mean percentage of endothelial cells (red), mural cells (blue), lining fibroblasts (tan), and sublining fibroblasts (green) of total stromal (CD45-) cells. g, Density plot of fibroblast composition in RA (orange) and OA (blue) synovia (n = 6) along positional axis. h, Archetypal analysis assigned each cell a probability distribution over the 6 biologically defined archetypes as well as the archetype representing ambient mRNA (background). Confusion matrix represents the probability that a cell with known type (y-axis) is assigned to one of 7 archetypes (x-axis). Each row is normalized to sum to 1. i, The positional identity (x-axis) of each fibroblast plotted against the probability that the cell was classified to the ambient RNA archetype (y-axis). j, heatmap visualization of transcriptional gradients along positional axis. k, pathway enrichment of 71 genes in the Intermediate group in (j), using GO biological process terms.

Figure 1.

Single-cell RNA-seq reveals fibroblast positional identity.

a, UMAP projection of 35,153 synovial cells from RA (n = 6) and OA (n = 6) patients. b, Trajectory analysis of synovial fibroblasts. Black line represents the trajectory and fibroblasts are coloured by their position from sublining (blue, position 0) to lining (grey, position 100). c, Expression of lining markers PRG4 (top) and sublining marker THY1 (bottom) along positional axis. d-f, Representative spatial analysis of synovial fibroblasts (n = 5 RA, n = 5 OA). d, Immunofluorescence microscopy showing VWF+ endothelial cells, CD146+ mural cells, PRG4+ lining and CD90 (THY1)+ sublining fibroblasts. e, Cells from d were segmented and abstracted spatially by their centroids (Methods). Endothelial cells were marked red and fibroblasts were colored by their CD90:PRG4 ratio (grey (low) to blue (high)). f, Fibroblast CD90:PRG4 ratio as a function of distance to the nearest endothelial cell. Cell aggregates were colored by their mean CD146 expression. In c and f, individual cells were binned by frequency (n = 100) along the x-axis and summarized by their mean (circle) and standard deviation (line). Significance determined by spearman correlation. See Supplementary Data for all images and spatial analysis performed.

Fibroblasts exhibit positional identity

Since the synovial lining and sublining are anatomically defined, we hypothesized that this gradual transition in fibroblast marker expression might capture a broader transcriptional gradient that encodes fine-grained positional identity. We performed an unbiased trajectory analysis[11] which models the linear transcriptional axis between sublining fibroblasts on one end, lining fibroblasts on the other, and transition states in the middle (Fig. 1b). We assigned each fibroblast a trajectory score along this axis, ranging from 0 at the sublining to 100 at the lining end (Fig. 1b). Along this trajectory axis, canonical markers PRG4 and THY1 show a gradual change in expression (Fig. 1c). Consistent with the expansion of THY1+ sublining fibroblasts in RA, we observed more cells from RA donors at the sublining pole of the gradient (Extended Data Fig. 1g). Using an archetype-based analysis, we confirmed that intermediate fibroblast states were not the result of technical artifact arising from mRNA contamination (Extended Data Fig. 1h, i). Moreover, we identified a group of 71 functionally related genes specific to the intermediate states (Extended Data Fig. 1j, k, Supplementary Table 2).

We next considered whether mural cells, a mesenchymal cell population positioned around vascular endothelial cells, could be consistently placed on this trajectory axis. Using genes associated with positional axis, mural cells consistently grouped towards the sublining end (score=20+−11) of the axis (Extended Data Fig. 1d). This placement suggests an anatomical endpoint for the sublining portion of the gradient adjacent to the endothelium. We thus hypothesize that the transcriptional gradient might correspond to a spatial, positional axis, spanning perivascular-positioned cells on one end to lining layer-positioned cells on the other.

To test this hypothesis directly, we used confocal microscopy to visualize the anatomical distribution of cells expressing positional markers in synovial tissues, including VWF, CD146(MCAM), CD90(THY1) and Proteoglycan4(PRG4) (Fig. 1d). We used immunofluorescence staining to label VWF+ endothelial cells and compute the CD90:PRG4 ratio in fibroblasts (Fig. 1e, Supplementary Table 3). In all images, we found a significant inverse correlation between distance and marker ratio (mean r=−0.33 +/− 0.12 SD, p<6.6×10−14, Fig. 1f, Supplementary data). Cells within 10–30 μm of an endothelial cell had the highest CD90:PRG4 expression, while cells further than 100 μm had the lowest (Fig. 1f). Consistent with the positional hypothesis, the CD90:PRG4 ratio decreased gradually with distance away from endothelial cells. As expected, we observed high expression of mural cell marker CD146 in cells nearest to endothelial cells (Fig. 1f). We repeated and confirmed this spatial trend using 3 additional positional markers identified in our trajectory analysis: sublining enriched Gamma-glutamyltransferase-5 (GGT5), lining enriched CD55 and Podoplanin (PDPN) (Extended Data Fig. 2, Supplementary Tables 4–6, Supplementary data). Together, single cell trajectory transcriptional analysis and imaging analysis support the existence of a smooth positional gradient between perivascular and lining fibroblasts.

Extended Data Figure 2.

Spatial analysis of synovial fibroblast positional markers.

a, Normalized expression (log counts-per-ten-thousand) of position-associated fibroblast markers CD55 (left), GGT5 (middle) and PDPN (right) along positional axis. b, Representative microscopic images of synovial tissues where CD55 (yellow), CD90(THY1) (green), CD146(blue), and VWF (red) are visualized by immunofluorescence staining. c, Cells labeled as endothelial cells were coloured in red and fibroblasts were coloured from low CD90:CD55 ratio (grey) to high CD90:CD55 ratio (blue). d, Spatial correlation between fibroblast CD90:CD55 ratio and distance from nearest endothelial cell. e, Representative microscopic image of synovial tissues whwere Podoplanin(PDPN) (red), Gamma-glutamyltransferase-5(GGT5) (green), and CD31 (blue) are visualized by immunofluorescence staining. f, Cells labeled as endothelial cells were coloured in red and fibroblasts were coloured from low GGT5:PDPN (grey) ratio to high GGT5:PDPN(blue) ratio. g, Spatial correlation between fibroblast GGT5:PDPN ratio and distance from nearest endothelial cells. For b-d, n = 4 RA and n= 5 OA synovial tissues were analyzed. For e-g, n = 7 RA synovial tissues were analyzed. See Supplementary data for individual images and spatial analysis. For d and g, Spearman correlation values and significance were computed with the base R cor.test function. Cells were binned by frequency into groups of 50 cells along the x-axis and summarized by their mean (dot) and standard deviation (line).

Endothelial cells provide positional cues for fibroblasts

Positional identity may represent a stable, cell-intrinsic program that reflects different developmental lineages[12] or anatomical sites[13]. Alternatively, positional identity may be determined by instructive signals derived from the local microenvironment which may be lost ex vivo. We tested if fibroblast positional identity is a stable feature by examining gene expression of CD90(THY1)+ and CD90(THY1)- fibroblasts over serial passages. By projecting transcriptomic profiles of fresh and in vitro passaged fibroblasts onto the position gradient trajectory[14], we found that after serial passages ex vivo, transcriptional profiles of CD90(THY1)+ and CD90(THY1)- fibroblasts converge, suggesting a loss of positional identity (Fig. 2a). Whereas freshly isolated CD90(THY1)+ fibroblasts mapped to the perivascular pole (position=28+−14) and CD90(THY1)- fibroblasts to the lining pole (position=73+−11) of the positional axis, passaged fibroblasts lost their positional identities and converged in between the perivascular and lining poles (Extended data Fig. 3a). This stark retreat from the poles reflects a loss of positional identity ex vivo after only two passages.

Figure 2.

Endothelial cells establish sublining fibroblast positional identity.

a, RNAseq profiles of CD90(THY1)- (blue) and CD90(THY1)+ (red) fibroblasts at indicated passage number (Fresh: n = 7 and n = 15, passage 1: n = 7 and n = 16, passage 2: n = 4 and n = 8) projected onto the scRNAseq trajectory embedding of fresh fibroblasts (black line). b, Apical view of a fibroblast organoid (top) and a fibroblast plus endothelial cell organoid (bottom). c, Representative hematoxylin and eosin staining of synovial organoids (n =4) demonstrating synovial lining (top) and sublining (bottom). d, Confocal microscopy images of organoids (n = 3) where fibroblasts were stained with PKH67 (green) and endothelial cells stained with PKH26 (red). Arrow indicate endothelial tubules with surrounding fibroblasts. e, Projection of scRNAseq profiles of fibroblast organoids (top, 4,336 cells) and fibroblasts plus endothelial cell organoids (bottom, 2,076 cells) onto the UMAP embedding of synovial tissue cells. Cells derived from organoids were coloured in orange and cells from synovial tissue (n = 35,153 cells) were coloured in grey.

Extended Data Figure 3.

Endothelial cells mediate differentiation of CD90 (THY1)+ fibroblasts.

a, Fibroblast positional identity scores of CD90(THY1)+ (red) and CD90(THY1)- (blue) fibroblasts at indicated passage number (Fresh: n = 7 and n = 15, passage 1: n = 7 and n = 16, passage 2: n = 4 and n = 8). Boxplots summarize the median, interquartile range, and 95% quantiles of the positional values. b, Positional identity of fibroblasts from the fibroblast-only (top, n = 4,336 cells) and fibroblast plus endothelial cell (Bottom, n = 2,076 cells) organoids. Cells are colored by inferred position along positional axis from 0 (perivascular pole) to 100 (lining pole). Perivascular fibroblasts (red) were defined as cells in position less than 20, while intermediate fibroblasts were defined as position between 20 and 80. c, Flow cytometric quantification of synovial CD90(THY1)+ fibroblasts and CD31+ endothelial cells from 22 synovial tissues. d, Percent of CD90(THY1)+ fibroblasts based on Doppler ultrasound scores of RA patients (n = 20) from the AMP-RA/SLE consortium. e, THY1 expression measured by RT-qPCR in fibroblasts after stimulation with indicated recombinant protein (n = 3 replicates for Wnt3a and Wnt5a, n = 4 replicates for other conditions). THY1 expression shown as fold-change over vehicle control for each condition. Significance determined by one-sample t-test. f-h, scRNAseq analysis of synovial organoids (n = 22,164 cells from 3 technical replicates). f, Normalized expression of defining marker gene for cell types: MKI67 for proliferating fibroblasts, VWF for endothelial cells, and THY1 for the THY1high and THY1low fibroblast groups. g, Representation of 4 major cell types, by proportion of cells in each organoid condition. h, NOTCH activation score in fibroblast derived from organoids, separated by organoid condition. ECs = endothelial cells. DAPT = organoids cultured in the presence of 10uM γ-secretase inhibitor DAPT. Significance determined by spearman correlation (c, d) and one-sample t-test (e).

Given the close anatomical link between CD90(THY1)+ fibroblasts and endothelium, we tested if culturing fibroblasts with endothelial cells could regenerate mural cell and/or sublining fibroblast positional identity in vitro. We developed a free-floating, 3-dimensional synovial tissue organoid system that supports stable co-culture of synovial fibroblasts and endothelial cells (Fig. 2b). Histologically, fibroblasts plus endothelial cell organoids formed two discrete zones: a fibroblast lining and sublining enriched with endothelial tubules (Fig. 2c). Confocal microscopy with fluorescent cell-labeling revealed a close association between endothelial tubules and fibroblasts in the mixed cell organoids (Fig. 2d), strikingly similar to vascularized sublining in RA (Fig. 1d). Next, we compared fibroblasts from mixed cell organoids to those in fibroblasts-only organoids by assaying gene expression with scRNAseq. We used Harmony[6] and k-nearest-neighbor based projection to map organoid cells to corresponding synovial tissue cells (Fig. 2e). Most cells (99.7%) from the fibroblast-only organoid were restricted to an intermediate region of the tissue fibroblasts reflecting neither strong sublining nor lining gene expression characteristics (Fig. 2e, Extended Data Fig. 3b). In contrast, 5.9% of fibroblasts from mixed fibroblasts plus endothelial cell organoids localized to the perivascular pole while 15.4% clustered with mural cells (Fig 2e, Extended Data Fig. 3b). This suggests that endothelial cell-derived signals are sufficient to differentiate fibroblasts into both perivascular fibroblasts and mural cells ex vivo. Consistent with endothelial cells as a source of positional cues for CD90(THY1)+ fibroblasts, we observed significant correlations between the number of CD90(THY1)+ fibroblasts and the degree of synovial tissue vascularization by flow cytometry (r=0.80, p=8.06×10−6) and joint doppler ultrasound (r=0.70, p=5.26×10−4) (Extended Data Fig. 3c, d).

NOTCH signaling controls fibroblast positional identity

Next, we sought to identify putative endothelial cell-derived signals that could confer perivascular fibroblast positional identity via an inductive programme. With ligand-receptor analysis in the synovial tissue and tissue organoid scRNAseq datasets (Supplementary Table 7), we prioritized NOTCH signaling (NOTCH1 and NOTCH3 receptors and JAG1, JAG2, and DLL4 ligands) as potential pathways for fibroblast-endothelial signaling (Fig. 3a). In vitro, stimulation with recombinant JAG1 and DLL4 differentiated fibroblasts towards a CD90(THY1)+ sublining phenotype (JAG1: 4.33 +/− 0.75 fold change, p=1.4×10−3. DLL4: 11.05 +/− 0.51 fold change, p=2.8×10−5, Extended Data Fig. 3e). In synovial tissue organoids, we observed a THY1high fibroblast cluster in fibroblast plus endothelial cell organoids, but not in fibroblast-only organoids (Fig. 3b, Extended Data Fig. 3f, g). Using a fibroblast-specific NOTCH gene signature defined by ligand activation (Methods), we found that THY1high fibroblasts were enriched in the NOTCH activation signature, compared to THY1 cells from the fibroblast-only organoid (Fig. 3c, Extended Data Fig. 3h). Consistent with a critical role of NOTCH signaling during differentiation of THY1high fibroblasts, addition of NOTCH inhibitor DAPT abrogated NOTCH activation and blocked THY1high fibroblast differentiation (Fig. 3c, Extended Data Fig. 3g, h).

Figure 3.

NOTCH signaling drives CD90/THY1+ fibroblast differentiation.

a, Ligand-receptor analysis in synovial tissue (top) and organoid (bottom) scRNAseq datasets. Black lines indicate highly expressed ligands in endothelial cells (cyan) and receptors in fibroblasts (orange). Ligand-receptor pairs identified in both datasets were highlighted in red. b-c scRNAseq analysis of synovial organoids (n = 22,164 cells from 3 replicates). b, UMAP projection of organoid cells where each cell is coloured by culture condition: fibroblast only (blue), fibroblast plus endothelial cells (ECs, orange), fibroblasts plus endothelial cells treated with 10uM DAPT (red). c, NOTCH activation score (grey (low) to red (high)) in organoid cells. Each organoid condition was projected separately. Circle = THY1high fibroblasts. d, Immunohistochemistry staining of NOTCH3 (purple) and Elastin (black) in synovial tissue (n = 5 RA and 5 OA). A = arterial endothelium, V = venous endothelium. e, NOTCH activation score of synovial tissue fibroblasts. Individual cells were binned by frequency (n =100) along the x-axis and summarized by their mean (circle) and standard deviation (line). f, Percentage of NOTCH3+ fibroblasts in OA (n = 5) and RA (n = 10) synovia determined by flow cytometry. g, Percentage of NOTCH-activated fibroblasts in OA and RA (n = 6). h, NOTCH activation scores from bulk RNAseq profiles of RA (n = 14) and OA (n = 12) fibroblasts. Significance determined by spearman correlation (e), two-tailed t-test (f, g), and two-tailed Wilcoxon rank sum test (h). All boxplots summarize the median, interquartile range, and 95% quantile range.

Since NOTCH3 is the Notch receptor more selectively expressed in synovial mural cells and fibroblasts (Extended Data Fig. 4a), we hypothesized that NOTCH3 was important in mediating Notch signaling in synovial fibroblasts. In synovial tissue, NOTCH3 mRNA and protein expression was highest in mural cells, including pericytes and vascular smooth muscle cells, followed by sublining fibroblasts (Extended Data Fig. 4b–d). Consistent with NOTCH3 receptor activation in mural cells and sublining fibroblasts, we observed increased NOTCH3 intracellular domain (NICD3)[15] in these cells as well as in vitro, JAG1- or DLL4- stimulated fibroblasts (Extended Data Fig. 4e, f),

Extended Data Figure 4.

NOTCH3 signaling in synovial fibroblasts and mural cells.

a-b, scRNAseq analysis in NOTCH1 and NOTCH3 expression in synovial stromal cells. a, Expression of NOTCH1 (top) and NOTCH3 (bottom) in lining fibroblasts (orange), mural cells (red), sublining fibroblasts (green), and endothelial cells (blue) from human primary synovial tissue scRNAseq data set (n = 35,153 cells from 6 RA and 6 OA synovial tissues). b, Mural cells were subdivided into 2 major mural subpopulations: vascular smooth muscle cells (orange) and pericytes (blue). No difference in expression of NOTCH3 observed (Wilcoxon rank sum test p=0.86) between pericytes and vascular smooth muscle cells. c, Representative synovial tissue sections (n = 6) showing RNAscope staining for NOTCH3 in purple, immunohistochemistry staining (IHC) for CD90 in brown and nuclei in blue. Arteriole marked by ∆. d, Summary of flow cytometric analysis of NOTCH3 mean fluorescence intensity (MFI) on synovial mural cells (n = 18), sublining fibroblasts (n = 40), lining fibroblasts (n = 20), endothelial cells (n = 15) and leukocytes (n = 19). Boxplot summarize the median, interquartile range, and 95% quantile range. Significance determined by spearman correlation. e, Representative immunohistochemistry staining of NOTCH3-intracellular domain (NOTCH3-ICD) in RA synovial tissue (n = 3). A = arterial endothelium. V = venous endothelium. f, Immunoblot of synovial fibroblast lysates probed with antibody against NOTCH3 intracellular domain (NICD3) (top) or GAPDH (bottom). Treatments: 10mmol/EDTA (Notch activation), plate-coated DLL4 (5ug/ml), or plate-coated JAG1 (5ug/ml), in the presence or absence of 10uM DAPT (γ-secretase inhibitor). n = 3 independent experiments. For gel source data, see Supplementary Data 1.

Given that Notch ligands JAG1 and DLL4 expression was higher (auROC>0.58) in synovial arterial (PODXL+) than venous (DARC+) endothelial cells (Extended Data 5a, b), we suggest that the source of fibroblast NOTCH3 signaling was derived from arterial endothelium. Consistent with this, NOTCH3 was detected on mural cells and fibroblasts surrounding arterial endothelium (Fig. 3d), and its expression correlated with distance from nearest PODXL+ arterial endothelial cells (Extended Data Fig. 5c–e, Supplementary Table 8). Further, NOTCH activation score correlated with perivascular fibroblast position (Spearman r=−0.71, p=0, Fig. 3e). Next, we explored how endothelial-derived NOTCH signaling can be propagated in fibroblasts to generate a positional gradient. In vitro, endothelial cells induced fibroblast expression of NOTCH3 and JAG1 in a NOTCH-dependent manner as addition of DAPT or siRNA-mediated NOTCH3 silencing in fibroblasts inhibited this induction (Extended Data Fig. 5f, g). Simultaneous upregulation of NOTCH3 and JAG1 in fibroblasts suggest that endothelial cells can induce a fibroblast positional gradient through a NOTCH3-JAG1 signal relay [16,17]. Consistent with this, we observed significant correlation between JAG1 and NOTCH3 expression in fibroblasts from synovial tissue (spearman r=0.17, p=7.25×10−285) and organoids (spearman r=0.26, p=7.21×10−105, Extended Data Fig. 5h).

Extended Data Figure 5.

Endothelial cells induce fibroblast NOTCH3 and JAG1 expression.

a, Fine-cluster analysis of synovial tissue endothelial cell scRNAseq: Arterial PODXL+ (orange) and venous DARC+ (blue) subtypes highlighted in UMAP projection. All non-vascular endothelial cells colored grey. n = 35,153 cells from 6 RA and 6 OA synovial tissues. b, Mean expression (colour) and percentage of cells with non-zero expression (size) of top gene markers distinguishing arterial and venous endothelial cells. c, Confocal microscopy images of synovial tissues where NOTCH3 (red), CD90(THY1) (green), and PODXL(blue) were visualized by immunofluorescence staining. d, Images from above (c) were then processed to segment cells and output their spatial location and mean intensity of all 3 markers. d, Cells labeled as arterial vascular endothelial cells were coloured in red. Non-arterial endothelial cells were coloured according to NOTCH3 expression, from low (grey) to high (blue). e, Spatial correlation between NOTCH3 expression and distance from nearest arterial vascular endothelial cell. Non-arterial cells in each image were used to analyze the function between arterial distance, measured by distance to nearest arterial vascular endothelial cells, and NOTCH3 expression. Spearman correlation values and significance were computed with the base R cor.test function. Cells were binned by frequency into groups of 50 cells along the x-axis. f-g, Flow cytometric analysis of fibroblasts and endothelial cell co-culture experiments. Cells are gated on CD31- CD90+ fibroblasts to exclude endothelial cells (n = 4 replicates). f, Representative flow plot of NOTCH3, CD90 (THY1), and JAG1 expression in fibroblast-only (left) and fibroblasts plus endothelial cells (ECs) co-culture (right). Colour indicates fibroblast expression level of JAG1. g, Mean fluorescence intensity (MFI) of fibroblast JAG1 expression in the presence and absence of endothelial cells (EC), gamma-secretase inhibitor DAPT, or pre-treatment with siRNA against control (si-CTRL) or NOTCH3 (si-NOTCH3). Boxplots summarize the median, interquartile range, and 95% quantile range. h, Biaxial plots of normalized JAG1 and NOTCH3 expression in CP10K from scRNAseq of synovial tissues (top, n = 12 donors) and synovial organoids (bottom, n = 3 organoids).

Fibroblast NOTCH activation increased in rheumatoid arthritis

In RA synovial tissues, we found a significant increase (2.9 fold, t-test p=0.005) in the proportion of NOTCH3+ fibroblasts (median 8.5%) compared to OA (median 2.9%, Fig. 3f). When we classified fibroblasts in our scRNAseq data as NOTCH-activated based on NOTCH activation score, we observed a significant increase in NOTCH-activated fibroblasts in RA versus OA (median 30.5% versus 4.6%, 6.6 fold, Wilcoxon p=0.015, Fig. 3g). We confirmed the enrichment of the NOTCH activation signature in RA synovial fibroblasts using an independent dataset generated from the Accelerating Medicines Partnership - Rheumatoid Arthritis (AMP-RA) cohort[3] (t-test p=1.2×10−4, Fig. 3h). Together, these data suggest a model in which endothelium-derived NOTCH ligands drive sublining fibroblast expansion through inductive NOTCH3 signaling in RA.

We next sought to determine whether NOTCH3 signaling could be specifically targeted in a mouse model of inflammatory arthritis[18]. We first examined Notch3 expression in mouse synovia by scRNAseq and determined that Notch3 expression is restricted to mural cells and fibroblasts (Fig. 4a, b). Co-staining of NOTCH3 protein and the leukocyte marker CD45 in arthritic mouse synovium confirmed perivascular NOTCH3 expression specifically in non-hematopoietic (CD45-) synovial cells (Extended Data Fig. 6a). Consistent with increased Notch receptor activation in perivascular cells, NOTCH activation score is significantly increased in these cells in arthritis compared to healthy state (p=0.0006, Fig. 4c). Adult Notch3−/− mice[19] exhibit normal joint architecture under steady-state conditions (Extended Data Fig. 6b). Upon arthritic serum transfer (K/BXN mice)[18,20], Notch3−/− mice exhibited significantly reduced arthritis disease index (p=3.01×10−16) and paw swelling (p=5.02×10−14) compared to wild-type mice (Fig. 4d, e). Further, twice weekly treatment with a NOTCH3 antagonizing antibody (anti-NRR3)[21] attenuated arthritis severity (p=3.47×10−8) and joint swelling (p=5.54×10−9) compared to treatment with an isotype control antibody (Fig. 4f, g).

Figure 4.

NOTCH3 blockade attenuates inflammatory arthritis.

a-c, scRNAseq analysis of mouse synovium. a, UMAP projection of 18,491 synovial cells where cells are colored by cell type (top) or Notch3 expression (bottom). b, Percentage of cells with non-zero Notch3 expression. c, NOTCH activation score in mural cells in mice of indicated group (n = 10 mice). Dotted grey line represents median NOTCH score from healthy wild-type mice. Clinical index (d) and paw swelling (e) in wild-type (n = 14) and Notch3−/− mice (n = 16) after serum transfer. Clinical index (f) and paw swelling (g) in isotype control antibody- (n = 20) or anti-NRR3- (n = 20) treated mice after serum transfer. Quantification (h) of cumulative joint histology score and representative hematoxylin and eosin staining (i) from wild-type (n = 14), Notch3−/− (n = 16), isotype control- (n = 10), or anti-NRR3- (n = 10) treated mice. S = synovium. B = bone. C = cartilage. Representative micro-CT image (j) and quantification of erosion scores (k) of joints from isotype control antibody- (n =10) or anti-NRR3- (n = 10) treated mice. Red arrows mark bone erosion. All boxplots summarize the median, interquartile range, and 95% quantile range. d-g, Mean and standard deviation (line) shown. Significance was determined by linear mixed effects models (d-g) and two-tailed t-test (h, k).

Extended Data Figure 6.

NOTCH3 expression in mouse synovia and the effects of NOTCH3 or NOTCH1 inhibition.

a, Serial sections of arthritic mouse synovia (n = 5 mice) showing representative hematoxylin and eosin staining (top) and immunofluorescence staining (bottom) of CD45 (red), NOTCH3 (green) and DAPI (blue). b, Representative hematoxylin and eosin staining of mouse joints from wild-type (top, n = 6) and Notch3−/− mice (bottom, n = 6). S = synovium. c, scRNAseq of wild-type, Notch3−/−, isotype control-, or anti-NRR3- treated mouse synovial cells (total of 18,491 cells from n = 5 mice per group). Cells were clustered and labeled in a joint analysis. In the 7 identified stromal populations, differential gene expression was performed for anti-NRR3 vs isotype control and for Notch3−/− vs wildtype. Each gene was plotted as a dot, representing the log fold-change. ECs = endothelial cells. Clinical index (d) and paw swelling (e) in IgG control antibody (dark blue, n = 20) or anti-NRR1 (light blue, n = 20) treated mice after K/BxN serum transfer. Significance of each treatment was determined by mixed effects linear models, controlling for time as a categorical fixed effect and mouse as a random effect. UMAP projections of normalized expression of (f) NOTCH1 in human primary tissue scRNAseq data (35,153 cells) and (g) Notch1 in mouse scRNAseq data (n = 18,491 cells)

Consistent with attenuation of NOTCH signaling, NOTCH activation scores in perivascular cells from Notch3−/− mice and anti-NRR3 treated mice were significantly reduced compared to wild-type (p=2.0×10−5) or isotype antibody-treated (p=0.0005) mice, respectively (Fig. 4c). Further, we found consistent, global transcriptomic changes in stromal cells between wild-type and either Notch3−/− mice and anti-NRR3 treated mice (Extended Data Fig. 6c). In both Notch3−/− mice and anti-NRR3 treated mice, we observed a reduction in pathological scores[20] when compared to wild-type mice (t-test p=0.0006) or treatment with isotype antibody (t-test p=0.0006, Fig. 4h, i). Further, anti-NRR3 treatment significantly reduced radiographic evidence of bone erosion (t-test p=0.0113, Fig. 4j, k).

Compared to anti-NOTCH3 blockade, antibody-mediated inhibition of NOTCH1 signaling with NRR1[22], the other Notch receptor identified in our receptor-ligand analysis, led to modest attenuation in arthritis severity (p=0.03) and paw swelling (p=0.03, Extended Data Fig. 6d, e), consistent with a role of NOTCH1 in synovial tissue inflammation[23-25]. The effect of anti-NRR1 on arthritis is likely not restricted to fibroblasts, as NOTCH1/Notch1 expression was detected in multiple lineages in both human and mouse synovia (Extended Data Fig. 6f, g).

Discussion

The function of individual cells and of tissues is dependent on their architectural design. Cellular zonation can reflect the division of labor among cells as a central design concept in tissue function. The correspondence between scRNAseq transcriptional gradients and cell positional identity has also been observed in other organ tissues[26,27]. Here, using scRNAseq we identified a transcriptional gradient among synovial fibroblasts that relates to their positional identity.

Fibroblasts play key roles in regulating immune responses in lymph nodes and in cancer[28]. In RA, synovial fibroblasts have long been considered an attractive therapeutic target, yet no therapy directly targeting synovial fibroblasts has been approved. Here, we found that NOTCH3 signaling contributes to mural cells and CD90(THY1)+ sublining fibroblast differentiation and that this step is required for the development of inflammatory arthritis. NOTCH3 signaling likely represents one of several important regulators of fibroblast identity. Indeed, NOTCH-activated fibroblasts in organoids partially recapitulated the intermediate states seen in native human and mouse synovial tissue fibroblasts (Extended Data Fig. 7a, b). Interestingly, the intermediate signature is enriched in a recently identified fibroblast subset[3] (Extended Data Fig. 7c). Future studies are required to determine if this subset represents a stable intermediate state or a transition state between lining and perivascular fibroblasts.

Extended Data Figure 7.

Characterization of intermediate fibroblast positional gene signature.

a, Gene enrichment scores for intermediate position-specific genes (x) vs NOTCH activation score (y) in synovial organoids (22,164 cells; n = 3 organoids per condition). ECs = endothelial cells. DAPT = organoids cultured in the presence of 10uM γ-secretase inhibitor DAPT. Fibroblast plus endothelial cell organoids (yellow) show evidence of increased NOTCH signaling but little enrichment in intermediate zone genes, as compared to fibroblasts that were cultured alone (blue) or co-cultured in the presence of DAPT (red). b, Enrichment of intermediate gene score in mouse synovial fibroblasts. Enrichment is lowest (red) in the previously identified lining and sublining zones and highest (blue) in the fibroblasts positioned in between the two zones, in the intermediate zone. c, Enrichment of pre-defined synovial fibroblast populations in the AMP-RA/SLE scRNAseq dataset (n = 1,844 cells).

Endothelium-derived NOTCH signaling in mural cell differentiation is well-appreciated in developing organ tissues[16,29] . Our studies underscore the importance and versatility of NOTCH signaling during pathological tissue remodeling. Our findings build on a body of knowledge that has implicated NOTCH signaling in RA pathogenesis[23,24,25,30] . While previous studies focused on the interaction of NOTCH1 signaling and cytokine activation in fibroblasts[23,24,30], our study advances previous understanding of NOTCH signaling in RA by identifying NOTCH3 as a critical receptor in synovial fibroblast differentiation and pathologic expansion in RA. These results provide a molecular basis by which stromal cells can be therapeutically targeted in RA by modulation of NOTCH3 signaling.

METHODS

Human research

Experiments involving human subjects were performed according to the Institutional Review Boards at Partners HealthCare. Patients with RA fulfilled the ACR 2010 Rheumatoid Arthritis classification criteria. For synovial tissue acquisition, patients from Brigham and Women’s Hospital or Hospital for Special Surgery undergoing joint replacement or synovectomy procedures were recruited under IRB-approved study. Synovial tissue were processed and cryopreserved as previously described [31] where each synovial tissue specimen was dissected into 1 to 5 mm fragments, then submerged in Cryostor CS10 (Sigma-Aldrich), then stored in −80 degrees prior to transferring to long-term storage in liquid-nitrogen. For RNA-seq and flow cytometric analyses, synovial tissues were thawed and disaggregated into single cell suspensions as previously described [31].

Flow cytometry and cell sorting

For human synovial tissue RNA-seq experiments, disaggregated human synovial cells were stained with antibodies against CD45 (HI30), CD235a (KC16), CD31 (WM59), THY1 (5E10), podoplanin (NZ1.3), CD146 (P1H12), CD34 (4H1) and in 1% BSA in Hepes-Buffered Saline (HBS, 20 mM HEPES, 137 mM NaCl, 3mM KCl, 1mM CaCl2) for 30 minutes. Propidium iodide (PI) or LIVE/DEAD (Invitrogen) viability dye was added to cell suspensions and cells were passed through a 100μm filter. For human scRNA-seq experiments, 20,000 synovial stromal cells (CD45-, CD235a-) were collected by fluorescence-activated cell sorting (BD FACSAria Fusion) into 0.4% BSA in PBS. For bulk RNA-seq, 1000 synovial fibroblasts were sorted based on surface markers as described directly into buffer TCL (Qiagen). For flow cytometric analysis of Notch pathway in cultured synovial fibroblasts, cells were stained with antibodies against NOTCH3 (MHN3–21) and JAG1 (MHJ1–152) in 1% BSA in Hepes-Buffered Saline (HBS, 20 mM HEPES, 137 mM NaCl, 3mM KCl, 1mM CaCl2) for 45 minutes at room temperature, followed by data acquisition (BD Fortessa) and analysis (Flowjo).

Single Cell RNA-sequencing

For scRNAseq experiments, viable cells in single cell suspension were resuspended in 0.4% BSA in PBS at a concentration of 1,000 cells per ul. 7,000 cells were loaded onto a single lane (Chromium chip, 10X Genomics) followed by encapsulation in lipid droplet (Single Cell 3’ kit, 10X Genomics) followed by cDNA and library generation per manufacturer protocol. For multiplexed scRNAseq experiments[32], cells with stained with cell hashing antibodies (TotalSeq, Biolegend) prior to cell capture. cDNA libraries were sequenced to an average of 50,000 reads per cell using Illumina Nextseq 500. scRNA-seq reads were processed with Cell Ranger v2.1, which demultiplexed cells from different samples and quantified transcript counts per putative cell. Quantification was performed using the STAR aligner against the hg19 transcriptome. For bulk RNA-sequencing experiments, full-length cDNA and sequencing libraries were performed using Illumina Smart-seq2 protocol [33]. Libraries were sequenced on MiSeq from Illumina to generate 35 base paired-end reads. Reads were mapped to the hg19 transcriptome using kallisto 0.42.4 and transcriptional levels of genes were quantified with the Log2(TPM+1) (Transcripts Per Kilobase Million) metric.

RNA-seq quality control and pre-processing

After filtering out low-quality cells (<1000 unique genes, >10% mitochondrial reads) and potential doublets (>12,000 UMIs), 35,153 cells from primary tissue and 6,412 cells from organoid experiments were further analyzed. With these high quality cells, we used a standard normalization procedure. We normalized each cell to 10,000 reads, by multiplicative scaling, then log scaled the normalized data. We then identified the top 1000 variable genes, ranked by coefficient of variation, within each individual 10X experiment. We pooled these genes to form the variable gene set of the analysis. Using only the variable genes, we mean centered and variance 1 scaled the genes across the cells. Note that this was done in the aggregate matrix, with all cells, rather than within each dataset separately. With these values, we performed truncated SVD keeping the top 100 eigenvalue-scaled eigenvectors. Visualization of single cells was performed using a UMAP[34] projection into 2 dimensions. We corrected for donor-specific effects in the primary-tissue only analysis using the Harmony algorithm [6] with default parameters. For bulk RNAseq, we removed 4 samples with low read depth (<400,000 estimated total gene counts).

Population identification in scRNAseq

We labeled major cell populations in the primary synovial tissue data using clustering and marker analyses. We clustered with the Leiden algorithm (resolution=0.005, all other parameters default) on the shared nearest neighbor graph. The graph was constructed with k=30 nearest neighbors with an exact nearest neighbor search using the RANN R package. Edges in the shared nearest neighbor graph with weight<1/15 were pruned. We identified markers for each cluster with the Mann-Whitney U-test, filtering for over-expressed genes in each cluster. We identified vascular endothelial cells as CD34+/VWF+, mural cells as MCAM+/ACTA2+, lining fibroblasts as PDFGRA+/PRG4+, and sublining fibroblasts as PDFGRA+/THY1+.

To identify subpopulations of endothelial cells, we subclustered the previously identified vascular endothelial cluster (Fig. 1a) into arterial and venous populations. To do this, we performed graph based clustering on the 7019 endothelial cells to get 5 fine grained clusters. For each cluster, we performed a gene set enrichment analysis on arterial-specific and venous-specific gene lists. For the enrichment analysis, we used the R package fgsea, run with 10,000 permutations. The gene lists were constructed in the following way. We downloaded differential expression results from Vanlandewijck et al26, in which the authors compare arterial and venous cells from mouse brain vasculature. We selected differentially expressed genes with p<0.001 and absolute log fold change>1. We then mapped mouse genes onto human orthologs. This filter yielded final genesets of 253 arterial genes and 181 venous genes. Finally, we performed differential expression on the 5 clusters, with the Wilcoxon rank sum test. We enriched these differential expression results with the 2 genesets to identify 1 venous cluster and 2 arterial clusters. The other 2 clusters, with a total of only 93 cells, were indeterminate. We did not assign a label to these clusters. The final assignments represent the 1 venous cluster and 2 merged arterial clusters.

Projection onto single cell data

We projected both single cells from the scRNAseq organoid experiments as well as bulk RNAseq samples from the ex vivo stability experiments onto the primary synovial tissue single cell datasets. For single cells, we first merged the organoid and synovial datasets using Harmony, correcting for both dataset and culture vs primary identity (theta=0.5 for both variables). Within this Harmonized space, we identified the 30 nearest synovial cells for each organoid cell, and used a radial basis function kernel (sigma=0.1) with cosine distance to build neighborhood probability distributions. With these distributions, we assigned, to each organoid cell, mean values of UMAP coordinates, DDRTree coordinates, and trajectory positional scores. We used the same approach to map mural cells from single cell synovial tissue datasets onto the DDRTree coordinates. For the bulk RNAseq samples, we adopted a similar strategy. Instead of using Harmonized coordinates to compute distance, we used the cosine distance in gene expression space. In particular, we used the intersection of scRNAseq variable genes and all available bulk genes. The expression profiles were log-counts-per-10,000 in the scRNAseq space and log TPM in the bulk RNAseq space.

Trajectory analysis

We defined the 1D transcriptional gradient of fibroblasts and mural cells using trajectory analysis on the primary human scRNAseq datasets. In particular, we inferred the trajectory using fibroblasts only, as mural cells appeared to form a distinct cluster, whereas fibroblasts appeared to have a more smooth structure. We then projected mural cells onto the trajectory using the single cell projection method described above. In previous analyses, we found that gene expression was influenced to a large degree by donor specific variation. To avoid biasing trajectory analysis with this variation, we used the donor corrected low dimensional embedding, learned by Harmony, as input into the trajectory algorithm. Specifically, we used DDRTree, the reversed graph embedding algorithm, implemented in the DDRTree R package. DDRTree is the core algorithm for trajectory analysis in the Monocle2 R package[11], a popular single cell analysis trajectory analysis pipeline. Instead of using gene expression values, as in Monocle2, we input the donor-corrected Harmony cell embeddings into DDRTree, using default DDRTree parameters. We attempted to infer pseudotime using the Monocle2 pseudotime algorithm. Unfortunately, this algorithm did not scale to the large number of cells in our analysis. We noted that the fibroblast trajectory did not involve branching events. Thus, principal curve analysis was sufficient to infer ordering along the inferred trajectory. We used the princurve R package to infer position (aka pseudotime) along the DDRTree trajectory.

scRNAseq background estimation

We developed a computational strategy to evaluate the contribution of non-cellular, background mRNA to each cell. At a high level, this strategy works in three steps. First, we define canonical cell types in the data, based on previous analysis, and compute a mean expression profile for each type. Second, we define a synthetic population that represents the background mRNA levels within a sample. We refer to the cluster-based and background expression profiles as archetypes. Third, we probabilistically assign each cell to one of the archetypes. We specify the details for each step below.

It is important to note that the background mRNA profile is different for each sample. Thus, all steps were performed separately within each sample. For step 1, we defined archetypes for lining fibroblasts, sublining fibroblasts, mural cells, blood vessel cells, lymph vessel cells, and (outlier) leukocytes. For fibroblasts, we excluded cells with intermediate positional identities, defined here as between 15 and 85. The transcriptional profile for each biological archetype was the mean normalized gene expression of cells in the corresponding cluster. For step 2, we followed the approach of Young and Behtaji (2018) and reasoned that background mRNA expression can be estimated by counting UMIs in empty droplets. Since these droplets contain no cells, the only way for them to contain mRNA molecules is through contamination of ambient molecules. We identified empty droplets as those with at most 10 UMIs. We summed gene counts over all empty droplets, normalized to 10,000 counts and performed a log transformation. The resulting gene expression profile is comparable to the normalized profiles of the cluster-based archetypes. For step 3, we used a cosine distance based radial basis function (sigma=0.1) to compute the similarity of each cell to each of the 7 archetypes. For each cell, we normalized the similarities to sum to 1, thus creating a unique assignment probability vector from cells to archetypes.

To evaluate the archetype assignment method, we used classification results in the analysis above to construct a confusion matrix. Here, we evaluated the probability of a cell of known type to be classified as one of the 7 archetypes. To have both positive and negative examples, we included samples that represented true cells (based on cluster-based labels) as well as samples that more likely represented empty droplets with background gene expression. To avoid overlap with step 2 of the analysis above, these test empty droplets were selected to have between 100 and 300 UMIs. While these samples may represent undersampled cells, they are more likely to be enriched in background mRNA. In the confusion matrix, we call these samples Background. The assignment probabilities of all known types was high: 93% for Background, 98% for blood vessel cells, 99% for leukocytes, 94% for lining fibroblasts, 96% for lymphatic vessel cells, 95% for mural cells, and 89% for sublining fibroblasts.

Intermediate trajectory signature

We analyzed non-linear expression patterns along the fibroblast trajectory to identify a core set of 71 genes associated with intermediate position along the trajectory. To do this, we divided the trajectory into two parts ([0,50] and [50, 100]) and filtered for genes correlated (Pearson |r|>0.1) with either or both sub-trajectories. For each of the resulting 949 genes, we binned (100 quantiles) expression along pseudotime and fit the mean normalized expression in each bin to trajectory position with natural cubic splines, with the R sm.ns function with 3 degrees of freedom. We then clustered the fit spline values with hierarchical clustering, using the R hclust method with default parameters, and plotted the clustered genes against trajectory position in a heatmap. The genes, cluster colors, and cluster labels are reported in Supplementary Table 2. We performed pathway enrichment of the 71 genes associated with intermediate position (green and yellow clusters) using the enrichR package and the GO Biological Processes 2018 genesets. Pathways with FDR<20% are reported. We used these 71 genes to identify gene enrichment signatures in 3 independent datasets: synovial organoid scRNAseq, mouse scRNAseq, and synovial fibroblast scRNAseq data from the recent Accelerating Medicines Partnership (AMP) publication[3]. In each case, the signatures were computed for individual cells as the sum of z-scores, computed by mean centering and standard deviation scaling of normalized (log CPM in AMP, log CP10K in organoid and mouse cells) expression data. Cluster assignments for the AMP data analysis were obtained directly from the authors.

Immunofluorescence staining and imaging analysis

Imaging of the OA and RA tissue was carried out in the University of Birmingham Medical School Imaging Suite. OA and RA synovial tissues were frozen and embedded in OCT (Fisher). 6um sections were cut using a Leica CM1950 onto X-tra Slides (Leica Biosystems) and fixed in acetone at 4°C for 20 minutes. GGT5 stains used FFPE sections antigen retrieved in pH9 Tris-EDTA 0.05% Tween 20 buffer for 1 hour at 96°C. Slides were rehydrated in PBS before blocking with 10% normal horse serum. For lubricin, Notch3 and CD55 fluorescent staining and THY1 staining was completed before the remaining antibodies were added. In brief, this entailed THY1 (sheep polyclonal, AF2067, Bio-Techne), followed by donkey anti goat FITC (Jackson Immuno Research) and then blocking with 10% normal sheep serum and purified sheep IgGs for 30 minutes. Other primary antibodies were PRG4 (IgG1, clone 9G3, Merck Millipore), MCAM (IgG2a, clone SHM-57, Biolegend), VWF (rabbit polyclonal, Dako), Notch3 (MHN3–21, BioLegend), PODXL (rabbit, clone EPR9518, Abcam), CD55 (IgG1, clone 143–30, Bio-Techne), GGT5 (HPA008121, Sigma), PDPN (NZ-1.3, eBioscience), THY1 (IgG2a, clone Thy-1A1 Bio-Techne), and CD31 (IgG1, JC70A, Dako) were then applied in 10% normal sheep serum. Primary antibodies for mouse synovial tissues were NOTCH3 (clone AF1408, R&D systems) and CD45 (clone D3F8Q, Cell Signaling). Appropriate isotype controls were used on a separate section. Secondary antibodies included goat anti rabbit Alexa Fluor 594, goat anti rabbit Alex Fluor 546, goat anti mouse IgG1 biotin, FITC or Alexa Fluor 546, goat anti mouse IgG2a or IgG1 Alexa Fluor 647 (all Thermo Fisher except from IgG1 FITC from Southern Biotech) were added following washing. Finally, goat anti FITC Alexa Fluor 488, Streptavidin-Alexa Fluor 546 and Hoescht 33258 were applied if required (all Thermo Fisher). Slides were mounted in prolong diamond (Thermo Fisher) and stored at −20°C before imaging. Images were obtained using a Zeiss LSM 880 in linear deconvolution or sequential mode and were analysed using Definiens developer XD. In brief, areas of tissue were outlined prior to user lead automatic identification of the lining layer if possible (as defined by PRG4 or CD55 staining), blood vessel structures, holes and sublining layer. Areas defined as holes were disregarded in further analysis. Nuclei were identified and cells grown out from those structures. The software then provided data on cellular location, intensity of staining and relations to the identified structures.

RNA scope was carried out on formalin fixed paraffin embedded tissues sections of OA and RA synovium. Assays for Notch 3 (558991) was used according to the manufacturer’s RNA scope 2.5 HD assay red kit instructions (322360, all Bio-Techne). Following the last wash for the RNAscope protocol, slides were transferred to distilled water for 5 minutes, then blocked with Bloxall (Vector Laboratories) for 10 minutes followed by 10% normal horse serum in Tris buffer for 10 minutes. Sheep anti THY1 (as above) was added overnight at 4°C followed by donkey anti-sheep HRP (713–035-147, Jackson Immuno Research) and ImmPACT® DAB Peroxidase (HRP) Substrate (Vector laboratories). Nuclei were visualised using Vector® Hematoxylin QS and slides were mounted using VectaMount® Permanent Mounting Medium (both Vector Laboratories). Images were obtained using the Zeiss Axio Scan and analyzed in Zen Blue (both Zeiss).

For downstream analysis, we identified individual vascular endothelial cells and computed the normalized CD90:PRG4 intensity ratio for all non-vascular endothelial cells. All analyses described in this section were performed within each image separately, thereby avoiding image-to-image variation in intensity levels. Raw intensities for each marker were log transformed and z-scored within each image. Finally, we thresholded the z-scores at −3 and 3 to dampen the impact of extreme outliers. Following this transformation, we computed the CD90:PRG4 intensity as the (thresholded) difference between the scaled CD90 and scaled PRG4 intensities. We labeled vascular endothelial cells by combining two pieces of information: manually curated blood vessel regions and the intensity of the canonical endothelial cell marker VWF. The reasoning was that per-cell VWF intensity was subject to imaging noise. On the other hand, not all cells in blood vessel regions are necessarily VWF+. Combining these two approaches decreased the chance of erroneously identifying non-endothelial cells. We labeled cells as VWF+ with 3-component Gaussian mixture models, implemented in the R package mclust. This model estimated densities of positive, negative, and intermediate VWF intensities. Cells assigned to the high VWF intensity cluster with at least 95% confidence were labeled VWF+. Finally, we computed the distance between non-endothelial cells to their nearest endothelial cell using cell-center to cell-center Euclidean distance.

For NOTCH3 intracellular domain (NICD3) staining, RA synovial tissue were stained with an antibody against NOTCH3 intracellular domain (V1662) as previously described[35]. For NOTCH3 total protein staining, all staining took place at RT. Slides were rehydrated in Tris buffered saline before blocking with Bloxall for 10 minutes (Vector labs) and 10% normal horse serum. Anti-NOTCH3 (clone MHN3–21 IgG1, BioLegend) was added for one hour before washing and addition of goat anti IgG1 biotin for 45 minutes. SA-HRP was added for 30 minutes prior to developing with ImmPACT DAB reagent (both Vector Laboratories). Sections were blocked with Bloxall again before staining with rat anti PDPN (NZ-1.3, Thermo Fisher). Goat anti rat alkaline phosphatase (Thermo Fisher) was applied prior to ImmPACT Vector Red and Haematoxylin QS (both Vector Laboratories). Slides were mounted in VectaMount (Vector Laboratories) prior to imaging with a Zeiss AxioScanZ.1 slide scanner.

Synovial tissue organoid studies

Early passage synovial fibroblasts cell lines (passage 3 to 5) were used for organoid experiments. Human Umbilical Vein Endothelial Cells (HUVECs) were purchased from Lonza and expanded in the presence of growth factors (EGM-2, Lonza). For confocal imaging analysis of organoids, fibroblasts and HUVECs were first labeled with fluorescent dyes PKH67 and PKH26 (both from Sigma-Aldrich), respectively. Fibroblast organoids were created using 200,000 cultured synovial fibroblasts suspended in a single droplet of Matrigel (Corning) as previously described [36,37]. For fibroblast plus endothelial cell organoids, 100,000 synovial fibroblasts and 100,000 HUVECs were combined in a single 35 ul droplet of Matrigel, seeded on polyHEMA (Sigma-Aldrich) coated plates, and cultured in EGM-2 media (Lonza) for 14–21 days. For scRNAseq analysis of synovial tissue organoids, organoid cells were isolated on day 14 or 21 through enzymatic dissociation using 1mg/mL of dispase II (Roche) in DMEM.

We performed scRNAseq on organoid experiments in 2 stages. In the first experiment, we pooled cells from fibroblast-only organoids into one scRNAseq batch and cells from the fibroblast with endothelial cell condition into a second batch. In an independent experiment, we used Cell Hashing[32] to multiplex 3 technical replicates from 3 conditions into a single scRNAseq experiment. The three biological conditions were fibroblasts cultured alone, fibroblasts cultured with endothelial cells, and fibroblasts cultured with endothelial cells and DAPT. Through this multiplexing strategy, we were able to recover the technical replicate and condition ID of each cell. This helped to minimize batch effects among the 9 samples.

To maximize statistical power, we combined cells from both stages into a single joint analysis. We harmonized over the 3 batches (scRNAseq runs) and accurately matched fibroblast-alone (blue) and fibroblasts with endothelial cells (orange) conditions from the two sets of experiments. Importantly, we did not need to harmonize over the 9 samples assayed in the second experiment, since they were all multiplexed into a single batch. In the joint analysis, we performed the standard steps described earlier for QC filtering, library size normalization, variable gene selection, PCA, Harmony, and UMAP. Through the combination of Cell Hashing and Harmony, we are able to attribute differences among culture conditions to real biological effects rather than confounding batch effects.

Cell culture studies

Synovial fibroblasts derived from RA synovial tissue were cultured in DMEM 10% fetal bovine serum as previously described [38]. For fibroblast Notch activation experiments, cell culture plates were coated with recombinant DLL4 or JAG1-FC (R&D systems) overnight at 4 degrees, at 5 μg/ml or at the concentration indicated for each experiment. Fibroblasts were then seeded on DLL4-, JAG1-FC- , or vehicle-coated plates for 72 to 120 hours. For cytokines and growth factors stimulation experiments, recombinant proteins were purchased (R&D systems) and reconstituted in DMSO or PBS per manufacturer’s instructions. Recombinant protein were diluted in media and then added to cultured fibroblasts at the following concentrations: TGFB1 10 ng/ml, BMP7 10 ng/ml, ACTIVIN A 10 ng/ml, WNT3A 100 ng/ml, WNT5A 100 ng/ml, EGF10 ng/ml, TNF 10 ng/ml, IFN-gamma or PDGF-BB 50 ng/ml. For flow cytometric analysis of Notch pathway in fibroblast-endothelial cell co-culture experiments, HUVECs (Lonza), passage 3 – 7, were cultured in a 1:1 ratio with synovial fibroblasts in EGM2 media for 72–96 hours in the presence or absence of 10 μM DAPT as indicated. For siRNA experiment, all siRNAs (Silencer Select) were purchased from Life Technologies. Fibroblasts were transfected with siRNA by reverse transfection at 30 nM using RNAi Max reagent (Life Technologies) as previously described [38].

NOTCH activation score

We used the NOTCH stimulation experiments described above to define a fibroblast-specific Notch activation signature. For each ligand, JAG1 and DLL4, we used limma to perform differential expression analysis against the PBS treated samples. Genes with a log fold change greater than 0 and a benjamini-hochberg corrected p value less than 0.05 were considered upregulated in response to ligand stimulation. We refer to these ligand-specific upregulated genes as the stimulation signatures. In both scRNAseq and bulk RNAseq data, NOTCH activation scores were computed as the sum of the expression of signature genes, either in log CP10K (scRNAseq) or log TPM (bulk RNAseq). We filtered on significantly (moderated t-test fdr<0.05) upregulated genes to define a combined gene signature of 118 genes, 45 genes from JAG1 and 111 from DLL4 stimulation (Supplementary Table 9).

Western blot

Protein lysates from fibroblasts were prepared using RIPA buffer, separated on Criterion Precast Gels (Bio-Rad) 4–15% Tris-HCl, followed by transfer to PVDF Membrane (Bio-Rad). Membranes were probed with the following primary antibodies in TBS-T and 5% milk for 12 to 16 hours: Anti-Notch1 (R&D AF33471), Anti-Notch3 (V1662) 1ug/ml, GAPDH (1:1000, Cell Signaling). Secondary antibodies (1:10,000) were from Jackson ImmunoResearch. Chemiluminescent detection reagents were from Bio-Rad.

Real-time quantitative PCR

cDNA was synthesized with QuantiTect Reverse Transcription kits (QIAGEN). qPCR was performed with Brilliant III Ultra-Fast SYBR Green qPCR master mix (Agilent Technologies) on a Mx3000 (Stratagene). The following primers were used in this study: GAPDH forward; GTCTCCTCTGACTTCAACAGCG, GAPDH reverse; GTCATGTTCAGGTCCAACTCGG, THY1 forward;GAAGGTCCTCTACTTATCCGCC, THY1 reverse; TGATGCCCTCACACTTGACCAG.

Clinical evaluation of synovitis by doppler ultrasound

The joint to be biopsied was assessed using ultrasound immediately prior to the procedure using a Siemens Acuson Antares scanner (Siemens PLC, Bracknell, UK) and multifrequency (5–13MHz) linear array transducers. Synovitis and power Doppler (PD) positivity were defined using consensus OMERACT definitions. Power Doppler ultrasound variables were graded on 0–3 semi-quantitative scales.

Animal studies

Experiments involving animal subjects was performed in accordance with BWH IRB-approved protocols. Notch3−/− mice (B6;129S1-Notch3tm1Grid/2J) were purchased from Jackson Laboratory and maintained as homozygous null lined. Adult, male, age-matched 129/B6 and C57Bl/6 mice were purchased from Jackson laboratory for serum-transfer arthritis studies. For antibody treatment studies, anti-NRR3 (20 mg/kg, Genentech), anti-NRR1 10 mg/kg, Genentech) or mouse IgG2a isotype control (20 mg/kg, Bio X-cell) were given intraperitoneally twice weekly for the duration of the study. For serum-transfer studies, arthritis was induced in 6–8 week old male mice by i.v. injection of 100 μL of pooled serum from K/BxN mice. Disease severity was assessed by daily clinical index and joint thickness measurements using a caliper. For single cell analysis of mouse synovia, synovial tissues were dissected from knee joints. Pooled synovial tissue from 10 mice from each treatment group were enzymatically dissociated as previously described[31]. Synovial cells were stained with LIVE/DEAD stain (ThermoFisher) and live synovial cells were isolated by flow cytometric sorting. 7,000 – 10,000 live synovial cells were then used for droplet-based single cell RNAseq (10X Genomics) as described above.

Histological analysis in mouse and organoid studies.

For histological analysis of animal joints, paws were fixed with 4% paraformaldehyde, decalcified, embedded in paraffin, and stained with hematoxylin/eosin as previously described [39]. Histopathological changes were scored according to an adapted method described previously [40]. The following five parameters were used and scored based on a scale from 0–3: synovitis, pannus formation, exudate (infiltration of inflammatory cells into the joint space), cartilage destruction, and bone erosion. Left and right hind paw scores for each of the five categories as well as left and right composite histopathological scores were averaged for each mouse. These scores were then averaged across each experimental group. For histological analysis of synovial organoids, organoids were fixed with 4% paraformaldehyde, embedded in paraffin, and stained with hematoxylin/eosin as previously described [41].

Micro–computed tomography

Fore paws were scanned using a Scanco mCT-35 instrument at an isotropic voxel size of 7mm, and 3-dimensional images were generated with software supplied by the manufacturer. The severity of periarticular erosions was determined blindly using the 3-dimensional images with a semiquantitative scale as described [42]. Four anatomical sites in the wrist joint were scored on a scale from 0–3: the distal ulnar epiphysis and the bases of the third, fourth, and fifth metacarpals. The left and right fore paw scores were averaged as the composite score for each mouse. Individual mouse scores were then averaged within experimental groups for final erosion scores.

Ligand-Receptor analysis

We used a heuristic analytical approach to prioritize pairs of ligands and receptors potentially involved in intercellular signaling between vascular endothelial and mesenchymal (mural cell and fibroblast) populations. Simply, we performed marker analysis, as described above in the section. We considered a set of previously curated ligand-receptor pairs from the Database of Interacting Proteins[43]. For each ligand-receptor pair, we required that the ligand be significantly (FDR<10−10) overexpressed in endothelial cells and the receptor overexpressed in mural cells. Additionally, ligands had to be expressed (>=1 UMI) in at least 10% of endothelial cells and receptors in at least 10% of mural cells.

Statistical analysis

Statistical comparisons were performed as indicated in the figure legends. P values <0.05 were considered significant after adjusting for multiple testing using the Bonferroni correction for ANOVA. In the analysis of the effect of Notch perturbation (i.e. treatment) in vivo on clinical outcomes (Fig. 4d–g), we used linear mixed effect modeling to account repeated measures across time. Specifically, we controlled for timepoint (as a categorical factor) as a fixed effect and animal ID as a random effect. To evaluate the significance of treatment, we used the likelihood ratio test, comparing the full linear mixed model (outcome ~ 1 + treatment + factor(timepoint) + (1|animal)) to the null model (outcome ~ 1 + factor(timepoint) + (1|animal)).