A single-cell and spatial genomics atlas of human skin fibroblasts reveals shared disease-related fibroblast subtypes across tissues

Abstract

Fibroblasts sculpt the architecture and cellular microenvironments of various tissues. Here we constructed a spatially resolved atlas of human skin fibroblasts from healthy skin and 23 skin diseases, with comparison to 14 cross-tissue diseases. We define six major skin fibroblast subtypes in health and three that are disease-specific. We characterize two fibroblast subtypes further as they are conserved across tissues and are immune-related. The first, F3: fibroblastic reticular cell-like fibroblast (CCL19⁺CD74⁺HLA-DRA⁺), is a fibroblastic reticular cell-like subtype that is predicted to maintain the superficial perivascular immune niche. The second, F6: inflammatory myofibroblasts (IL11⁺MMP1⁺CXCL8⁺IL7R⁺), characterizes early human skin wounds, inflammatory diseases with scarring risk and cancer. F6: inflammatory myofibroblasts were predicted to recruit neutrophils, monocytes and B cells across multiple human tissues. Our study provides a harmonized nomenclature for skin fibroblasts in health and disease, contextualized with cross-tissue findings and clinical skin disease profiles.

Fig. 1. Identification of fibroblast subtypes in healthy skin.

a, Overview of study methodology, including skin atlas integration to delineate fibroblasts, construction of a healthy/nonlesional reference, mapping of 23 diseases to the healthy reference atlas and downstream analysis for cross-tissue comparison. b, Uniform Manifold Approximation and Projection (UMAP) of healthy and nonlesional skin fibroblasts colored by fibroblast subtype. DS, dermal sheath; DP, dermal papilla. c, Dotplot of marker gene expression for healthy fibroblasts. ‘All’ indicates a marker for a general population, but which contains subtypes. Supplementary Data Fig. 1a provides additional differentially expressed genes for fibroblast subtypes. d, Summary of skin fibroblast subtypes in healthy steady-state tissue. Illustrations in a and d were partly created using BioRender.com.

Fig. 2. Skin fibroblasts occupy unique spatial and functional niches.

a, Spatial location of fibroblast subtypes in microenvironments (cell2location abundance predictions (10x Genomics Visium)) in a single section of healthy human skin (left). Histopathological annotation of tissue microenvironments (right). b, Spatial location of fibroblasts at single-cell resolution (10x Genomics Xenium 5000-gene panel) for skin sections from nonlesional skin of atopic dermatitis (noninflamed (left) and noninflamed post-treatment (right)), colored by cell type. c, Summary of fibroblast niches: Xenium cell types overlying H&E-stained image (manual approximations).

Fig. 3. Prototype meta-learning to identify disease-adapted and disease-specific populations.

a, Overview of reference-mapping approach used for integration, where lesional/diseased data were mapped using a pretrained model. b, UMAP of scPoli embeddings colored by predicted cell-type labels and relabeled populations after re-clustering. c, Dotplot of marker gene expression for disease-adapted and disease-specific fibroblast populations. Supplementary Data Fig. 1c provides additional differentially expressed genes for disease-associated fibroblast subtypes. d, Density of cells in embedding by site status. e, Gene expression in F1 and F3 fibroblasts from health and disease, including differentially expressed genes in lesional/diseased states. f, Summary of disease-adapted and disease-specific populations. g, Feature maps for genes associated with myofibroblast subtypes. Color bars indicate expression (log₁P norm). h, PROGENy pathway scores for fibroblasts from lesional and healthy samples. Illustrations in f were partly created using BioRender.com.

Fig. 4. Fibroblast compositional signatures characterize the stroma of distinct skin diseases and scarring risk categories.

a, Proportion of fibroblast populations by individual disease. Labels overlying each bar indicate the disease category. b, Proportion of disease-adapted and disease-specific fibroblast subtypes by disease category (mean ± s.e.m.). Scarring risk group was based on clinical profiles (Methods). c, Immunofluorescence of LRRC15 (green) and ADAM12 (magenta) showing myofibroblast populations only in inflamed hidradenitis suppurativa skin (right-most) (from two representative atopic dermatitis and hidradenitis suppurativa inflamed and noninflamed samples). Scale bar, 100 µm. d, Xenium 5k data for lesional/inflamed atopic dermatitis skin, with cells colored by cell type. e, Xenium 5k data for cutaneous melanoma, with cells colored by cell type. f, Proportion of fibroblast populations by disease status for Xenium 5k data. g, Gene module scores for each disease-associated fibroblast subtype across diseases with row normalization (0–1). VE, vascular endothelium; SCLE, subacute cutaneous lupus erythematosus; DLE, discoid lupus erythematosus.

Fig. 5. Origin of skin disease-specific fibroblast subtypes.

a–c, Velocity pseudotime (a), directed PAGA overlaid on UMAP (b) and velocity kernel from CellRank2 for lesional fibroblasts (c). For further details see Methods. d, UMAP visualization of fibroblast subtypes from human skin wounds data colored by cell type (left) and MKI67 (encodes Ki-67) expression (bottom right). Proportions of fibroblast populations by time point (top right), where each bar represents a donor at a given time point. e, Schematic of predicted trajectories. Dashed arrows indicate predictions with multiple lines of evidence. Fibroblast populations are colored by the predominant scarring/fibrosis risk observed in an earlier analysis: green (prevalent in low-risk scarring stroma), orange (prevalent in scarring risk stroma and cancer), red (prevalent in established scarring/fibrotic disorders). Gray boxes indicate signaling pathways identified in our gene expression/pathway analysis. Schematic in e was partly created using BioRender.com.

Fig. 6. Human cross-tissue disease fibroblast populations.

a, Human tissues previously included in cross-tissue fibroblast studies and the fibroblast subtypes they have identified (above heatmap), colored by study (Buechler et al. in green, Korsunsky et al. in orange and Gao et al. in blue). Heatmap shows gene expression of marker genes previously reported for cross-tissue fibroblast populations in our lesional skin fibroblast subtypes. Immediately above the heatmap we show the skin fibroblast subtype with most similar gene expression to reported cross-tissue populations. b, UMAP visualization for cross-tissue integration (left) and for fibroblasts specifically (right), colored by tissue. Color bars indicate expression (log₁P norm). c, UMAP visualization for fibroblasts colored by re-annotated clusters (Methods). d, Dotplots of expression of marker genes we previously used for skin fibroblasts in cross-tissue atlas clusters by tissue type. Note that not all genes were available as the endometrial dataset contained ~17,000 genes. Illustrations in a were partly created using BioRender.com.

Fig. 7. Cross-tissue F3: FRC-like fibroblasts and F6: inflammatory myofibroblasts regulate skin immune niches.

a, Proportion of disease-associated F3: FRC-like fibroblasts, F6: inflammatory myofibroblasts and F7: myofibroblasts in non-skin tissues by disease (left) and disease category (right). b, Dotplot of F3: FRC-like expression across diseases (skin and non-skin) (left). Dotplot of F6: inflammatory myofibroblasts gene expression across diseases (skin and non-skin) (right). Diseases with a minimum of 50 cells. c, H&E slide of lesional atopic dermatitis skin with annotation of perivascular infiltrate regions (top left). Niche identification and proportion of cells in the perivascular superficial niche (bottom left). Composition of the perivascular superficial niche in 10x Genomics Xenium. Insert: zoomed in version of perivascular niche cluster. d, Cell–cell communication analysis for F6: inflammatory myofibroblasts and skin immune cells (Methods). T_CM, T central memory. e, Proportion of F6: inflammatory myofibroblasts in IBD by intestinal tissue inflammation status and linear regression of proportion with inflammation scores with 95% CI. f, Schematic summary of cell–cell interactions for F3: FRC-like fibroblasts and F6: inflammatory myofibroblasts. Schematic in f were created using BioRender.com.

Fig. 8. Adult skin F3: FRC-like fibroblasts potentially arise from prenatal skin LTo cells.

a, UMAP visualization for prenatal human skin and adult skin colored by cell type and (insert) age group. b, Heatmap of adult skin gene module scores applied to prenatal skin fibroblasts. c, UMAP visualization for prenatal human skin and intestine. Insert: cluster of LTo-like cells. (m)LTo: (mesenchymal) lymphoid tissue organizer. Dotplot of marker gene expression for LTo-like cells by tissue. d, Dotplot of expression of F3: FRC-like fibroblasts and LTo-associated marker genes in healthy adult skin F3: FRC-like fibroblasts, prenatal skin CCL19⁺ fibroblasts and intestinal mLTo cells. e, Dotplot of expression of marker genes for adult skin F3: FRC-like fibroblasts and LTo-associated marker genes in diseased adult skin, by disease, for diseases with a minimum of 100 F3: FRC-like fibroblasts. DRESS, drug reaction with eosinophilia and systemic symptoms. f, Dotplot of expression of F3: FRC-like fibroblast marker genes in mouse steady-state cross-tissue atlas (Buechler et al.). g, Proportion of Ccl19⁺ fibroblasts (labels from original study) in mouse in different healthy tissues (including mouse flank skin) and of F3: FRC-like fibroblasts in healthy human skin. Illustrations in f and g were partly created using BioRender.com.

Extended Data Fig. 1. Metadata and integration with Schwann cells, pericytes, and fascial cells.

a) Marker genes for fibroblast selection. b) Proportion of fibroblasts by site status, c) sex, d) dataset, and e) sample. f) Anatomical location latent space generated by CellDISECT colored by anatomic location and cell subtype (see Methods). Proportion of fibroblasts by anatomic location. g) Differentially expressed HOX genes and other differentially expressed genes by anatomic site using CellDISECT anatomic site clusters. h) Heatmap of reported marker genes for fibroblastic reticular cell (FRC) subsets by skin fibroblast subtype. TRC: T reticular cell. TBRC: T-B border reticular cells. FDC: Follicular dendritic cell. PRC: Perivascular reticular cell. i) Dotplot of marker genes for fascial fibroblasts and expression of F2 and F2/3 markers. j) UMAP of additional integration including Schwann cells, fascial fibroblasts, and pericytes. Inlet: UMAP color by expression of normalized PDGFRA expression (general fibroblast marker). k) Dotplot of Schwann cell, pericyte, and F5: Schwann-like fibroblast gene expression. l) Dotplot of expression of CGRP receptor components by fibroblast subtype.

Extended Data Fig. 2. Pathway analysis and transcription factor enrichment for healthy fibroblast populations.

a) Gene set enrichment analysis for healthy/nonlesional fibroblasts using GSEAPY. b) Heatmap of transcription factor activity inference.

Extended Data Fig. 3. Fibroblast location by deconvolution of fibroblast spots.

a) H&E stain of skin slides used for 10x Genomics Visium, cell type abundances predicted by cell2location, and histopathological annotations. H&E: Hematoxylin and eosin b-d) Close-up regions of H&E stained image for regions annotated in panel a.

Extended Data Fig. 4. Additional Xenium data.

a) Heatmap of marker gene expression in 10x Genomics Xenium data for fibroblast and non-fibroblast populations. b) Neighborhood enrichment scores for each tissue section shown in Fig. 2b. c) 10x Genomics Xenium section of nonlesional atopic dermatitis skin highlighting the location of F2/3: Perivascular and F5: Schwann-like fibroblasts at higher zoom. d) Xenium data colored by cell type annotation and H&E stain of skin indicating position of F4: TNN + COCH + , superior to the insertion of the arrector pili muscle.

Extended Data Fig. 5. Disease-adapted and disease-specific fibroblast additional data.

a) Dotplot of expression of marker genes used for healthy fibroblasts in lesional/diseased skin fibroblasts. b) Heatmap of genes that were used as marker genes for healthy fibroblasts but that were less specific when diseased fibroblasts were included. c) Gene set enrichment analysis for disease-associated fibroblast populations. d) Heatmap of gene expression for select collagens, common fibroblast activated genes, and additional markers for skin myofibroblasts. e) UMAP visualization of transcription factor activity scores (RUNX2, SCX).

Extended Data Fig. 6. Disease stroma profiles.

a) Fibroblast proportions in neurofibroma and dotplot of expression of F2/3 and F5 marker genes for fibroblasts in neurofibroma. b) Methodology for random forest classifier. Results: F1-score for scarring risk category classification and importance scores by fibroblast subtype. c) Heatmap of marker genes for fibroblast and non-fibroblast populations in Xenium 5k data for lesional/inflamed atopic dermatitis skin. d) Cell2location predictions and histopathological annotation for lesional (inflamed) atopic dermatitis skin profiled by 10x Visium. Same scale used for fibroblast subtypes as shown in Extended Data Fig. 3a. e,f) Zoomed in regions of perivascular inflammatory infiltrate (predicted location of F3: FRC-like fibroblasts). g) Heatmap of marker genes for fibroblast populations in Xenium 5k data for cutaneous melanoma. CAF: cancer-associated fibroblast. h) Heatmap of expression of skin CAF genes from Forsthuber et al. in our skin fibroblast subtypes (scRNA-seq data). i) Dotplot of expression of inflammatory CAF marker genes and our reported F6: Inflammatory myofibroblast marker genes in F6: Inflammatory myofibroblasts from different disease categories.

Extended Data Fig. 7. Additional trajectory inference.

a) UMAP and partition-based graph abstraction (PAGA) for all fibroblast subtypes. b) Trajectory inference from Monocle 3 (using F2: Universal as the root state) (left), pseudotime (middle), and RNA velocity from scVelo (right). c) Marker gene expression for fibroblasts identified from human wound data. d) Dotplot of MKI67 expression by time point in F6: Inflammatory myofibroblasts. e) PAGA analysis of human wound data.

Extended Data Fig. 8. Additional cross-tissue populations and human skin-intestine fibroblast integration.

a) Heatmap of marker genes for additional reported shared/universal fibroblast populations from Gao et al. b) Proportion of fibroblasts by tissue in renamed fibroblast clusters from cross-tissue atlas. c) Dotplot of skin fibroblast marker genes for cross-tissue fibroblast clusters (excluding skin fibroblasts). Skin fibroblasts were excluded to avoid skin fibroblasts driving the gene expression signature for each cluster. d) Expression of F1 marker genes by tissue in fibroblasts annotated as F1: Superficial fibroblasts in our cross-tissue atlas. e) Dotplot of F3 marker gene expression in clustered HLCA data and bar plot showing original annotations for fibroblasts that we re-annotated as ‘F3-like’. HLCA: Human Lung Cell Atlas. f) Dotplot of skin marker genes in HLCA data (original annotations from HLCA). g) UMAP visualization for joint human skin-intestine fibroblast integration using labeled data from the gut atlas (Elmentaite et al). Dotplot of marker genes for F6: Inflammatory myofibroblasts (from skin) in F6: Inflammatory myofibroblasts (skin) and Stromal 4 (MMP1+) (intestine) populations; dotplot of marker genes for F3: FRC-like (from skin) in F3: FRC-like (skin) and T reticular cells (intestine) populations (same scale for both plots); and dotplot of marker genes for F2/3: Perivascular in F2/3: Perivascular (skin) and various intestinal populations (same scale for all plots).

Extended Data Fig. 9. F3: FRC-like fibroblast cell–cell communication.

a) Cell–cell communication results from CellphoneDB, including only significant interactions and genes identified as F3: FRC-like marker genes previously. b) Dotplot of marker genes for non-fibroblast cells in inflamed atopic dermatitis skin (superficial perivascular niche).

Extended Data Fig. 10. Prenatal and cross-species F3 comparison.

a) Proportion of fibroblasts in lymph node tissue that were annotated as the Ccl19+ subtype in mouse cross-tissue fibroblast atlas. b) Top 20 differential expressed genes (DEGs) for adult and prenatal skin fibroblasts. c) Dotplot of expression of adult skin fibroblast marker genes in prenatal skin. d) Heatmap of adult skin gene module scores applied to prenatal skin fibroblasts using top 1000 genes.