Author response: A single-cell transcriptome atlas of pig skin characterizes anatomical positional heterogeneity

Abstract

Full text Figures and data Side by side Abstract Editor's evaluation Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Different anatomical locations of the body skin show differences in their gene expression patterns depending on different origins, and the inherent heterogeneous information can be maintained in adults. However, highly resolvable cellular specialization is less well characterized in different anatomical regions of the skin. Pig is regarded as an excellent model animal for human skin research in view of its similar physiology to human. In this study, single-cell RNA sequencing was performed on pig skin tissues from six different anatomical regions of Chenghua (CH) pigs, with a superior skin thickness trait, and the back site of large white (LW) pigs. We obtained 233,715 cells, representing seven cell types, among which we primarily characterized the heterogeneity of the top three cell types, including smooth muscle cells (SMCs), endothelial cells (ECs), and fibroblasts (FBs). Then, we further identified several subtypes of SMCs, ECs, and FBs, and discovered the expression patterns of site-specific genes involved in some important pathways such as the immune response and extracellular matrix (ECM) synthesis in different anatomical regions. By comparing differentially expressed genes of skin FBs among different anatomical regions, we considered TNN, COL11A1, and INHBA as candidate genes for facilitating ECM accumulation. These findings of heterogeneity in the main skin cell types from different anatomical sites will contribute to a better understanding of inherent skin information and place the potential focus on skin generation, transmission, and transplantation, paving the foundation for human skin priming. Editor's evaluation This valuable manuscript provides a single-cell RNA sequencing analysis of adult pig skin from different species and anatomical regions. The evidence supporting the conclusions is compelling, with identification of molecular and cellular differences in pig skin, including analysis of regional and species-specific gene signatures. https://doi.org/10.7554/eLife.86504.sa0 Decision letter Reviews on Sciety eLife's review process Introduction The issue of how inherent information contributes to anatomical site-specific differences has inspired extensive exploration, and the pattern formation of spatial arrangement depends on the expression control of specific genes with a cell type (Rinn et al., 2006). The cellular specialization of anatomical site-specific pattern is determined in the embryo, and the inherent information could also be maintained throughout adulthood along with continual self-renewal tissues (Rinn et al., 2006). Some cellular special information of anatomical site-specific patterns in other tissues has been discovered, such as in the heart (Litviňuková et al., 2020) and muscle (De Micheli et al., 2020), but the inherent information of cellular specialization is less well understood in physiologically different anatomical skin regions. Skin is the largest organ, providing a physical, chemical, and biological barrier for the body. It consists of the upper epidermis and the lower dermis layers separated by the basement membrane, with unambiguous spatial patterns of morphological and functional specialization (Simpson et al., 2011). Embryological studies have shown that anatomical positional-specific pattern is provided by the stroma, which is composed of extracellular matrix (ECM) and mesenchymal or dermal cells during embryogenesis (Rinn et al., 2006). Pioneering studies discovered that the different anatomical regions of the body skin dermis arose from different origins. The dorsum dermis originates from the dermato-myotome, while the ventral and face dermis derive from the lateral plate mesoderm and the neural crest, respectively (Jinno et al., 2010; Ohtola et al., 2008; Wong et al., 2006). In adults, dermal cells confer the positional identity and memory for skin patterning and function (Driskell and Watt, 2015), raising the issue of what regional discrepancy could be maintained against plentiful cellular turnover in skin. The dermis is mainly composed of resident fibroblasts (FBs), smooth muscle cells (SMCs), endothelial cells (ECs), and immune cells, and these skin cells provide structure, strength, flexibility, and defense to the skin (Driskell and Watt, 2015). FBs, the main cell type in the dermis, are responsible for the collagen deposits and elastic fiber formation of the ECM (Parsonage et al., 2005), participating in skin morphogenesis, homeostasis, and various physiological and pathological mechanisms, including skin development, aging, healing, and fibrosis (Auxenfans et al., 2009; Driskell et al., 2013; Driskell and Watt, 2015). SMCs, which form blood vessels and arrector pili muscle (APM), play a critical role in controlling blood distribution as well as maintaining the structural integrity of the blood vessels and APM in skin (Driskell et al., 2013; Liu and Gomez, 2019). ECs organize the vascular plexus, which plays a predominant role in vascular remodeling, metabolism, and the immune response in the dermis, and EC metabolism is tightly connected to barrier integrity, immune and cellular crosstalk with SMCs (Cantelmo et al., 2016; Miyagawa et al., 2019; Tombor et al., 2021). Moreover, during formation and development of the skin, the dermal reaction is realized by cell–cell communication, and dynamic interactions between cells and ECM, as well as regulatory factors (Driskell and Watt, 2015). The pig is used as a model animal to research human skin biology because of its similar pathological and physiological skin attributes to those of human skin (Khiao In et al., 2019). The Chenghua (CH) pig, a novel Chinese indigenous population with superior skin thickness and strong FBs activity (Zou et al., 2022; Zou et al., 2023), is considered as a potential model animal for researching mammalian skin biology. Here, to reveal the anatomical positional heterogeneity of the skin, single-cell RNA sequencing was performed on pig skin tissues from six anatomical regions of CH pigs and the back site of large white (LW) pigs. We obtained a well-resolved single-cell transcriptome atlas of 233,715 cells and identified seven cell types with unique gene expression signatures. In our datasets, we focused on the top three cell types, including SMCs, ECs, and FBs. SMCs revealed the signature of contractile SMCs, mesenchymal-like phenotype, and macrophage-like phenotype and presented the expression patterns of site-specific genes related to ECM-integrins and immune response in different skin anatomical sites. ECs were classified into four EC phenotypes, and the gene expression patterns, which are related to integrins, immunity, and metabolism, were explored across different skin anatomical sites. Moreover, based on these comparative differentially expressed genes (DEGs) of FBs, we identified three subtypes among different regions and found that TNN, COL11A1, and INHBA might be candidate genes for ECM accumulation. Taken together, the present data offer a comprehensive understanding of the single-cell atlas that displays the cellular inherent information of anatomical site-specific patterns in skin, supporting future exploration as a baseline for healthy and morbid human skin. Results Single-cell transcriptome profiling identified different skin anatomical sites in CH pig To characterize an overview of single-cell transcriptomic atlas of pig skin tissues from different anatomical sites, we sampled CH pig skin tissues on the head, ear, shoulder, back, abdomen, and leg from three female 180-day-old individuals, then applied scRNA-seq of 18 isolated skin cell samples (Figure 1A). After stringent cell filtration, we obtained a total of 215,274 cells, which were globally visualized with 21 cell clusters in the t-SNE plot (Figure 1B). On average, 956 genes and 2687 unique molecular identifiers (UMIs) per cell were detected (Figure 1—figure supplement 1A and B). Twenty-one cell clusters were identified according to the expression matrix of marker genes for each cluster and were shown in the heatmap basing on the top 12 marker genes for each cluster (Figure 1C). The 21 cell clusters constituted seven cell types, of which the SMCs (clusters 0, 2, 5, 6, and 13) were marked by MYH11 and ACTA2, ECs (clusters 3, 4, 7,10, and 11) were marked by PECAM1 and APOA1, FBs (clusters 1, 8, 9, and 12) were expressed by LUM and POSTN, myeloid dendritic cells (MDCs) (clusters 14, 16, and 18) were labeled by BCL2A1 and CXCL8, T cells (TCs) (cluster 15) were highly expressed by RHOH and SAMSN1, keratinocytes (KEs) (cluster 17) were tabbed by KRT5 and S100A2, and epidermal stem cells (ESCs) (clusters 19 and 20) were stamped by TOP2A and EGFL8 (Figure 1D and E and Figure 1—figure supplement 1C). Figure 1 with 3 supplements see all Download asset Open asset Single-cell transcriptome profiling of different skin anatomical sites in CH pig. (A) Flowchart overview of skin single-cell RNA sequencing from different anatomical skin regions of CH pigs. (B) The t-SNE plot visualization showing 21 clusters of annotated cell types from CH pig skin. (C) Heatmap showing the top 12 highly expressed genes from each cluster. Each column represents a cluster, each row represents a gene. Light yellow shows the maximum expression level of genes, and deep green shows no expression. (D) Dot plot showing the two representative genes for each cell type. Color indicates the log2 value, and circle size indicates gene expression level. (E) The marker genes for each cell type are distributed on the t-SNE plot. Color indicates gene expression. (F) The distribution ratio of cell types for total cells and six different anatomical skin areas. (G) The most enriched Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for each cell type. CH, Chenghua; SMC, smooth muscle cell; EC, endothelial cell; FB, fibroblast; MDC, myeloid dendritic cell; TC, T cell; KE, keratinocyte; and ESC, epidermal stem cell. Figure 1—source data 1 Source data of marker genes for each cluster in Figure 1C. https://cdn.elifesciences.org/articles/86504/elife-86504-fig1-data1-v1.xlsx Download elife-86504-fig1-data1-v1.xlsx The distribution ratio of these cell types was visualized among total skin data consisting of 42.9% SMCs, 28.1% ECs, 24.6% FBs, 2.5% MDCs, 0.9% TCs, 0.6% KEs, and 0.3% ESCs, which were similar to the distribution ratio for the main cell types in different skin regions (Figure 1F). In addition, the cell number and types among the six anatomical skin sites were comparable, which indicated that the cell types displayed subtle differences, but cell number per cell type varied significantly (Figure 1—figure supplement 2). The marker genes for each cell type revealed the dominant transcriptional features and enriched pathways relevant to their distinct physical functions. Significant examples of Gene Ontology (GO) function terms were involved in ECM structural constituent or collagen binding for FBs, actin binding or structural constituent of muscle for SMCs, and ECM structural constituent or collagen binding for ECs (Figure 1G). Meanwhile, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were prominently attributed to FBs such as protein digestion and absorption or ECM–receptor interaction, ECs involved in cell adhesion molecules or the Rap1 signaling pathway, and SMCs including the NF-kappa B signaling pathway or the TNF signaling pathway (Figure 1G). Moreover, given the potential cross-species comparisons, we implemented overlapping skin cell atlases among pig, human, and mouse using a t-SNE plot (Figure 1—figure supplement 3A). The captured gene and UMI counts were more advantageous for human skin cells (Figure 1—figure supplement 3B). The cell types were similar for the three species, while the percentage of cell types was different such as SMCs, ECs, or KEs (Figure 1—figure supplement 3A and C). Some marker genes of skin tissue were shown on the heatmap and dot plots, which examined the shared or species-specific genes in all cell types among the three species (Figure 1—figure supplement 3D and E). When discounting the unique skin thickness of CH pig resulting in this ratio discrepancy of the cell types such as excessive cell number of SMCs and ECs, dominantly originated from the vessel bed, we suggested that the pig skin tissue could be considered as the human skin model at single-cell levels for research purposes. Heterogeneity of skin SMCs in different anatomical sites SMCs play critical roles in forming blood vessels and APM in skin tissues (Driskell et al., 2013; Liu and Gomez, 2019); however, no study has uncharacterized skin SMCs at single-cell resolution. Here, we interrogated the heterogeneity and function of cutaneous SMCs. The t-SNE analysis divided SMCs into five subpopulations (clusters 0, 2, 5, 6, and 13) (Figure 2A), in which the MYH11 and ACTA2 marker genes were used for the immunohistochemistry staining of skin sections to validate the SMCs’ microanatomical sight (Figure 2B). Meanwhile, GO functional analysis was performed on the highly expressed genes for each cluster (Figure 2C). Clusters 0 and 13 predominantly took part in structural constituent of muscle, acting filament binding, and acting binding. The engagement of main inflammatory response and chemokine activity belonged to clusters 2, 5, and 6, of which cluster 2 was also involved in collagen binding and metallopeptidase activity. These results implied that SMCs play important roles in blood vessel homeostasis and function, partial collagen binding, and immune responses in skin tissue. Figure 2 Download asset Open asset Heterogeneity of skin smooth muscle cells (SMCs) in different anatomical sites. (A) The t-SNE plot visualization of SMCs including clusters 0, 2, 5, 6, and 13. (B) Confocal images showing immunofluorescence staining of ACTA2 (green) and MYH11 (red) in back skin sections, representative markers of SMCs. Scale bar = 50 μm. n = 3. (C) The enriched Gene Ontology (GO) terms of biological process for each SMC subpopulation were sorted by q-value. (D) Pseudotime ordering of SMC subpopulations using Monocle 2. (E) Heatmap illustrating the dynamics of representative differentially expressed genes among SMCs phenotypes, in which the important GO terms relating to biological processes were described. (F) These genes’ expression along pseudotime in SMC subpopulations. Then, we differentiated skin SMCs into other-like cell types. This evidence, combined with GO function analysis and the expression levels of conventional marker genes, such as MYH11 and ACTA2 for SMCs, GUCY1A2, CCL19, FGF7, and ASPN for mesenchymal cells (MECs), and LPL, CCL2, IL6, and CXCL2 for macrophages (MACs), presumed that cluster 2 might be mesenchymal-like phenotype while clusters 5 and 6 might represent macrophage phenotype. To further validate the topography of SMCs phenotypes, the analysis of pseudotime trajectory was performed by Monocle algorithm (Figure 2D). The trajectory demonstrated that SMCs experienced a dynamic transition from SMCs to mesenchymal-like phenotype and mesenchymal-like phenotype to macrophage-like phenotype. The sequential dynamics of gene expression with all branches were visualized and showed five gene sets along with expression patterns, which primarily deciphered three cell states (Figure 2E). These results discovered that gene sets 1 and 3 showed high expression levels of CTGF, LGR4, FABP4, CCL2, CCL19, and FGF7, and enriched GO terms associated with negative regulation of cell death, intestinal stem cell homeostasis, long-chain fatty acid transport, and immune response, which conformed well to the mesenchymal-like cells. Meanwhile, gene sets 2 and 5 showed high expression levels of MYH11, MYOM1, TPM1, TPM2, SQLE, BTG2, ADIRF, and TGFB3, and enriched GO terms involved in muscle contraction, actin filament organization, and ‘de novo’ action filament nucleation, which was greatly similar to contractile SMCs. In addition, gene set 4 showed high expression levels of CXCL10, CXCL2, ICAM1, LPL, and IL6, which were mainly gathered in GO terms of cellular response to lipopolysaccharide, cell chemotaxis, and defense response that may represent macrophage-like cells. Additionally, the expression levels of some cell type-specific marker genes in five SMCs clusters were presented (Figure 2F). The results showed the high expression levels of MECs-specific genes (GUCY1A2, FGF7, and CCL19) in cluster 2, MACs-specific LPL gene in clusters 5 and 6, and SMCs-specific MYH11 gene in clusters 0 and 13. These results proved our hypothesis that cluster 2 was mesenchymal-like phenotype while clusters 5 and 6 were macrophage-like phenotype. The cell number of skin SMCs showed a significant difference in different anatomical sites, while the distribution ratio of SMC subpopulations displayed a similar trend (Figure 3A). To decode the transcriptomic changes of skin SMCs in different anatomical sites, the DEGs were presented among 15 compared groups by pairwise comparison method (Figure 3B). GO analysis showed that the significant enriched terms for upregulated genes between compared groups primarily referred to extracellular region, collagen-containing ECM, and long-chain fatty acid transport, while these downregulated genes between compared groups took part in cytokine activity, CXCR chemokine receptor binding, and positive regulation of T cell migration (Figure 3C). The majority of upregulated genes subsisted in back skin compared to other locations, so we implemented KEGG analysis, which was involved in PI3K-Akt signaling pathway, MAPK signaling pathway, immune response, and integration (Figure 3D). Then, we chose some genes of related ECM-integrins and immune response to present their expression levels in different skin anatomical sites, which showed that immune response correlated closely with shoulder skin region or ECM-integrins tightly linked to skin locations on the head, back, and shoulder (Figure 3E). Moreover, to gain insights into gene expression regulation, we investigated the key transcription factors (TFs), which regulated the DEGs among various compared groups, using single-cell regulatory network inference and clustering (SCENIC). The SCENIC algorithm demonstrated a series of main regulons such as EGR1, ATF3, NFKB1, PRDM1, and REL, and their related target genes (Figure 3F). TFs, especially ATF3 and EGR1, primarily regulated their target genes at back skin. These results provide good insights into the inherent heterogeneity of skin SMCs in different anatomical sites. Figure 3 Download asset Open asset Heterogeneity of skin smooth muscle cells (SMCs) in different anatomical sites. (A) The cell number of SMC subpopulations in different skin regions. (B) Heatmap showing the differentially expressed genes of SMCs in multiple compared groups. Red represents upregulated genes, blue represents downregulated genes, and the number of differentially expressed genes is indicated. (C) The enriched Gene Ontology (GO) terms of multiple compared groups. Color indicates q-value, and circle indicates gene counts. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis for upregulated genes of back skin compared to other locations. (E) The expression level of genes involved in extracellular matrix (ECM)-integrins and immune response pathways in different skin regions. Red represents high expression of genes. (F) Transcriptional regulatory network of differentially expressed genes for SMCs in multiple compared groups. Blue nodes represent regulators and green nodes represent the target genes of regulators. Figure 3—source data 1 Source data of the differentially expressed genes of smooth muscle cells (SMCs) in multiple compared groups in Figure 3B. https://cdn.elifesciences.org/articles/86504/elife-86504-fig3-data1-v1.xlsx Download elife-86504-fig3-data1-v1.xlsx Heterogeneity of skin ECs in different anatomical sites Previous studies discovered that ECs underlie the vascular systems and primarily participate in blood and skin homeostasis (Kalucka et al., 2020). Here, ECs were captured from six different anatomical sites and were classified into five subpopulations, which were visualized with the t-SNE plot (Figure 4A). GO functional terms analysis was carried out according to the enriched expression genes for each cluster, which were closely related to some functional terms, including angiogenesis, immune response, response to viruses, cell migration, cell adhesion, and regulation of catalytic activity (Figure 4—figure supplement 1A). To validate the spatial position of ECs in dermis, we detected the expression levels of representative PECAM1 and APOA1 genes via the immunofluorescence of skin section (Figure 4B). Figure 4 with 2 supplements see all Download asset Open asset Heterogeneity of skin endothelial cells (ECs) in different anatomical sites. (A) The t-SNE plot visualization of ECs. (B) Immunofluorescence staining of APOA1 (green) and PECAM1 (red) in back skin sections, representative markers of ECs. Scale bar = 50 μm. n = 3. (C) Dot plot representing marker genes of ECs phenotypes. Color indicates gene expression, and circle indicates the log2FC value. (D) Pseudotime ordering of ECs subpopulations using monocle 2. (E) Heatmap showing the gene expression of metabolic pathways in ECs subpopulations. (F) Heatmap of gene expression of metabolic pathways in ECs subpopulations of different skin regions. (G) Heatmap of differentially expressed genes (DEGs) for ECs in multiple compared groups. Red represents upregulated genes, and blue represents downregulated genes. (H) The significantly enriched Gene Ontology (GO) terms of ECs in multiple compared groups. (I) Regulatory network of DEGs for ECs of different skin regions. Blue nodes represent regulators, and green nodes represent the target genes of regulators. Figure 4—source data 1 Source data of the differentially expressed genes of endothelial cells (ECs) in multiple compared groups in Figure 4G. https://cdn.elifesciences.org/articles/86504/elife-86504-fig4-data1-v1.xlsx Download elife-86504-fig4-data1-v1.xlsx Further, we classified EC phenotypes according to previous reported methods (Li et al., 2021; Wang et al., 2022) and found that ECs were composed of arteriole ECs expressing markers SEMA3G and MECOM (clusters 7 and 10), capillary ECs expressing marker PLVAP (cluster 3), venule ECs expressing markers SELE and ACK1 (cluster 4), and lymphatic ECs expressing markers LYVE1 and PROX1 (cluster 11) in dermis (Figure 4C). The pseudotime trajectory analysis of EC phenotypes showed an organized axis of blood ECs starting from arteriole and ending at venule, and it also formed an arteriovenous anastomosis tendency (Figure 4D). Due to exhibiting diverse molecules and functions on EC phenotypes, we further explored the expression levels of these EC phenotype-related genes, which were involved in integrins (focal adhesion), immune (cell adhesion molecules, chemokine signaling pathway, antigen processing and presentation, leukocyte transendothelial migration, and Th1 and Th2 cell differentiation), and metabolism (inositol phosphate, mucin type o-glycan biosynthesis, ether lipid, sphingolipid, and glycerolipid). ECs-related metabolism in our dataset was considerably active in arteriole ECs, especially cluster 10 involving ACER3, which controlled the homeostasis of ceramides, and LCLAT1, a lysocardiolipin acyltransferase-regulating activation of mitophagy (Figure 4E). The focal adhesion genes were more significantly upregulated in arteriole ECs and lymphatic ECs compared to other phenotypes, including ACTG1, BIRC3, and THBS1 (Figure 4—figure supplement 1B). In cell adhesion molecules, PTPRM and CDH5, mainly responsible for intercellular adhesion between ECs, were highly enriched in arteriole ECs; meanwhile, PECAM1, SELE, and SELP were enriched in venule ECs (Figure 4—figure supplement 1C). Other immune pathways showed that different EC phenotypes significantly highly expressed diverse genes, such as CXCL14 (involved in monocyte and recruitment) in capillary ECs, CCL26, CXCL19, and CCL26 in venule ECs (Figure 4—figure supplement 1D). The observed functional diversity of EC phenotypes proved the degree of ECs heterogeneity. Here, we found that the cell number of EC phenotypes was different among different anatomical sites, with the back skin holding the most arteriole ECs and minimal lymphatic ECs (Figure 4—figure supplement 1E and F). To further confirm the heterogeneity of EC phenotypes in the six different anatomical sites, we compared the expression levels of these gene-related integrins, immune, and metabolism pathways (Figure 4F and Figure 4—figure supplement 2A–C). For example, for the metabolism pathway, compared to other sites, the activity of capillary ECs, venule ECs, and arteriole ECs (cluster 7 not including cluster 10) was depressed in shoulder skin, while high activity in capillary ECs, venule ECs, and arteriole ECs was shown in leg skin, including ACER3 and PIK3C2A, which enhanced cell viability (Gulluni et al., 2021), and high activity in lymphatic ECs, arteriole ECs (cluster 10 not including cluster 7), and venule ECs was presented in ear skin, including ACER3, GALNT10, and PIK3C2B, a member of class II PI3Ks controlling cellular proliferation, survival, and migration. The obtained abundant results on the gene expressions for EC phenotype-related pathways in different sites showed the heterogeneity of skin ECs for different anatomical sites. To uncover the underlying molecular mechanisms driving the differential skin sites of ECs, we compared the DEGs with differentially compared groups among different anatomical sites (Figure 4G) and GO terms were implemented for upregulated and downregulated genes (Figure 4H). Enriched terms, which were related to long-chain fatty acid transport, lipoprotein particle binding, and ECM, were shown in upregulated differential genes, while downregulated differential genes mainly existed in terms involving regulation of catalytic activity, acting binding, and molecular adaptor activity. Of note, the CD36, a multifunctional fatty acid transporter, was related to the metabolic state of fibroblasts for ECM regulation (Zhao et al., 2019). Here, we found that CD36 was upregulated in the back compared with others, except head, which was enriched in metabolic terms such as long-chain fatty acid transport and regulation of nitric oxide. FABP4, fatty acid-binding protein 4, was significantly differentially expressed in nine pairs compared groups. A pioneering study showed that FABP4 was strongly expressed in subcutaneous adipocytes and adipose ECs (Wang et al., 2022). Combining data indicated that skin thickness might have a positive correlation with subcutaneous fat deposits. Additionally, we constructed single-cell transcription-factor regulatory networks with all DEGs of ECs (Figure 4I). The analysis predicted the following main transcriptional factors: ATF3, EGR1, ERG, FLI1, PRDM1, and NFKB1. The expression levels of ATF3, EGR1, and ERG were predominantly regulated in back and leg sites. With these findings, we presented the heterogeneity of skin ECs in different anatomical sites. Heterogeneity of skin FBs in different anatomical sites The dermal FBs synthesize the ECM that forms the connective tissue of skin dermis to maintain the skin morphology such as thickness and homeostasis (Zhao et al., 2019). Phenotypic data showed that the CH skin thickness of differential anatomical sites showed striking difference such as back skin thickness on average at 5.48 mm and that of the ear at 1.52 mm (Figure 5—figure supplement 1A). In terms of overall skin section, the skin histomorphology of different anatomical sites exhibited some difference in sparsity of collagen fibers or the number of appendages, and dermal thickness descended from the back, head, shoulder, leg, abdomen, to the ear (Figure 5A and Figure 5—figure supplement 1B). Curiously, we inquired whether the discrepancy in ECM accumulation in different skin sites was caused by FBs heterogeneity. Next, the FBs single-cell maps were presented from six different skin anatomical sites using the t-SNE plot, which was established by four clusters (clusters 1, 8, 9, and 12) (Figure 5B), and the cell number of clusters was estimated (Figure 5C). Similar to previous reports (Philippeos et al., 2018; Solé-Boldo et al., 2020), here FBs in cluster 1 highly expressed MGP and MFAP5, known markers of reticular FBs, the most representative markers of COL6A5, WIF1, and APCDD1 of papillary FBs belonged to clusters 8 and 9, and the mesenchymal subpopulation signature was typically characterized by enriched expressed CRABP1, TNN, and SFRP1 in cluster 12. GO analysis showed that the functions of these four FBs clusters were closely related to ECM organization, collagen fibril organization, and cell adhesion (Figure 5—figure supplement 1C). Likewise, the label-LUM and POSTN genes were marked on FBs of skin section via immunofluorescence (Figure 5D). Figure 5 with 2 supplements see all Download asset Open asset Heterogeneity of skin fibroblasts (FBs) in different anatomical sites. (A) Skin section with HE staining (left) and dermal thickness of six different sites (right) including head, ear, back, shoulder, abdomen, and leg. Scale bar = 100 μm. n = 3. (B) The t-SNE plot showing FBs populations. (C) The cell number of FBs populations in different skin regions. (D) Images showing immunofluorescence stain