10× genomics scRNA-seq data of regenerative [GSM4213633; large full-thickness excision (1 cm²) allows de novo follicle regeneration] and fibrotic (GSM4213632; large wounds lead to hairless scars) wound-induced hair neogenesis (WIHN) wounds of adult 6- or 7-week-old C57Bl/6j mice were curated from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds/). The accession number was GSE141814 (11). Regenerative wounds were defined as hair neogenesis, decreased contraction, decreased Wnt and TGF-β signaling activity, and decreased collagen production, while fibrotic wounds were defined as decreased hair neogenesis, increased contraction, increased Wnt and TGF-β signaling activity, and increased collagen production. This dataset was based on the platform of GPL21103 Illumina HiSeq 4000 (Mus musculus).

The DropletUtils package (v 3.13) was adopted to read unique molecular identifiers (UMI) count matrix, identify cells from empty droplets, remove barcode-swapped pseudo-cells, and downsample the count matrix (12). The calculateQCMetrics function of the Scater package was used for counting the expression of genes in cells (13). Cells with proportions of mitochondrial genes ≤ 10% and ribosomal genes ≥ 10% were determined for further analysis.

The expression matrix was normalized with the NormalizeData function of the Seurat package (14). The top 2,000 highly variable genes were screened by the FindVariableFeatures function. Then, expression data were linearly scaled utilizing the ScaleData function. Finally, principal component analysis (PCA) was performed with the RunPCA function based on the 2,000 genes.

The principal components with large standard deviations were selected. Then, cell clustering analysis was performed using the FindNeighbors and FindClusters function of the Seurat package. With the RunUMAP function, Uniform Manifold Approximation and Projection (UMAP) was carried out for dimension reduction. Cell types were annotated on the basis of the known marker genes.

To calculate the differentially expressed genes between each cluster and all other cells, the FindAllMarkers function of the Seurat package was used and novel marker genes were identified according to the following criteria: |log fold change (FC)| ≥ 0.1, the minimum expression ratio of cell population = 0.25, and p-value ≤ 0.05.

Based on the ligand–receptor pairs from the previous literature (15), the relationship pairs of receptors and ligands were analyzed based on the marker genes of various cells. Then, a cell–cell communication network was conducted and visualized with the Cytoscape software (16). The core cell population was identified according to the largest number of receptor–ligand pairs in the network. Moreover, the receptor and ligand genes were extracted.

Function enrichment analysis of the indicated genes was carried out utilizing the clusterProfiler package, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (17). GO categories contain biological process, cellular component, and molecular function. Terms with p < 0.05 were considered significantly enriched.

The Search Tool for the Retrieval of Interacting Genes (STRING) database (version 11.0; https://string-db.org/) was utilized for exploring the functional interactions between marker gene-encoded proteins (18). Then, PPI networks were constructed and the top 20 hub genes were identified.

Pseudotime analysis was carried out with the Monocle 3 tool (19). Firstly, genes that were expressed in at least 5% of the cells were selected. Then, the reduceDimension function was utilized to perform dimensionality reduction analysis, followed by cell cluster with the clusterCells function. Afterwards, the differentialGeneTest function was adopted to determine candidate genes with differences between the clusters with p < 0.05. The dimensionality reduction analysis of the cells was carried out using the DDRTree approach and the reduceDimension function based on the candidate genes. Through the orderCells function, the cells along the quasi-chronological trajectory were sorted and visualized.

The single-sample gene set enrichment analysis (ssGSEA) function of the Gene Set Variation Analysis (GSVA) package was utilized for comparisons of the differences in GO and KEGG terms between groups (20).

C57BL/6 male mice (8–10 weeks old; Sankyo) were used for fibroblast isolation. Briefly, mice were sacrificed by cervical dislocation. The trunk skin was separated in the ultra-clean bench, immersed in 75% ethanol for disinfection, and then cut into small pieces. Blood was removed by rinsing with PBS buffer and transferred evenly to cell culture dishes. DMEM complete medium (Wako) was added to submerge the tissue block that was placed in a constant temperature incubator to fully cultivate. After 24 h, DMEM complete medium was added, which was replaced every 3 days. The mouse skin fibroblasts were purified by the differential adhesion method and were used for subsequent experiments. Our study was approved by the Animal Ethics Committee of Keio University School of Medicine [12090(5)].

Using the TransIT-TKO Transfection Reagent (Mirus), siRNA-Engrailed-1 (horizon) and siRNA-control were transfected into fibroblasts in a constant-temperature incubator. Forty-eight hours later, the knockdown effect of siRNA was confirmed by real-time quantitative polymerase-chain reaction (RT-qPCR).

Total RNA was extracted from fibroblasts using the Isogen reagent (Nippon Gene) following the manufacturer’s instructions. cDNA synthesis was achieved based on the cDNA Synthesis System (Bio-Rad). RT-qPCR was carried out utilizing SYBR Qpcr Mix (Toyobo) on a 7500 Real-Time PCR system (Applied Biosystems). The primer sequences were as follows: EN1, 5’-ACACAACCCTGCGATCCTACT-3’(forward) and 5’-GGACGGTCCGAATAGCGTG-3’ (reverse); ACTB, 5’-GGC TGTATTCCCCTCCATCG-3’(forward) and 5’-CCAGTTGGTAACAATGCCATGT-3’ (reverse). The relative expressions were calculated with the 2^−ΔΔCt method.

Fibroblasts were plated onto a 6-well plate (about 3 × 10⁵ cells/well). When the confluence reached 100%, the fibroblast monolayer was scratched with a 1000-μl pipette tip. Additionally, detached fibroblasts were removed with serum-free medium. At 0 h and 24 h, the wounded area was photographed.

All statistical analysis was performed using the R language (version 3.6.1) and R Bioconductor packages. p < 0.05 indicated statistical significance.

Herein, we collected scRNA data of dermal cells from large skin wounds on day 18 with two distinct healing fates (fibrosis: GSM4213632 or regeneration: GSM4213633) from the GSE141814 dataset. Before analysis, we presented quality control of scRNA data. Barcode rank plots separately depicted the distribution of barcodes in total UMI count for fibrotic and regenerative wound dermal cells (Supplementary Figures 1A, B). Knee and inflection points in the barcode rank plots indicated the transition of the total UMI count distribution, which reflected the difference between empty droplets and cell droplets. After filtrating empty droplets, we counted the expression of genes in each cell (Supplementary Figures 1C, D). Afterwards, we filtrated out cells with proportions of mitochondrial genes > 10% and ribosomal genes < 10% (Supplementary Figures 1E, F).

After normalizing scRNA data, we screened the top 2,000 highly variable genes across fibrotic and regenerative wound dermal cells (Figure 1A). Then, scRNA data were linearly scaled and analyzed by dimensionality reduction with PCA. Here, we screened the top two principal components for subsequent analysis (Figure 1B). PCA results uncovered the prominent difference between fibrotic and regenerative wound dermal cells (Figure 1C). According to the elbow point, we identified the optimal principal components as 8 (Figure 1D). Heatmaps depicted the top 20 marker genes in each principal component (Figure 1E). With the UMAP method, dermal cells were clustered into 15 clusters (Figure 1F). The top ten marker genes of each cell cluster are presented in Figure 1G.

This study attempted to identify cell types across fibrotic and regenerative wound dermal cells. Based on the known marker genes, six cell types were annotated, as follows: EN1-negative myofibroblasts (n = 6,392), EN1-positive myofibroblasts (n = 2,219), hematopoietic cells (n = 3,774), macrophages (n = 1,461), pericytes (n = 1,493), and endothelial cells (n = 303; Figure 2A). Table 1 lists the cell ratio of each cell type. In particular, we noticed the differences in ratios of EN1-negative and -positive myofibroblasts between fibrotic and regenerative wound dermal cells (Figure 2B). With |logFC| ≥ 0.1, the minimum expression ratio of cell population = 0.25, and p-value ≤ 0.05, we identified novel marker genes in each cell type (Supplementary Table 1). The top ten marker genes in each cell type were visualized, as follows: EN1-negative myofibroblasts (Aebp1, Col1a1, Col1a2, Col3a1, Col8a1, Dcn, Eln, Mfap2, Mfap4, and Sparc), hematopoietic cells (AW112010, Cd3d, Cd3g, Cd52, Hcst, Ltb, Ptprcap, Rac2, Srgn, and Trbc2), macrophages (Apoe, C1qb, Ccl9, Cd74, Ctss, Fcer1g, H2-Eb1, Lyz2, Ms4a6c, and Tyrobp), pericytes (Acta2, Col4a1, Col4a2, Gm13889, Higd1b, Myl9, Mylk, Rgs5, Sparcl1, and Tagln), EN1-positive myofibroblasts (Birc5, Pclaf, Stnm1, Ube2c, Hist1h2ap, Col5a3, Cks2, Aqp1, Tnfaip6, and Timp1), and endothelia cells (Egfl7, Cldn5, Cdh5, Ramp2, Ecscr, Pecam1, Cd200, Ltbp4, Aqp1, and Hist1h2ap) (Figure 2C). Furthermore, we detected the expression levels of the known marker genes that were used for annotating cell types, as follows: endothelial cells (Cldn5, Pecam1, and Cd74), EN1-negative and -positive myofibroblasts (En1, Col1a1, Dcn, Sfrp4, Fndc1, and Lum), macrophages (Cd14, Cd68, and Csf1r), and hematopoietic cells (Ptprc, Cd69, Acta2, and Rgs5) (Figures 2D–J).

Wound healing is a complex process that necessitates the collaborative efforts of diverse cell lineages (21). Cell-to-cell communications across diverse cell types thoroughly govern appropriate functions of metazoans as well as widely rely on interactions between secreted ligands and cell-surface receptors. Based on the marker genes, ligand–receptor interactions were matched. The number of ligands/receptors for myofibroblasts, pericytes, endothelial cells, macrophages, and hematopoietic cells was 114, 91, 32, 28 and 17, respectively (Figure 3A). According to the number of intercellular receptor–ligand pairs, we screened out myofibroblasts as the core cell population.

We further evaluated the biological functions of ligand and receptor genes between myofibroblasts and macrophages. Our results demonstrated that ligand and receptor genes between myofibroblasts and macrophages were mainly involved in tube morphogenesis and development, regulation of cell migration, and motility (Figure 3B). Moreover, we found that the TGF-β signaling pathway was markedly enriched by these ligand and receptor genes between myofibroblasts and macrophages (Figure 3C).

We further verified the effects of EN1 on the migration of fibroblasts. Firstly, siRNA against EN1 was designed and transected into fibroblasts. RT-qPCR demonstrated that EN1 mRNA expression was distinctly reduced following siRNA-EN1 transfection (Figure 3D). According to wound healing results, EN1-knockout fibroblasts displayed significantly enhanced migration capacity (Figures 3E, F). Hence, EN1 suppression enabled to facilitate fibroblast migration.

With the cutoffs of |FC| > 1.2 and p < 0.05, we identified 546 up- and 481 downregulated specific genes in regenerative compared to fibrotic myofibroblasts (Figures 4A–C). Table 2 lists the first 20 up- and downregulated specific genes between regenerative and fibrotic myofibroblasts. As depicted in Figure 4D, we observed that the specific genes markedly participated in collagen-containing extracellular matrix, posttranscriptional regulation of gene expression, positive regulation of cell migration, mRNA metabolic process, and apoptotic signaling pathway. Moreover, ribosome and thermogenesis were prominently enriched by the specific genes (Figure 4E).

With the cutoffs of |FC| > 1.2 and p < 0.05, we found that 100 specific genes were significantly upregulated while 197 specific genes were significantly downregulated in regenerative compared to fibrotic macrophages (Figures 5A–C). Table 3 lists the first 20 up- and downregulated specific genes between fibrotic and regenerative macrophages. GO enrichment analysis uncovered that the specific genes were markedly involved in the negative regulation of programmed cell death, the regulation of cell migration, innate immune response and apoptotic signaling pathway, collagen-containing extracellular matrix, the positive regulation of T cell activation, and response to interferon γ (Figure 5D). Moreover, we observed that antigen processing and presentation, pathways in cancer, phagosome, ribosome, and tuberculosis were prominently enriched by the specific genes (Figure 5E).

With the STRING tool, we probed the interactions between myofibroblast- or macrophage-specific gene-encoded proteins. In Figure 6A, there were 616 nodes in the PPI network of myofibroblasts, reflecting the close interactions of myofibroblast-specific gene-encoded proteins. According to degree, the top 20 nodes were identified as hub genes, including Rps27a, Rps11, Rps23, Rps3, Rps5, Rps15a, Rps6, Rps9, Rps13, Rps14, Rps25, Rps3a1, Rps27, Rps8, Rps19, Rps28, Rps7, Rpl8, Rps18, Rpl26, Rpl32, and Rps16, indicating that the above genes were the core of the network. Figure 6B depicts the interactions between macrophage-specific gene-encoded proteins. The 20 hub genes were as follows: Uba52, Rps9, Gnb2l1, Rpl27, Rpl38, Rps13, Rps15a, Fau, Rpl18, Rpl30, Rpl35a, Rpl7, Rplp2, Rps24, Rpl13a, Rpl4, Rps10, Rps12, Rps27rt, and Rps2. The above genes deserve in-depth explorations.

To investigate the underlying evolution among myofibroblasts and macrophages, this study adopted the Monocle tool to reveal a pseudotemporal ordering for the similarity of cell clusters with developmental lineages. For myofibroblasts, the results clearly demonstrated the uniform development of myofibroblasts from cluster 6 to cluster 10 (Figure 7A). The trends of pseudotime‐dependent genes along the pseudo‐timeline were divided into six cell clusters of myofibroblasts with diverse expression dynamics. Furthermore, we observed that macrophage under fibrotic conditions was in the beginning position of the differentiation process and was sequentially transformed into macrophage under regenerative conditions (Figure 7B).

According to the results of pseudotime analysis of myofibroblasts, we carried out GSVA between the initially differentiated cluster 6 and the final differentiated cluster 10. Compared with cluster 10 of myofibroblasts in fibrotic and regenerative dermal cells, biological processes such as the metabolic process significantly activated cluster 6 of myofibroblasts in fibrotic and regenerative dermal cells (Figure 8A). As depicted in Figure 8B, we noticed the prominent activation of cellular components such as mitochondria in cluster 6 of fibrotic and regenerative myofibroblasts in comparison to those in cluster 10. Moreover, we observed that fibrotic and regenerative myofibroblasts in cluster 6 had significantly activated molecular functions like oxidoreductase activity compared with fibrotic and regenerative myofibroblasts in cluster 10 (Figure 8C). We also compared the differences in KEGG pathways between clusters. Diverse signaling pathways like metabolic pathways, RNA transport, spliceosome, thermogenesis, oxidative phosphorylation, carbon metabolism, ribosome, cell cycle, protein processing in the endoplasmic reticulum, and biosynthesis of amino acids were prominently activated in fibrotic and regenerative myofibroblasts in cluster 6 compared to those in cluster 10 (Figure 8D).

GSVA was also presented between fibrotic and regenerative macrophages. In Figure 9A, we determined that biological processes such as the metabolic process and immune response were markedly activated in fibrotic macrophages compared to regenerative macrophages. The significantly activated cellular components such as the spliceosomal complex, catalytic complex, ribonucleoprotein complex, nuclear lumen, nucleoplasm, nucleolus, cytosol, nucleus, catalytic step 2 spliceosome, chromosome, and protein-containing complex were found in fibrotic macrophages compared with regenerative macrophages (Figure 9B). As shown in Figure 9C, we investigated the marked activation of molecular functions like RNA binding, ATP binding, mRNA binding, adenyl ribonucleotide binding, adenyl nucleotide binding, drug binding, nucleic acid binding, heterocyclic compound binding, organic cyclic compound binding, and ATPase activity in fibrotic macrophages in comparison to regenerative macrophages. Moreover, our results showed that KEGG pathways such as spliceosome, NOD-like receptor signaling pathway, Fc gamma R-mediated phagocytosis, antigen processing and presentation, endocytosis, necroptosis, and natural killer cell-mediated cytotoxicity displayed marked activation in fibrotic macrophages compared to regenerative macrophages (Figure 9D).