Single-Cell RNA-seq Analysis Reveals Cellular Functional Heterogeneity in Dermis Between Fibrotic and Regenerative Wound Healing Fates

Background Fibrotic scars are common in both human and mouse skin wounds. However, wound-induced hair neogenesis in the murine wounding models often results in regenerative repair response. Herein, we aimed to uncover cellular functional heterogeneity in dermis between fibrotic and regenerative wound healing fates. Methods The expression matrix of single-cell RNA sequencing (scRNA-seq) data of fibrotic and regenerative wound dermal cells was filtered, normalized, and scaled; underwent principal components analysis; and further analyzed by Uniform Manifold Approximation and Projection (UMAP) for dimension reduction with the Seurat package. Cell types were annotated, and cell–cell communications were analyzed. The core cell population myofibroblast was identified and the biological functions of ligand and receptor genes between myofibroblast and macrophage were evaluated. Specific genes between fibrotic and regenerative myofibroblast and macrophage were identified. Temporal dynamics of myofibroblast and macrophage were reconstructed with the Monocle tool. Results Across dermal cells, there were six cell types, namely, EN1-negative myofibroblasts, EN1-positive myofibroblasts, hematopoietic cells, macrophages, pericytes, and endothelial cells. Ligand and receptor genes between myofibroblasts and macrophages mainly modulated cell proliferation and migration, tube development, and the TGF-β pathway. Specific genes that were differentially expressed in fibrotic compared to regenerative myofibroblasts or macrophages were separately identified. Specific genes between fibrotic and regenerative myofibroblasts were involved in the mRNA metabolic process and organelle organization. Specific genes between fibrotic and regenerative macrophages participated in regulating immunity and phagocytosis. We then observed the underlying evolution of myofibroblasts or macrophages. Conclusion Collectively, our findings reveal that myofibroblasts and macrophages may alter the skin wound healing fate through modulating critical signaling pathways.


INTRODUCTION
The skin is the organ with the largest surface area in the human body that provides an efficient protective barrier against mechanical injury, microbial pathogens, and trauma (1). The skin's immune system is divided into two structural compartments: epidermis and dermis, both of which contain a plethora of immunocompetent cell types (2). The epidermis is home to the main skin-resident immune cells, Langerhans cells, and melanocytes. Meanwhile, immune-specialized cells like dendritic cells, macrophages, and T cells reside in the dermis (3). The communications within immune populations and the skin environment are critical to the effectiveness of the skin immune system (4). Wound healing is a complex process in the human body, where numerous cell populations with different functions are involved in the stages of hemostasis, inflammatory response, growth, re-epithelialization, and remodeling (5). It is essential to repair the skin after damage (6). Skin wound healing involves three primary phases: inflammation, re-epithelialization, and tissue remodeling (7). Nevertheless, effective therapeutic strategies of accelerating healing and decreasing scarring remain lacking. Single-cell RNA sequencing (scRNA-seq) technology has emerged as an indispensable tool for elucidating cellular phenotype and functional heterogeneity (8). Deciphering the role of each cell type and interactions within cells is of importance to understand the mechanism of normal wound closure (9). Alterations in the microenvironment may influence cellular recruitment or activation, resulting in damaged states of wound healing. ScRNA-seq can be applied for deciphering the cellular changes in chronic wounds and hypertrophic scarring, thereby promoting the development of more effective therapeutic solutions for healing wounds (10). Moreover, in-depth understanding of the differences between fibrotic and regenerative wound healing fates is a prerequisite for developing more effective therapeutic interventions (2). Here, the purpose of this study was to reveal cellular functional heterogeneity in the dermis between fibrotic and regenerative wound healing fates.

Quality Control
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.

Data Preprocessing and Principal Component 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.

Cell Cluster and Annotation
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.

Identification of Novel 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 pvalue ≤ 0.05.

Ligand-Receptor Network Analysis
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
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. exploring the functional interactions between marker geneencoded proteins (18). Then, PPI networks were constructed and the top 20 hub genes were identified.

Pseudotime Analysis
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.

Gene Set Variation Analysis
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).

Isolation and Culture of Fibroblasts
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)].

Transfection
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).

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

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

Quality Control of scRNA-seq Data of Fibrotic and Regenerative Wound Dermal Cells
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).

Cell Cluster of Fibrotic and Regenerative Wound Dermal Cells
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.

Cell-Cell Interactions Based on Ligand-Receptor Interactions
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.

Biological Functions of Ligand and Receptor Genes Between Myofibroblasts and Macrophages
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-b signaling pathway was markedly enriched by these ligand and receptor genes between myofibroblasts and macrophages ( Figure 3C).

Knockdown of EN1 Facilitates Fibroblast Migration
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.

Identification of Specific Genes Between Fibrotic and Regenerative Myofibroblasts and Their Biological Functions
With the cutoffs of |FC| > 1.2 and p < 0.05, we identified 546 upand 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).

Identification of Specific Genes Between Fibrotic and Regenerative Macrophages and Their Biological Functions
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 g ( 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).

PPI Network Analysis of Specific Genes Between Fibrotic and Regenerative Myofibroblasts or Macrophages
With the STRING tool, we probed the interactions between myofibroblast-or macrophage-specific gene-encoded proteins.

Reconstruction of the Temporal Dynamics of Myofibroblast and Macrophage
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 pseudotimedependent 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).

GSVA Between Clusters 6 and 10 of Fibrotic and Regenerative Myofibroblasts
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 Between Fibrotic and Regenerative Macrophages
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).

DISCUSSION
Skin wound healing involves complicated coordinated interactions within cells. Through scRNA-seq data, this study identified six cell   populations, namely, EN1-negative myofibroblasts, EN1-positive myofibroblasts, hematopoietic cells, macrophages, pericytes, and endothelial cells, across the dermis. Evidence suggests that EN1positive fibroblasts are known to function in scarring, and EN1negative fibroblasts yield wound regeneration. Thus, we used EN1 as a marker to divide the subgroups. Dynamic cellular events after skin injury rely on bidirectional cell-cell communications against effective wound healing (22). Our results demonstrated the crosstalks between myofibroblasts, hematopoietic cells, macrophages, pericytes, and endothelial cells in the dermis based on the ligandreceptor interactions. As per previous studies, CX3CR1 may mediate the recruitment of bone marrow-derived monocytes or macrophages in skin wound healing, thereby releasing profibrotic as well as angiogenic mediators (23). Moreover, macrophages support proliferation and heterogeneity of myofibroblasts in skin repair (24). Serum endothelial cell-derived extracellular vesicles facilitate diabetic wound healing via enhancing myofibroblast proliferation and decreasing senescence (25). Intradermal adipocytes modulate the recruitment of myofibroblasts in skin wound healing (26). Fibroblasts promote NG 2+ pericyte populations in murine skin development as well as repair (27). On the basis of the above lines of evidence, there were remarkable interplays between diverse cell types during dermis progression. According to the number of ligands and receptors, we identified myofibroblasts as the core cell population. Our function enrichment analyses uncovered that the ligand and receptor genes between myofibroblasts and macrophages were mainly involved in regulating cell proliferation and migration, tube development, and the TGF-b pathway. The TGF-b signaling pathway plays an important role in the formation of collagen in fibroblasts and myofibroblasts (28). Cytokine TGF-b may induce dermal dendritic cells to express IL-31, thereby activating sensory neurons as well as stimulating wound itching during skin would healing (29). Hence, targeting the TGF-b pathway is the promising therapeutic intervention to reduce abnormal skin scar formation.
To explore the differences in molecular mechanisms involving myofibroblasts between fibrotic and regenerative wound healing fates, we identified 546 up-and 481 downregulated specific genes in regenerative compared to fibrotic myofibroblasts. This revealed the heterogeneity of myofibroblasts between fibrotic and regenerative wound healing. Our GO and KEGG enrichment analysis uncovered the key biological functions involving the specific genes between fibrotic and regenerative myofibroblasts. As a result, these specific genes between fibrotic and regenerative myofibroblasts prominently participated in the mRNA metabolic process and organelle organization. Extracellular matrix of connective tissues is synthesized by myofibroblasts that play a critical role in sustaining the structural integrity of various tissues (30).
Skin wound macrophage is an important regulator of skin repair, and its dysfunction may cause chronic and non-healing skin wounds (31). Further analysis identified that 100 specific genes were significantly upregulated while 197 specific genes were significantly downregulated in regenerative compared to fibrotic macrophages. Functional enrichment analysis uncovered that these specific genes between fibrotic and regenerative macrophages primarily participated in regulating inflammatory response, immunity, and phagocytosis. Immunity is the most important function of the skin, which can prevent harmful exposure from the external and internal environment (32). Furthermore, late wound macrophage phagocytosis of the Wnt inhibitor may induce chronic Wnt activity during fibrotic skin healing (11). Collectively, our findings revealed that the heterogeneity of myofibroblasts or macrophages might determine wound healing fate as regenerative or fibrotic.

CONCLUSION
Taken together, this study uncovered cellular functional heterogeneity in dermis between fibrotic and regenerative wound healing fates. Moreover, myofibroblasts and macrophages may change the skin wound healing fates by modulating critical signaling pathways. Therefore, our data provided an insight into the development of more effective therapeutic interventions for improving healing fates.

ETHICS STATEMENT
Ethical review and approval were not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements. The animal study was reviewed and approved by Keio University School of Medicine. Written informed consent was not obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.