We obtained control punch skin biopsies from non-diseased skin of two patients undergoing skin surgery at the Department of Dermatology, University Hospital Zurich, Switzerland. Additionally, single punch skin biopsies were collected from six patients with systemic sclerosis (SSc) who were admitted to the Department of Rheumatology, University Hospital Zurich, Switzerland for diagnostic procedures. Clinical characteristics of patients are provided in Supplementary Table S1. All patients signed the informed consent for participation in the study prior to study enrollment. The ethical committee of the Canton Zürich (ethical approvals No. EK-687, EK-800, PB-2020-00035) approved the skin collection and analysis of patients’ data and skin tissues.

For each control subject, we collected two paired 4 mm punch biopsies which represented technical duplicates for tissue digest. These biopsies were placed in complete RPMI (10% FCS) medium, shipped at 4°C and freshly dissociated into a single cell suspension for scRNA-seq within 2 h of skin collection (see step by step protocol). Punch skin biopsies (4 mm) from patients with SSc (n = 6) were collected from nonfibrotic skin of the proximal dorsal forearm under local anesthesia with 1% Lidocaine. Biopsies were placed into complete RPMI (10% FCS) medium, shipped at 4°C, cut into small fragments with a scalpel (see step by step protocol) within 2 h of skin collection and cultured ex vivo for 24 h. Explanted skin fragments were then digested (n = 6) into single skin cell suspension (see step by step protocol) and used for scRNA-seq (n = 4).

scRNA-seq libraries from skin cells were prepared using the single-cell 3′ v3.1 protocols and reagents (10× Genomics) according to the manufacturer’s instructions. 6,000 single skin cells were targeted for the encapsulation. The quality and quantity of cDNA and scRNA-seq libraries were assessed with Agilent Bioanalyzer (Agilent Technologies). Diluted 10 nM libraries were pooled in equimolar ratios, and the library pool was sequenced on the Illumina NovaSeq6000 sequencer (paired-end reads, R1 = 28, i7 = 8, R2 = 91, sequencing depth min. 50,000 reads/cell).

Cell Ranger (v2.0.2, 10× Genomics) was used to demultiplex, align the reads to Ensemble reference build GRCh38.p13 and collapse unique molecular identifiers (UMIs). The fastq files were processed into transcript count tables using Cell Ranger count (version 6.0.0) with the reference Genome GRCh38.p13 from the Genome Reference Consortium (Schneider et al., 2017). The subsequent analysis was carried out on the filtered Cell Ranger count matrix in which empty droplets had already been removed. This was done using R (version 4.1.0) and Bioconductor (version 3.13) packages. Further potential doublets were identified in the R analysis using scDblFinder (Germain et al., 2021). We ran the same processing pipeline on the fresh skin scRNA-seq data and the combined fresh and cultured scRNA-seq datasets.

The data was first filtered for cells with at least 150 expressed genes and for genes expressed in more than 10 cells. Quality control (QC) was applied based on 1) library size, 2) number of genes, 3) mitochondrial gene counts as a percentage of the library size of a cell. Cells were excluded if the library size or number of genes deviated by two median absolute deviations (MADs) below the sample median for 1) and 2), and above two MADs for 3). Furthermore, explicit filters were applied for a library size of at least 1,500 counts, a gene count of at least 500, and cells where no more than 25% of the library size was made up by mitochondrial counts. The data was normalized by library size and subsequently the top 2,000 highly variable genes were identified for dimensionality reduction. The data was batch corrected using the Harmony package (Korsunsky et al., 2019) and clustered using the Walktrap algorithm in the iGraph and bluster package (Lun, 2021). Marker gene analysis was run on the resulting clusters using Scran (Lun et al., 2016) to annotate clusters and identify cell types. As a reference for subsequent cluster annotation, cells were annotated using AUCell (Aibar et al., 2017) with skin cell markers identified by Xue et al. (2020), He et al. (2020), Tabib et al. (2021), Raynolds et al. (2021) as well as our selected marker genes for expected cell types. Cell types were then identified by annotation of the clusters. In the fresh skin dataset one of the clusters was found to contain Melanocytes and Schwann cells. Cells in this cluster were assigned to the respective cell types using a threshold of the expression of a marker gene which was one of several identified in iSEE to be capable of clearly distinguishing between the two cell types. This resulted in ten cell populations identified in the fresh samples and twelve cell populations identified in the fresh and cultured samples combined. The identified major cell populations in the fresh skin where then further subclustered and annotated for cell subtypes. For the subclustering of some cell types, the hyperparameters in the clustering method were varied to cluster cells which were identifiable as clear subgroups in feature plots.

We integrated our scRNA-seq data from freshly dissociated skin with skin scRNA-seq datasets from He et al. (2020), Solé-Boldo et al. (2020), and Tabib et al. (2021) using a similar pipeline as described above with some minor deviations. Only the healthy skin tissues from these datasets were included in our integration analysis. Authorized access to the data from He et al. (2020), Solé-Boldo et al. (2020), and Tabib et al. (2021) was obtained via the linked data on the Gene Expression Omnibus (GEO) accession numbers linked to the respective papers (GSE147424, GSE138669, GSE130973). As only the raw reads were available for Tabib et al., we started the data analysis with the raw Cell Ranger counts for all samples to allow a fair comparison. The four scRNA-seq datasets were subset to a list of 15,190 genes that were found to occur in at least 10 cells of each protocol and then merged to one dataset. For the QC, the MAD cut-offs for the relative filters were relaxed to 3 MADs for filters 1 to 3 (see above). Furthermore, the manually set filters for the minimum library size and gene count were lowered from 1,500 to 750 and from 750 to 500, respectively, to allow inclusion of sufficient cell numbers from all datasets. As the datasets exhibited strong heterogeneity, scTransform (Hafemeister and Satija, 2019) was chosen as a more robust normalisation method. Highly variable genes were identified by the deviance of their variance from the modelled trend using Scran (Lun et al., 2016), using a false discovery rate of 5% resulting in 1,290 genes being identified. Due to the size of the integrated dataset, the clustering method was adapted to a more efficient two-step clustering using K-means for centroid identification first and then Walktrap clustering subsequently. Ten clusters were identified which were manually annotated after analysis of the marker genes as mentioned above. Mast cells were not identified by the unsupervised clustering method but were apparent when inspecting marker genes as well as feature plots and were manually annotated based on these marker genes. Additionally, to allow for a fair comparison of detected genes across samples, we downsampled our and published datasets in a separate analysis to the sequencing depth of the dataset with the lowest coverage (Supplementary Figure S5). All datasets were downsampled using the downsampleBatches function from the scuttle R package, such that the average per-cell total count is the same across batches, i.e., in this case the samples. This analysis does not just shrink the counts by a factor across cells and samples but instead, it randomly samples counts from the data to leave all cells with a pre-specified count—in this case of the sample with the lowest coverage—thereby simulating a stronger dropout effect for samples with comparatively higher coverage.

In this step, we describe pre-processing steps to prepare fresh and cultured punch skin biopsies for subsequent tissues digestion.

Timing: 10 min

Pre-processing fresh skin

1) Start processing fresh skin, stored in complete RPMI (10% FCS) medium, within 2 h of skin biopsy acquisition.

Skin explant culture

2) To prepare skin biopsy for explant culture, proceed with skin processing within 2 h of biopsy acquision.

In this step, we dissociate fresh or cultured skin tissue fragments into a suspension of single cells and partially digested tissue fragments using the combined enzymatic-mechanical tissue dissociation protocol involving dispase II. Dispase is a protease that hydrolyzes the N-terminal peptide bonds of non-polar amino acid residues (Volovitz et al., 2016). It cleaves extracellular matrix components such as fibronectin, collagen IV and to a lesser extent collagen I and is thus used to gently dissociate skin tissue.

Timing: 150 min

1) Prepare dispase digestion mix in a 15 ml conical tube.

In this step, dispase pre-digested skin tissue fragments are further dissociated into a suspension of single cells using the combined enzymatic-mechanical tissue dissociation protocol which includes collagenase IV. Collagenase IV is designed to be especially low in tryptic activity to limit damage to membrane proteins and receptors. It cleaves the bond between a neutral amino acid (X) and glycin, abundant in collagen, thus it is able to break the peptide bonds in collagen, facilitating extracellular matrix digestion and releasing cells into suspension (Seltzer et al., 1990).

Timing: 80 min

1) Prepare collagenase digestion mix in a conical 15 ml polystyrene tube.

In this step, the remaining tissue fragments are additionally digested with trypsin, to release any remaining cells and maximize cell yield. The tissue debris and cell aggregates are removed. In our experience this step contributed to further cell enrichment increasing cell counts from 16% to 44% with median cell viability of 88%.

Timing: 30 min

1) Pre-warm 3 ml/well of 0.25% trypsin in a 6-well plate kept in cell culture incubator at 37°C.

In this step, cell yield and viability are determined, potential cell debris and aggregates removed by additional filtering, and cell suspension are diluted or concentrated to be ready for starting the Gel Bead-in-Emulsion (GEM) generation.

Timing: 15 min

1) Filter the combined single cell suspension through the 35 μm strainer cap of Falcon 5 ml tubes (Catalog number #100-0087) into a 1.5 ml Eppendorf tube, repeat filtering two times.

To generate a single cell atlas of the fresh human skin, we analyzed scRNA-seq data from freshly dissociated skin samples. Specifically, 6,192 freshly dissociated cells passed QC (see Methods). Detailed QC analysis of scRNA-seq datasets from freshly dissociated samples showed a similar distribution of library sizes and identified genes per cell with average of 2,767 and 2,538 genes detected per sample, respectively. Furthermore, both samples exhibited a sufficiently high number of viable cells and did not exhibit any aberrant behavior. Filtered high quality cells, that passed QC, segregated into ten major clusters (Figure 2A), which displayed distinct expression profile of cell marker genes (Figure 2B). These cell clusters represented ten main cell types including COL1A1+ DCN+ fibroblasts, KRT5+/KRT10+ keratinocytes, TAGLN+ pericytes/vascular smooth muscle cells (VSMCs), PECAM1+ vascular endothelial cells, LYVE1+ lymphatic endothelial cells, LYZ+ macrophages/dendritic cells (DCs), CXCR4+ T cells, TPSB2+ mast cells, PMEL+ melanocytes and CDH19+ Schwann cells (Figures 2B,C). Macrophages (cluster 1, Figure 2D) and DC (cluster 2, Figure 2D) were distinguished by elevated expression of previously described markers of dermal macrophages (CD14, CD68, CD163) (Fuentes-Duculan et al., 2010), while DCs expressed CCR7 and CCL22, which are markers of mature DCs.

Fibroblasts, keratinocytes and pericytes/VSMC were among the largest resident skin cell populations (Figure 2C) and were further subclustered to identify their subsets.

By subclustering the fibroblasts in our fresh skin cell atlas, we identified five distinct fibroblast subpopulations (Figure 3A, Supplementary Figure S2A). Dermal fibroblasts show large diversity; for example, different body sites and skin compartments, such as papillary and reticular dermis, are populated by transcriptionally and functionally specialized fibroblast subsets (Janson et al., 2012). To determine whether our five fibroblast clusters represent spatially-specialized dermal fibroblast populations, we analyzed the expression of known markers of papillary and reticular dermal fibroblasts across clusters (Supplementary Figure S3). A papillary gene expression signature (Figure 3B) was mostly restricted to fibroblast clusters 3, 4, and 5, with the highest expression in clusters 3 (FOS+ JUN+ fibroblasts) and 4 (HMOX1+ TNFAIP6+ fibroblasts). In contrast, cluster 1 showed a prominent reticular gene signature, including MGP, ANGTPL1, FGF7, EFEMP1, MFAP5, suggesting that these cells represent secretory reticular dermal fibroblasts. Therefore, our subclustering analysis confirmed the presence of main spatially-specialized fibroblast subtypes with an additional distinct fibroblast cluster (cluster 2), defined by a mixed papillary/-reticular signature and a prominent inflammatory gene expression signature (CCL19, C3) (Figure 3A). Looking closer in papillary dermal fibroblast subsets, we uncovered enriched transcriptional regulation signature in cluster 3 and anti-oxidative/immune-regulatory gene signature in fibroblast cluster 4. Specifically, cluster 4 fibroblasts expressed genes involved in anti-oxidative response (HMOX1, SOD2, FTH, SQSTM) and immunomodulatory (SQSTM1, TNFAIP6, TNFAIP3, NFKB1) cell functions, while cluster 3 fibroblasts (FOS/JUN^high) resembled pro-inflammatory cluster 2 fibroblasts, but with higher collagen gene expression, confirming their papillary characteristics.

Further, we identified four distinct keratinocyte clusters that corresponded to keratinocyte differentiation states and epidermal location (Figure 3C, Supplementary Figure S2B). Keratinocytes in cluster 1 were enriched for known basal keratinocyte markers (KRT5, KRT14, KRT15), while cluster 3 keratinocytes, characterized by KRT1 and KRT10 expression, showed a gene signature of suprabasal spinous keratinocytes. KRT1 and KRT10 were also expressed in cluster 4 cells with a signature of differentiated granular keratinocytes (KRTDAP+). An additional keratinocyte population, which clustered separately (cluster 2), represented a transitional state between basal and suprabasal keratinocytes. A differentiated granular keratinocyte cluster (cluster 4) expressed the known differentiation gene markers (KRT2, IVL, SBSN, CALML5, DSC1) and contained a small subset of LOR+ FLG+ cells which could represent terminally differentiated cells of the cornified envelope (Hughes et al., 2020).

Pericyte/VSMCs (Figure 3D, Supplementary Figure S2C) separated into two different clusters. Based on the expression of established lineage markers, we identified one cluster of VMSC expressing ACTA2, TAGLN, MYH11, CNN1, TPM2 cells and one cluster of PDGFRB+ pericytes.