scRNA-seq analysis of sun-protected human skin

The anatomy of the skin can vary considerably depending on a number of endogenous and environmental factors28,29. In addition to the dermal changes that occur as a result of intrinsic, chronological aging, there are the effects of photoaging caused by chronic exposure of the skin to low, non-extreme doses of UV radiation, the most important cause of extrinsic skin aging25,26,30. To minimize confounding effects of photoaging, our scRNA-seq analysis was based on five independent whole-skin samples that were obtained specifically from the sun-protected inguinoiliac region of male donors. Since we also sought to address the effects of intrinsic aging on dermal fibroblasts, samples were obtained from two younger (25 and 27 y/o) and three older (53, 69, and 70 y/o) donors. After enzymatically and mechanically disrupting the tissue, dead cells were thoroughly removed and samples were subjected to scRNA-seq using the 10X Genomics platform (v2 chemistry). This commercial version of the high-throughput Drop-seq protocol31 identifies cell populations by analyzing the expression of highly expressed genes in a high number of cells.

In an initial analysis, we obtained an overview of the diverse skin populations by integrating the cells from all five samples. Data analysis of the 15,457 cells that passed our quality controls (Supplementary Table 1, see Methods for details) resulted in a uniform manifold approximation and projection (UMAP)32 plot displaying 17 clusters with distinct expression profiles (Fig. 1a). Importantly, all identified clusters contained cells from all donors (Supplementary Fig. 1a). Comparing known markers with the most representative expressed genes of each cluster (Fig. 1b and Supplementary Data 1) revealed the identity of the 17 cell clusters, all of which are known constituents of the human skin and represent nine main cell types (Fig. 1c and Supplementary Figs. 1b–d). Seven clusters comprised the two key cell types of the skin, keratinocytes and fibroblasts. Keratinocytes were detected in three clusters (#5, #7 and #15) and their diversity was mainly due to their degree of differentiation. While epidermal stem cells (EpSC) and other undifferentiated progenitors (#7 and #15) expressed markers such as KRT5, KRT14, TP63, ITGA6, and ITGB1, differentiated keratinocytes (#5) were defined by KRT1, KRT10, SBSN, and KRTDAP expression33 (Fig. 1c and Supplementary Fig. 1c). Fibroblasts were identified by their archetypal markers LUM, DCN, VIM, PDGFRA, and COL1A227, constituted the most abundant skin cell type (5948 cells in total), and were represented by four clusters (#1, #2, #3, and #9, Fig. 1c and Supplementary Fig. 1c). Despite the lower amount of cells, the individual analysis of each sample generated a similar number of clusters and identified the same major cell types (Supplementary Fig. 2).

Fig. 1: Single-cell RNA sequencing analysis of sun-protected whole human skin identifies seventeen distinct cell populations. a Uniform manifold approximation and projection (UMAP) plot depicting single-cell transcriptomes from whole human skin (n = 5). Each dot represents a single cell (n = 15,457). Coloring is according to the unsupervised clustering performed by Seurat. b Heatmap showing the five most differentially expressed genes of each cell cluster, as provided by Seurat. Each column represents a single cell, each row represents an individual gene. Two marker genes per cluster are shown on the right. Yellow indicates maximum gene expression and purple indicates no expression in scaled log-normalized UMI counts. c Average expression of 3–5 well-established cell type markers was projected on the UMAP plot to identify all cell populations (see Methods for details). Red indicates maximum gene expression, while blue indicates low or no expression of a particular set of genes in log-normalized UMI counts. DC dendritic cells, EC endothelial cells. Full size image

Functional and spatial signatures of fibroblast subpopulations

To investigate whether specific functions could be assigned to the different fibroblast subpopulations, we performed gene ontology (GO) analyses using the most representative markers of each cluster. Since it is well established that skin and its fibroblasts undergo specific changes upon aging34,35,36, we only used the expression profiles of the fibroblasts from young samples (1792 cells) for this analysis (Fig. 2a and Supplementary Data 2). Classical fibroblast functions related to collagen or ECM production and organization were strongly enriched for three of the clusters (#1, #3, and #9, Fig. 2b). GO analyses also assigned mesenchymal functions such as skeletal system development, ossification or osteoblast differentiation to the cells belonging to clusters #3 or #9 (Fig. 2b). Interestingly, our results also showed a strong enrichment for functions related to inflammation specifically in the fibroblasts of cluster #2. Significant examples include inflammatory response, cell chemotaxis or negative regulation of cell proliferation, necessary for the final anchoring of leukocytes (Fig. 2b). Functions that are typically attributed to fibroblasts, such as collagen or ECM production and organization, did not appear among the most statistically significant categories for the cells of this cluster. These findings provide a first illustration for the functional heterogeneity of the fibroblasts in our samples.

Fig. 2: Dermal fibroblast subpopulations display specific spatial and functional transcriptomic signatures. a Left: UMAP plot displaying dermal fibroblasts from young donors (n = 2). Each dot represents a single cell (n = 1792). Coloring is according to the original unsupervised clustering performed by Seurat. Right: Bar plots indicating the percentage of fibroblasts corresponding to each subpopulation and donor. b Top 8 enriched Gene Ontology (GO) terms in each fibroblast subpopulation, sorted by p-value. c Heatmap showing the expression of all collagen genes in the distinct fibroblast subpopulations. Each column represents a single cell and each row an individual collagen gene. Yellow indicates maximum gene expression while purple indicates no expression in scaled log-normalized UMI counts. d Average expression of the genes constituting the papillary and reticular gene signatures for predicting dermal localization of the fibroblasts from the four clusters. In all UMAP gene expression projections, red indicates maximum expression and blue indicates low or no expression of each particular set of genes in log-normalized UMI counts. In the violin plots, X-axes depict cell cluster number and Y-axes represent average expression of each set of genes in log-normalized UMI counts. Statistical significance of the expression changes in the gene signatures between cell clusters is indicated below through the p-values of the corresponding Wilcoxon rank sum tests. Full size image

The expression of some specific collagens has also been linked to particular fibroblast functions37. We therefore analyzed the four fibroblast clusters at the level of collagen expression patterns. In agreement with our GO analysis, the results clearly indicated lower global collagen expression levels in cluster #2 (Fig. 2c). The analysis also revealed that the collagen genes COL11A1 and COL24A1, associated with cartilage and bone development38,39, respectively, were specifically expressed in the fibroblasts of cluster #9 (Fig. 2c), suggesting a stronger mesenchymal component for this cell subpopulation. For the remaining collagen and ECM secreting clusters (#1 and #3), the data showed a bias in the production of collagens that have previously been linked to specific dermal locations (Fig. 2c). In particular, fibroblasts in cluster #3 express COL13A1 and COL23A1, two known markers of papillary fibroblasts27,40,41,42. High expression levels of another epidermal-dermal junction collagen gene, COL18A1, also supported the location of these fibroblasts within the papillary layer27,40,43,44.

To better predict the potential localization of the four observed fibroblast subpopulations within the dermis, we next studied the expression of sets of genes that have previously been related to papillary or reticular fibroblasts. While the most representative markers of the papillary fibroblasts comprise APCDD1, AXIN2, COLEC12, PTGDS, and COL18A1, the reticular fibroblast signature is typically defined by a group of ten genes, including MGP or MFAP511,27,40,44 (Supplementary Fig. 3). In agreement with the collagen expression data, our results showed that the papillary gene expression signature is mostly restricted to fibroblasts in cluster #3 (Fig. 2d). In contrast, fibroblasts in clusters #1, #2, and #9 showed more prominent expression of the reticular signature, with the highest reticular (and lowest papillary) expression levels observed in cluster #1 (Fig. 2d).

Taken together, these results suggest that the functional annotation differences between subpopulations reflect their priming for distinct functional roles. Our results thus predict two subpopulations with prominent roles in the generation of structural collagen and ECM organization, one located in the reticular dermis (cluster #1, ‘secretory-reticular fibroblasts’) and the other in the papillary dermis (cluster #3, ‘secretory-papillary fibroblasts’). A third subpopulation, which was predicted to have a more reticular localization, showed a greater mesenchymal potential (cluster #9, ‘mesenchymal fibroblasts’). The fourth and final subpopulation, also predicted to have a mostly reticular localization, was characterized by pro-inflammatory functions (cluster #2, ‘pro-inflammatory fibroblasts’).

The subdivision of fibroblasts into four subpopulations was also confirmed by a second-level clustering of the 1792 fibroblasts obtained from the young samples. This approach also identified both secretory and mesenchymal subpopulations, while further subdividing the pro-inflammatory subpopulation into two closely related subclusters that were defined by the differential expression of a subset of genes (Supplementary Fig. 4a and Supplementary Data 3). Second-level analyses of the fibroblasts from each individual sample (young or old) also often subdivided the pro-inflammatory subpopulation and additionally separated the well-established dermal papilla-associated fibroblasts, which are characterized by high CRABP1 and TNN expression levels9,24, from the mesenchymal subpopulation (Supplementary Figs. 4b, c). However, these additional subclusters are clearly related to the four main fibroblast subpopulations and were therefore not considered separately.

Validation of fibroblast subpopulations in skin sections

To further characterize and validate the four fibroblast subpopulations from our initial analysis, we identified the most representative markers for each subpopulation according to their expression in the specific cell clusters (Table 1 and Supplementary Fig. 4a). Since no cell-surface markers were found specific enough for all subpopulations, to assess the microanatomical distribution of the characterized subpopulations we then performed RNA FISH on independent, formalin-fixed paraffin-embedded (FFPE) skin sections from young (28–37 y/o) and old (54-86 y/o) donors. Collagen Triple Helix Repeat-Containing 1 (CTHRC1), and Adenomatosis polyposis coli downregulated 1 (APCDD1) were selected to detect the secretory-reticular and secretory-papillary fibroblasts, respectively. C-C motif chemokine ligand 19 (CCL19) and Apolipoprotein E (APOE), were used to detect the pro-inflammatory fibroblasts, and Asporin (ASPN) was chosen as a marker for the mesenchymal subpopulation. Platelet-derived growth factor receptor alpha (PDGFRA) was included as a pan-fibroblast control45. RNA FISH experiments confirmed the location of each secretory subpopulation within the papillary and the reticular dermal layers, respectively (Fig. 3a and Supplementary Fig. 5a). The locations were also confirmed by immunofluorescence staining of Tetraspanin 8 and the Collagen alpha-1(XVIII) chain, two additional markers of the secretory subpopulations (Table 1, Supplementary Figs. 4a and 6a, b). The pro-inflammatory fibroblasts showed a more widespread distribution and a preferential association with the vasculature (Fig. 3b and Supplementary Fig. 5b). Finally, the mesenchymal subpopulation was localized mostly to the reticular dermis, particularly in the vicinity of hair follicles (Fig. 3c and Supplementary Fig. 5c). This localization was also confirmed by immunofluorescence staining of Periostin, which is another marker for this subpopulation (Table 1, Supplementary Figs. 4a and 6c). Together, these results provide important confirmation for our findings obtained by single-cell transcriptomics and establish markers for the detection of specific fibroblast subpopulations.

Table 1 Representative marker genes of each fibroblast subpopulation. Full size table

Fig. 3: RNA FISH detection of fibroblast subpopulations in young skin. a Representative confocal images showing mRNA expression of CTHRC1 (green) and APCDD1 (red), selected markers for the secretory-reticular and secretory-papillary fibroblast subpopulations, respectively. Details from the papillary and reticular regions of the images above are shown in the lower panels (left and center, respectively), and percentage of positive cells for each gene and per dermal region are shown in the lower right panel. b Representative confocal images showing mRNA expression of CCL19 (green) and APOE (red), selected markers for the pro-inflammatory fibroblast subpopulation. A detail of a vessel of the images above is shown in the lower panel. c Representative confocal images showing mRNA expression of ASPN (green), selected marker for the mesenchymal fibroblast subpopulation. A detail of the hair follicle bulb of the images above is shown in the lower panel. Dashed lines in a and b denote the papillary dermis regions while in c denote the dermal papilla. Nuclei were counterstained with DAPI. Each assay was performed in three independent young FFPE skin sections (28–37 y/o). Images are shown at ×40 original magnification. Scale bar: 50 μm for main images and 10 μm for detail images. Pap papillary dermis, Ret reticular dermis, Deep ret deep reticular dermis, HF hair follicle, DP dermal papilla. Statistical analyses were performed using a two-way ANOVA test (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001); error bars represent the standard deviation. Full size image

Aging leads to loss of dermal fibroblast priming

We also investigated the effect(s) of aging at the level of dermal fibroblast subpopulations. In relative terms, our results suggest an apparent reduction in the number of mesenchymal fibroblasts in the old samples (Figs. 2a and 4a). However, this observation could not be experimentally demonstrated due to the low amount of hair follicles present in the available tissue sections. Next, we compared the fibroblast transcriptomes from old donors with their young counterparts. In agreement with the reduced proliferative capacity of aged cells, the expression profiles from old fibroblasts indicate a significant delay at the G1/S transition of the cell cycle in the pro-inflammatory and secretory-papillary subpopulations, as well as a similar tendency for the secretory-reticular cells (Fig. 4b). Interestingly, GO analyses of the most representative genes from the old subpopulations (Supplementary Data 4) suggested a considerable age-dependent loss of the functional annotations for each cluster. In comparison to young fibroblast subpopulations, the aged counterparts showed fewer function-related terms, and/or substantially reduced p-values, consistent with fewer genes supporting the terms (Fig. 4c). A similar effect was also seen at the level of collagen gene transcription, which became particularly decreased in the secretory clusters (Fig. 4d and Supplementary Figs. 7 and 8a). Finally, aging also changed the previously defined spatial gene expression signatures, as old papillary fibroblasts presented less papillary and more reticular gene expression signatures, while reticular fibroblasts presented a less pronounced reticular gene expression signature (Fig. 4e and Supplementary Fig. 8b). These findings are in agreement with an age-related loss of fibroblast priming.

Fig. 4: Aging leads to loss of dermal fibroblast priming. a Left: UMAP plot displaying dermal fibroblasts from old donors (n = 3). Each dot represents a single cell (n = 4156). Coloring is according to the original unsupervised clustering performed by Seurat. Right: Bar plots indicate the percentage of fibroblasts corresponding to each subpopulation and donor. b Percentage of fibroblasts of each subpopulation that were in the G1, S or G2/M phase of the cell cycle in young and old skin samples, respectively. c Top 8 enriched Gene Ontology (GO) terms in each old fibroblast subpopulation sorted by p-value. Coloring is according to the unsupervised clustering performed by Seurat. d UMAP and violin plots displaying the average expression of all collagen genes in the fibroblasts of all subpopulations, for young and old skin. e UMAP and violin plots displaying the expression of the papillary and reticular gene signatures in the fibroblasts of all subpopulations, for young and old skin. In all UMAP gene expression projections, red indicates maximum expression and blue indicates low or no expression of each particular set of genes in log-normalized UMI counts. In the violin plots, X-axes depict fibroblast subpopulations and Y-axes represent average expression of each set of genes in log-normalized UMI counts. For comparing the ratio of G1 cells between young and old subpopulations, a two-sided two-proportion z-test was used (b). Statistical analyses in d and e were performed using the Wilcoxon Rank Sum test. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. Y young, O old, S.R Secretory-reticular, INF Pro-inflammatory, S.P Secretory-papillary, MES Mesenchymal. Full size image

Aging effects on SAASP profiles and cell–cell interactions

Finally, we also analyzed whether age-related transcriptomic changes in fibroblast subpopulations could explain age-related skin phenotypes. For example, it is known that aged fibroblasts become more susceptible to the accumulation of reactive oxygen species46. Consistently, GO analyses performed with the most downregulated genes of each aged cluster (Supplementary Data 5) showed that genes related to hydrogen peroxide metabolism were decreased in three of the four fibroblast subpopulations (Supplementary Fig. 9). Furthermore, old skin is known to acquire a chronic, low-grade inflammatory phenotype47,48. This could also be detected at the level of fibroblast subpopulations, as we found immune response among the main enriched terms in our GO analyses of the most up-regulated genes of each aged cluster (Supplementary Figs. 9 and 10 and Supplementary Data 5). Old fibroblasts also showed changes in the expression of genes that were previously detected in the aging dermal fibroblast secretome49. Furthermore, differential expression profiles of SAASP genes could also be observed in the different old fibroblast subpopulations (Fig. 5a and Supplementary Fig. 11). In conclusion, our scRNA-seq data thus recapitulates known phenotypes associated with skin aging.

Fig. 5: Other age-related changes in dermal fibroblast subpopulations. a Expression of genes encoding skin aging-associated secreted proteins (SAASP) (rows) that are differentially (fold-change > 1.25) expressed between young and old fibroblasts in at least one subpopulation (columns). The heatmap shows the mean relative expression by cluster. b Bar plots showing the number of ligand-receptor interactions predicted for the four observed fibroblast subpopulations with the rest of the cell types identified in human skin, in the two young (≤27 y/o) (up) and the two oldest (≥69 y/o) (down) samples. The medium-aged sampled (53 y/o) showed an intermediate phenotype in this analysis and was therefore omitted. Coloring and numbering are according to the original unsupervised clustering performed by Seurat. c Summary of the top four exclusive interactions lost between each fibroblast subpopulation and undifferentiated keratinocytes, sorted by p-value. The table shows interactions in both directions for each pair. Y young, O old, S.R Secretory-reticular, INF Pro-inflammatory, S.P Secretory-papillary, MES Mesenchymal. Full size image

Finally, fibroblasts are known to establish interactions with many other skin cell types during homeostasis6. Importantly, scRNA-seq also provides novel opportunities to identify communicating pairs of cells based on the expression of cell-surface receptors and their interacting ligands50. Interestingly, our results indicate that a high number of the interactions predicted for young fibroblasts are lost during intrinsic skin aging. This effect was particularly pronounced in the two oldest samples (≥69 y/o) and for interactions involving undifferentiated keratinocytes (Fig. 5b, c and Supplementary Fig. 12). These findings suggest that the loss of interactions between fibroblasts and their communicating cell types represent a previously unrecognized molecular phenotype of the aging human skin.