Fibroblast Subpopulations in Systemic Sclerosis: Functional Implications of Individual Subpopulations and Correlations with Clinical Features

Abstract

Fibroblasts constitute a heterogeneous population of cells. In this study, we integrated single-cell RNA-sequencing and bulk RNA-sequencing data as well as clinical information to study the role of individual fibroblast populations in systemic sclerosis (SSc). SSc skin demonstrated an increased abundance of COMP+, COL11A1+, MYOC+, CCL19+, SFRP4/SFRP2+, and PRSS23/SFRP2+ fibroblasts signatures and decreased proportions of CXCL12+ and PI16+ fibroblast signatures in the Prospective Registry of Early Systemic Sclerosis and Genetics versus Environment in Scleroderma Outcome Study cohorts. Numerical differences were confirmed by multicolor immunofluorescence for selected fibroblast populations. COMP+, COL11A1+, SFRP4/SFRP2+, PRSS23/SFRP2+, and PI16+ fibroblasts were similarly altered between normal wound healing and patients with SSc. The proportions of profibrotic COMP+, COL11A1+, SFRP4/SFRP2+, and PRSS23/SFRP2+ and proinflammatory CCL19+ fibroblast signatures were positively correlated with clinical and histopathological parameters of skin fibrosis, whereas signatures of CXCL12+ and PI16+ fibroblasts were inversely correlated. Incorporating the proportions of COMP+, COL11A1+, SFRP4/SFRP2+, and PRSS23/SFRP2+ fibroblast signatures into machine learning models improved the classification of patients with SSc into those with progressive versus stable skin fibrosis. In summary, the profound imbalance of fibroblast subpopulations in SSc may drive the progression of skin fibrosis. Specific targeting of disease-relevant fibroblast populations may offer opportunities for the treatment of SSc and other fibrotic diseases.

Keywords: Dermal fibrosis; Disease outcomes; Fibroblasts; SSc; Subpopulations.

Figure 1:. The abundance of individual fibroblast subpopulation signatures in the skin of systemic sclerosis (SSc) patients and matched healthy controls (HC)

a: Relative proportion of individual fibroblast subpopulation in patients with diffuse cutaneous SSc (dSSc) (n=48) and HC (n=33) in the PRESS cohort. **** P<0.0001 by unpaired t-test or Wilcoxon test. b: Fibroblast subpopulation signatures in patients with dSSc (n=43), limited cutaneous SSc (lSSc) (n=18) and HC (n=36) from GENISOS cohort. * P<0.05, ** P<0.01, *** P<0.001, **** P<0.0001 by unpaired t-test or Wilcoxon test. c: Geneset enrichment analysis (GSEA) plot of SFRP4/ SFRP2+ and PRSS23/SFRP2+ fibroblast gene sets in the skin of patients with dSSc or lSSc compared to HC, and dSSc compared to lSSc. **** P<0.0001. d: Proportion of fibroblast subpopulation signatures according to the intrinsic molecular subsets in the GENISOS cohort: normal-like SSc subset (n=21), inflammatory SSc subset (n=22), fibroproliferative SSc subset (n=18) as compared to normal-like HC. * P<0.05, ** P<0.01, *** P<0.001, **** P<0.0001 by unpaired t-test or Wilcoxon test. e: GSEA plot of the SFRP4/SFRP2+ and PRSS23/SFRP2+ fibroblast gene sets in the SSc intrinsic molecular subsets as compared to HC. *P<0.05, ****P<0.0001. f: Relative ratios of fibroblast subpopulations in the PRESS cohort, the GENISOS cohort and SSc scRNA–seq datasets.

Figure 2:. Transcriptome-based functional annotation of individual fibroblast subpopulations.

a: Functions of individual fibroblast subpopulations as predicted by GO enrichment for biological process (BP). GO terms with p-values <0.05 were considered statistically significant. Circles are color coded according to the p-values based on Fisher’s exact test. The size of the circles is proportional to the percentage of genes in the GO terms regulated in this subtype. b: Violin plots showing the mRNA levels of collagen genes (COL1A1, COL1A2, COL3A1) and myofibroblasts marker genes (ACTA2, TAGLN, CTGF, ASPN, ADAM12, FNDC1 and PDGFRA).

Figure 3:. Cell-cell interactions between CCL19+ fibroblasts and immune cell populations and ANGPTL7+ fibroblasts and endothelial cells

a: Circle plot of ligand-receptor interactions between CCL19+ fibroblasts and immune cells, as predicted using CellChat analysis. The lines inside the circle plot indicate the relative signal strength of the ligand and receptor, thicker line indicates a stronger signal, the number represents the weights of the interaction strength. Circle sizes are proportional to the number of cells in each cell group. b: Bubble plots showing the ligands-receptor interactions between CCL19+ fibroblasts and different immune cells. Circles are color coded according to the communication probability; the size of circles represents the respective p-values. c: Circle plot of ligand-receptor interactions between ANGPTL+ fibroblast and ACKR1+ endothelial cells. The lines inside the circle plot indicate the relative signal strength of the ligand and receptor, thicker line indicates a stronger signal, the number represents the weights of the interaction strength. Circle sizes are proportional to the number of cells in each cell group. d: Bubble plots showing the ligands-receptor interactions between ANGPTL7+ fibroblasts and endothelial cells. Circles are color coded according to the communication probability; the size of circles represents the respective p-values.

Figure 4:. Associations between the proportion of individual fibroblast subpopulation signatures and clinical features of SSc patients.

a: Correlation of the proportions of individual fibroblast subpopulation signatures with the local skin score at the site of biopsy, the total mRSS, and the skin thickness progression (STPR) at the time of biopsy in the PRESS cohort (left). Correlation of the proportions of individual fibroblast subpopulation signatures with clinical parameters in the GENISOS cohort (right). Dots indicate the coefficient; error bars represent the confidence interval. Colors indicate the individual fibroblast subpopulations. * P<0.05, ** P<0.01, *** P<0.001 by Pearson (for normally distributed variables) or Spearman’s Rank correlation (for non-normally distributed variables) coefficient. b: Correlation of the proportion of fibroblast subpopulation signatures with a previously established fibroblast signature score (Skaug et al., 2020) and histological parameters such as myofibroblast counts or relative collagen thickness(Skaug et al., 2020). Dotted lines separate SFRP4/SFRP2+ and PRSS23/SFRP2+ fibroblasts calculated by ssGSEA from the other populations determined by CIBERSORTx. * P<0.05, ** P<0.01, *** P<0.001 by Pearson or Spearman’s Rank correlation coefficient.

Figure 5:. Machine learning models to classify SSc patients based on fibroblast subpopulations abundance and clinical indices.

a: Boxplot showing the accuracy score in the different machine learning models to predict progression of skin fibrosis in SSc patients. b: The Mean Decrease Gini represents the relative importance of variables to differentiate SSc patients with skin progression from those with stable diseases computed by the feature selection technique in the random forest model. c: MDS plots displayed that random forest model assigned the patients into two clear separation clusters, one of which was rather exclusively composed of patients without progression of skin fibrosis.

Figure 6:. Summary of the changes of individual fibroblast subpopulations in SSc skin and their main functions.

Illustration assistance provided by Katie Vicari, Katie Vicari Scientific Illustration, Bethpage, New York.