Single-cell analysis reveals fibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds

1Skin forms the outermost layer of the body, and principally consists of a stratified epidermis residing on top of a collagen-rich dermis. While epidermis endows skin with its barrier function, dermis provides mechanical strength and houses numerous epidermal appendages, principally hair follicles and sweat glands. Hair follicles are complex epithelial–mesenchymal mini-organs that are rich in stem cells and regenerate cyclically. When fully grown, hair follicles span the entire dermis and part of the dermal white adipose tissue (dWAT), where they engage in signaling crosstalk. As a result of this crosstalk, hair follicles induce adipocyte progenitor proliferation and adipocyte hypertrophy. Reciprocally, dWAT modulates hair stem cell quiescence and activation,.

2Upon significant injury, such as full-thickness excisional wounding, skin undergoes repair. While small wounds, <1 cm2, typically repair by forming scar devoid of epidermal appendages and fat, large wounds, larger than 1 cm2, can regenerate de novo hair follicles and adipocytes in their center. Large wounds in mice heal primarily by contraction, while the uncontracted portion closes by re-epithelialization and forms new connective tissue, rich in fibroblasts. In our model, wounds close in two weeks, and then new hair follicles regenerate in the central region by week three,, followed by new adipocytes during the fourth week. The process of de novo hair follicle regeneration, termed wound-induced hair neogenesis (WIHN), involves reactivation of embryonic hair development programs. Similarly, the process of de novo fat regeneration involves reactivation of an embryonic adipose lineage formation program (Supplementary Figure 1). It remains unclear why regeneration is limited to the wound center. Beyond laboratory mice,,, WIHN is observed in rodents from the genus Acomys and in rabbits. WIHN also likely occurs in sheep, but not in laboratory rats. While typically wounds in humans heal without regeneration, de novo hair follicles in facial skin wounds have been reported.

3Over the last decade, many signaling pathways necessary for WIHN have been partially elucidated. Activation of canonical WNT signaling is necessary for WIHN,,, and both epidermal and dermal wound cells secrete and respond to WNT ligands at distinct phases. Production of WNT ligands by dermal wound cells is initiated by FGF9, secreted by γδ T cells. Activation of Hedgehog signaling downstream of WNT is necessary for dermal papilla (DP) fate specification by wound fibroblasts and for successful WIHN. Macrophages also promote WIHN by secreting TNFα, which, in turn, activates p-AKT/p-β-catenin signaling. Activation of TLR3 signaling by double-stranded RNA augments IL6 production and STAT3 activation, both of which positively impact WIHN efficiency. In contrast, prostaglandin PDG2 inhibits hair follicle neogenesis. Furthermore, WIHN is modulated by several transcriptional regulators, including homeobox factor MSX2, zinc finger protein CXXC5, and RNA-binding protein MSI2. De novo fat regeneration is driven by BMP ligands produced by neogenic hair follicles. Responding to BMP signals, wound myofibroblasts activate ZFP423, a transcriptional regulator that drives adipogenic lineage commitment. Hair follicles are critical for fat neogenesis and no adipocytes regenerate in hairless wounds.

4Less is known about the lineage origin of cells during WIHN. Fate mapping experiments by Ito, Yang showed that progeny of the pre-existing KRT15-positive hair follicle bulge stem cells do not give rise to de novo hair follicles. Instead, progeny of non-bulge LGR6-positive and LGR5-positive epithelial stem cells can contribute toward neogenic hair follicles. Little is known about the cells that regenerate DP, the principal mesenchymal component of hair follicles. Lineage studies suggest that myofibroblasts are the source for neogenic DPs,. Tracing experiments on CD133-positive DP cells of preexisting hair follicles indicate that they do not mobilize upon wounding. Other fate mapping studies suggest that multiple skin fibroblast lineages mobilize to repair small wounds–. Using lineage tracing with En1 (Engrailed 1)-Cre, Rinkevich, Walmsley identified En1-positive and En1-negative dermal fibroblasts. After wounding, En1-positive fibroblasts become major contributors toward wound repair, and genetically depleting them leads to scar reduction. The contribution of En1-positive fibroblasts toward wound repair was confirmed by Shook, Wasko. Driskell, Lichtenberger showed that dermal fibroblasts contribute to wound repair in two waves. Progeny of lower dermal fibroblasts contribute early, while progeny of upper fibroblasts migrate into the wound later. Contribution from distinct dermal fibroblasts toward regenerating DPs in the WIHN model requires further investigation, yet the lack of clear lineage master-regulators complicates functional validation of DP fate mapping.

5Unlike DPs, fat tissue with its well-established master-regulators, including ZFP423, CEBPs, and PPARγ, provides a tractable model system for studying de novo regeneration. Recently, we showed that de novo adipocytes regenerate from Sm22 (aka Tagln) and Sma (aka Acta2) positive myofibroblasts. Myofibroblast-specific ablation of Pparγ or BMP receptor 1a (Bmpr1a) largely prevented adipocyte regeneration in otherwise hair-bearing wounds. However, the degree of wound myofibroblast heterogeneity and their competency for adipogenic reprogramming remains unclear. The advent of single-cell RNA-sequencing (scRNA-seq) enables profiling of cellular heterogeneity in tissues with poorly characterized cell types. In this study, using a scRNA-seq approach, we identify and characterize multiple distinct fibroblast populations in regenerating mouse wounds. We show that major populations co-exist in wounds across the time course of regeneration. Furthermore, we identify bone marrow-derived adipocytes and a rare subset of wound fibroblasts with myeloid characteristics that undergo fat regeneration.

6We performed scRNA-seq on unsorted cells from wound dermis 12 days post-wounding (PW) (Fig. 1a). This time point coincides with completion of wound re-epithelialization and strong SMA expression. Approximately 21,819 sequenced cells met quality control metrics (Supplementary Figure 2) and were analyzed. Unsupervised clustering using the Seurat package identified 13 cell clusters (Fig. 1b, left). Using the differentially expressed gene signatures, we attributed clusters to their putative identities (Fig. 1b, right) and hierarchical similarities (Fig. 1c; Supplementary Figure 3a). Figure 1d provides a summary diagram of identified cell types. Figure 1e–g show selected differentially expressed genes in the form of a heatmap (Fig. 1e), bar charts (Fig. 1f), and feature plots (Fig. 1g). Several clusters contained immune cells. The most abundant of them, representing ~16% of all cells, was cluster C3. It was enriched for myeloid markers, including C1qb, Cd14, Cd68, Lyz2, Mafb, and Pf4 (Supplementary Figure 3b; Supplementary Data 1). Cluster C7 cells were classified as T lymphocytes (~4%) and they expressed Cd3d, Cd81, Icos, Maf, Nkg7, and Thy1 (Supplementary Figure 4). Cluster C8 cells were identified as B lymphocytes (~3%) and C12 as dendritic cells (~1%). Two other distinct cell clusters were C5 (~9%) and C13 (~1%). Cluster C5 cells were enriched for endothelial markers Cav1, Cd34, Cd93, Ly6e, Ly6c1, and Pecam1, while cluster C13 cells—for lymphatic endothelial markers Ccl21a, Lyve1, Pdpn, and Prox1. Cluster C10 cells (~1%) were classified as Schwann cells based on expression of Cadm4, Dbi, Gpm6b, Itih5, and Plp1. The remaining five clusters—C1, C2, C4, C6 and C9, representing ~65% of all analyzed wound cells, were collectively characterized as fibroblasts. They were highly enriched for collagens, including Col1a1, Col3a1, Col5a1, Col12a1, and extracellular matrix protein genes, including Dcn (Decorin) and Fbln2 (Fibulin 2) (Supplementary Figure 3b; Supplementary Data 1). Many of these cells expressed high levels of contractile myofibroblast factors: Cald1, Myl9, Tagln, and Acta2 (Supplementary Figure 3b).

7Next, we performed unsupervised clustering on all wound fibroblasts and observed further heterogeneity with 12 subclusters, sC1 through sC12 (Fig. 2a, b; Supplementary Figure 5), each containing unique marker gene profiles (Fig. 2c; Supplementary Figure 6; Supplementary Data 2). Considering that transcription factors (TFs) commonly regulate cellular characteristics, we examined their expression patterns. All fibroblasts shared the following 20 TFs: Cebpb, Egr1, Fosb, Fosl2, Hif1a, Klf2, Klf4, Klf6, Klf9, Nfat5, Nfatc1, Nfkb1, Nr4a1, Nr4a2, Pbx1, Prrx1, Runx1, Stat3, Tcf4, and Zeb2 (Supplementary Figure 7a). These can be defined as a common wound fibroblast TF signature. Among these factors were Runx1, Tcf4, and Zeb2 (Supplementary Figure 8), previously implicated in myofibroblast differentiation. The Ebf1high/Id3high/Zeb2high/En1low/Nfixlow/Prrx2low/Sox9off signature marked fibroblasts from sC3/9/11 subclusters (Supplementary Figure 7b, c), which are hierarchically distinct (Fig. 2b; Supplementary Figure 5). Among the remaining nine subclusters, fibroblasts in six subclusters, sC1 and sC4–sC8, had low expression of Id2 and Id3. Among these six En1high/Id2low/Id3low subclusters, sC4 was Sox11high, sC5–Twist2high, sC7–Twist1high/Twist2high/Foxp1low and sC8–Nfiahigh (Supplementary Figure 7b, c). Other subclusters also had their own, albeit complex TF expression signatures.

8Next, we examined signaling pathway markers, receptors and ligands. Hierarchically distinct subclusters sC3/9/11 had the following receptor signature: Mcamhigh/Pdgfrbhigh/Fgfr1low/Tgfbr2low/Tgfbr3low/Ncam1off/Pdgfralow, and ligand signature: Il6high/Pdgfahigh/Igf1low/Igfbp3low/Mdklow/Dkk3off. Nine remaining subclusters were primarily differentiated from subclusters sC3/9/11 based on high Pdgfra expression, while having low Il6, Pdgfa, Pdgfrb and high Igf1, Mdk, Tgfbr2, Tgfbr3 expression (Supplementary Figure 7c). Additionally, among these nine subclusters, fibroblasts in sC2 were Angptl1high, sC5 –Ccl8high/Cxcl5high/Grem1high/Spp1high, sC6–Il1bhigh, sC7–Ccl8high/Igfbp3low, sC10–Angptl1high/Il1bhigh/Pf4high and sC12–Angptl1high/Fsthigh (Supplementary Figure 7b, c; Supplementary Data 2). We also profiled cell cycle state. Intriguingly, sC4 and sC11 subclusters were prominently enriched for G2/M markers (Supplementary Figure 9). Considering hierarchical similarities (Fig. 2b), sC4 and sC11 likely represent actively cycling subsets of sC1 and sC3 populations, respectively.

9We also spatially mapped a portion of fibroblast subclusters to day 12 wound regions using immunostaining. We chose markers based on their high, yet specific expression on scRNA-seq and their membranous or intracellular localization (Fig. 2b). Fibroblasts in clusters sC1/2/4-8/10/12 were marked based on PDGFRA expression. Among these, cells in subclusters sC1/4/5/7/8 were differentiated from sC2/6/10/12 cells based on CRABP1 expression. Intriguingly, PDGFRA+/CRABP1+ cells were enriched under the epidermis, while PDGFRA+/CRABP1neg cells primarily localized to the lower dermis (Fig. 2d, e). Among CRABP1+ cells, rare sC7 cells (~3%) were also CYP26B1+ and showed scattered presence in the dermis (Fig. 2f). Among PDGFRA+/CRABP1neg cells, rare sC6 cells (~4%) were G0S2+ and were seldom found in the lower dermis (Fig. 2g), while MEST2+ sC2/10 cells, which were abundant on scRNA-seq, showed broad presence throughout the lower dermis (Fig. 2h). sC3/9/11 fibroblasts were distinguished based on RGS5 expression. RGS5+ cells were present throughout wound dermis (Fig. 2i, j), and majority of them localized outside of the blood vessels.

10Taken together, our analyses suggest that, upon completion of re-epithelialization on day 12, large wounds may contain two major fibroblast populations. One population, representing ~24% of wound fibroblasts, consists of three subclusters sC3/9/C11, which express low levels of TGFβ receptors Tgfbr2, Tgfbr3, and PDGF receptor Pdgfra and high levels of Pdgfrb. They also stain with RGS5 and localize throughout wound dermis. The second and more heterogeneous population, representing ~76% of wound fibroblasts, consists of nine subclusters and expresses intermediate to high levels of Tgfbr2, Tgfbr3, high levels of Pdgfra, but not Pdgfrb. Notably, PDGFRA signaling is an established driver of fibrosis in multiple tissues, including fat and skeletal muscle. Pdgfrahigh cells also express higher levels of receptors for the TGFβ pathway, another well-established driver of fibrosis. Pdgfralow cells are also Pdgfrbhigh, and Pdgfrb-expressing perivascular cells were identified as the precursors for new adipocytes in visceral fat in obesity. Additionally, we identified a high level of previously unappreciated heterogeneity within both fibroblast populations. For instance, Pdgfrahigh fibroblasts further subdivide into upper and lower dermal cells both anatomically and by gene expression.

11While t-SNE analysis helped to reveal heterogeneity among wound fibroblasts, we also wondered if they share common differentiation trajectories. Indeed, upon wounding, many fibroblasts differentiate into myofibroblasts and on scRNA-seq contractile markers were detected across multiple fibroblast subclusters. In large wounds, myofibroblasts also serve as the principal progenitors for fat regeneration. Ordering of cells in pseudotime arranged most of fibroblasts into a major trajectory, with two minor bifurcations (Fig. 3a, top). Fibroblasts from different subclusters distributed broadly across the pseudotime space, with Pdgfralow sC3/9/11 cells primarily occupying the right half of the major trajectory, with the remaining half consisting of Pdgfrahigh cells (Fig. 3a, bottom). Next, we performed RNA velocity analysis, which considers both spliced and unspliced mRNA counts to predict potential directionality and speed of cell state transitions. It represents such transitions as vectors, with long vectors marking rapid differentiation events. RNA velocity distinguished three sets of vectors, defined as paths. Path1 and Path2 vectors appear to represent differentiation of the Crabp1neg fraction of Pdgfrahigh cells and a fraction of Pdgfralow cells, respectively, toward a shared state surrounding minor bifurcations, that is primarily composed of Crabp1+ Pdgfrahigh fibroblasts. Path3 vectors appear to represent differentiation of mainly Pdgfralow fibroblasts toward a distinct state (Fig. 3b). Cells expressing collagen Col14a1 and mesoderm-specific transcript Mest preferentially distributed at the beginning of all paths (Fig. 3d; Supplementary Data 3), while the contractile genes Acta2 and Tagln increased in density and expression levels toward the ends of these paths (Fig. 3e). We posit that all three RNA velocity paths may represent putative fibroblast differentiation trajectories toward myofibroblast states. We also examined pseudotime dynamics using scEpath, which identifies pseudotime-dependent genes and arranges them into clusters. This analysis revealed five gene clusters (Fig. 3c), including multiple signaling factors and TFs (Fig. 3f, g; Supplementary Figure 10). Overlaying RNA velocity vectors onto these clusters suggests that cluster pC1 contains genes expressed at the beginning points of Path1 and Path2, cluster pC2 primarily contains Path1 differentiation program genes, cluster pC3—mainly Path2 differentiation program genes, whereas clusters pC4 and pC5 are largely composed of Path3 differentiation program genes. Taken together, pseudotime, RNA velocity, and scEpath analyses establish a basis for exploring signaling and transcriptional regulators of wound myofibroblast differentiation programs.

12We noted that, on scRNA-seq, many fibroblasts across all 12 subclusters expressed hematopoietic markers. Specifically, many fibroblasts expressed the myeloid-specific marker Lyz2 (Lysozyme 2, aka LysM) (Fig. 4a). We examined fibroblasts that co-expressed Lyz2, collagen Col12a1 and contractile markers Acta2 and Tagln (Fig. 4b, c). This identified Lyz2+/Acta2+/Tagln+/Col12a1+ quadruple-positive myofibroblasts, representing ~11% of all wound fibroblasts. They were present in all subclusters with sC9 and sC10 being the most enriched, containing ~32% and ~21% of such cells, respectively (Fig. 4c). We also used two additional markers, hematopoietic progenitor-associated Cd34 and myeloid-specific Cd14 to identify presence of Cd34+/Acta2+/Tagln+/Col12a1+ and Cd14+/Acta2+/Tagln+/Col12a1+ quadruple-positive myofibroblasts (Fig. 4d).

13To examine if Lyz2+ myofibroblasts might derive from wound-resident myeloid cells, we performed pseudotime analysis jointly on all Lyz2-expressing wound fibroblasts and myeloid cells (cells from cluster C3 on Fig. 1b). This analysis arranged all cells into a major trajectory, with several minor bifurcations (Fig. 4e). Fibroblasts (orange) and myeloid cells (black) distributed differentially across the pseudotime space; however, a prominent overlap region was observed in the middle of the trajectory (Fig. 4f). Within this region, both cell types displayed hybrid expression patterns, simultaneously enriched for fibroblast-specific and myeloid-specific genes (Fig. 4g; Supplementary Figure 11a). Acta2 and Tagln were enriched in Lyz2+ fibroblasts, but were also expressed in Lyz2+ myeloid cells, including in the overlap region (Fig. 4h). scEpath identified five pseudotime-dependent gene clusters (Fig. 4i; Supplementary Figure 11b–e; Supplementary Data 4), while RNA velocity distinguished four vector paths (Fig. 4j). Importantly, Path3 vectors, which encompass the overlap region, appear to represent putative conversion of Lyz2+ myeloid cells into fibroblasts (red on Fig. 4j), while Path4 vectors appear to mark differentiation of wound-resident myeloid cells (green on Fig. 4j). Overlaying RNA velocity vectors onto scEpath data indicates that genes in cluster rC5 mark myeloid differentiation program along Path4, while genes in clusters rC4 and rC3 mark early and late stages of putative myeloid-to-fibroblast conversion pathway, respectively.

14Next, we performed scRNA-seq on genetically marked contractile cells, isolated from day 12 wounds of Sm22-Cre;tdTomato mice (Fig. 5a). A microfluidic-based capture platform was used (Supplementary Figure 12a). Although it captured substantially fewer cells, 116 in total, compared to the droplet-based capture approach, it yielded a higher depth of sequencing, 5830 vs. 1101 median genes/cell (Supplementary Figure 12b vs. 2b), and eliminated the possibility of cell doublets. On t-SNE analysis, cells segregated into three clusters—fC1, fC2, and fC3 (Fig. 5b), each with its distinct gene markers (Supplementary Figure 13; Supplementary Data 5). Pearson correlation showed that fC1 cells were the most similar to Pdgfrahigh fibroblasts from subclusters sC1/4/5/8/10. Cells from fC2 were the most similar to Pdgfralow fibroblasts from subclusters sC3/11 and expressed high levels of Pdgfrb. Cluster fC3 had the highest correlation to sC9 subcluster (Fig. 5d), enriched for Lyz2+/Acta2+/Tagln+/Col12a1+ quadruple-positive myofibroblasts. Consistently, this cluster was also the most enriched for hematopoietic markers, including pan-hematopoietic Cd45 (aka Ptprc), and myeloid-specific Lyz2 and Ccl6 (Fig. 5c; Supplementary Figure 13). To confirm Lyz2+ myofibroblasts at the protein level, we performed single-cell western blot on unsorted cells from day 12 Sm22-Cre;tdTomato wounds and stained for LYZ and tdTomato (Fig. 5e; Supplementary Figure 14). Indeed, ~6% of tdTomato-expressing wound cells were LYZ+/tdTomato+ double-positive (77 out of 1293 cells) (Fig. 5f). We also co-stained day 12 wounds for LYZ and fibroblast markers PDGFRA (Fig. 5g) and SMA (Fig. 5h), and in both cases detected double-positive cells.

15To determine if myeloid marker-expressing fibroblasts persist in healed wounds, we performed scRNA-seq on tdTomatohi cells isolated from Sm22-Cre;tdTomato wounds at two additional time points, 15 days and 21 days PW (Supplementary Figure 15a). Combined analysis on all cells from 12 days, 15 days, and 21 days PW showed that they form four clusters, tC1 through tC4 (Supplementary Figure 15b; Supplementary Data 6). Cells from all three time points contributed to all four clusters (Supplementary Figure 16a), although relative cluster contribution was time point-dependent (Supplementary Figures 15c and 16b). Cluster tC4 was closely related to cluster fC3 from 12 days PW analysis (Supplementary Figure 15d) and also contained Cd45 and Lyz2 expressing cells (Supplementary Figure 16c). Taken together, we identified a population of myofibroblasts with myeloid features that persist in wounds over the course of repair and regeneration and potentially form via myeloid cell reprograming.

16Following our scRNA-seq findings, we noted that previous work indicated that myeloid cells can convert into myofibroblasts, but the extent of this conversion is organ and injury context-specific–. Considering that in large wounds, adipocytes regenerate predominantly from myofibroblasts, we asked to what extent myeloid cells contribute to their population. We performed wounding experiments in bone marrow transplantation (BMT) mouse models. In some BMT models, we reconstituted lethally irradiated mice with GFP-expressing cells, while in others with cells expressing lacZ under the control of lineage-specific promoters (Fig. 6a). To distinguish between hematopoietic and non-hematopoietic cells, we performed BMT using the sorted CD45+ bone marrow (BM) fraction. In CD45+ BMT models, we achieved high levels of peripheral blood chimerism, compatible to BMT mice reconstituted with whole BM (Supplementary Figure 17). We also generated BMT mice using multipotent hematopoietic stem cells (HSCs) purified based on the SLAM marker signature: Lineageneg, SCA1+, c-KIT+, CD150+, CD48neg (Supplementary Figure 17e). We transplanted between 2300 and 4400 HSCs per recipient mouse and in all cases achieved successful reconstitution of the hematopoietic lineage, confirmed by high levels of peripheral blood chimerism (Supplementary Figure 17f, g) and bone marrow fluorescence (Supplementary Figure 18a). Following wounding, healed tissue in GFP+ HSCs BMT mice showed consistently high contribution from hematopoietic lineage on 28 days PW, with many GFP+ cells surrounding neogenic hair follicles (n = 18 animals) (Fig. 7a; Supplementary Figure 19). Flow cytometry on wound tissue confirmed that long-term contribution from the hematopoietic lineage constituted ~30% at both 28 days (n = 3 biologically independent samples) and 2 months PW (n = 3 biologically independent samples) (Fig. 6b–d). In contrast, BMT mice that received GFP+ CD45neg non-hematopoietic BM fraction, which includes mesenchymal stromal cells, had no GFP+ contribution to the wound (Fig. 7b) despite showing BM repopulation (n = 11 animals; these BMT mice also received wild type whole BM support cells) (Supplementary Figure 18b).

17Next, we wounded BMT mice reconstituted with Sm22-Cre;R26R HSCs. In this model, we expected lacZ expression to mark contractile cells that derived from reconstituted hematopoietic lineage cells (but not from transplanted HSCs directly), that activated Sm22 promoter-controlled Cre recombinase. Indeed, consistent with the possibility of long-lasting hematopoietic contribution to wound myofibroblasts we observed lacZ-positive cells in wounds from Sm22-Cre;R26R HSCs BMT mice (n = 9 animals) 6 months after wounding (Fig. 6e–h). This was contrasted by the absence of lacZ-positive cells in the wounds of BMT mice reconstituted with Sm22-Cre;R26R CD45neg BM cells (n = 4 animals) (Fig. 6i). We also found lacZ-expressing cells along many blood vessels (Fig. 6e, f), suggesting that upon BMT hematopoietic lineage cells may also participate in reconstituting perivascular mural cells.

18Because fibroblast and fat lineages are closely related developmentally, we hypothesized that some hematopoietic lineage cells that initially convert into myofibroblasts might then become de novo adipocytes. Although contested,, several studies have reported that, in principle, myeloid cells can convert into adipocytes upon integration into pre-existing fat tissue–. Contribution of myeloid cells towards fat is variable and is dependent on gender and anatomical site, with visceral fat in female mice experiencing up to 10% contribution.

19We employed both BMT and non-BMT mouse models to evaluate this possibility. Indeed, wounds in GFP+ HSCs BMT mice contain many GFP+ cells that co-stained for the adipocyte marker FABP4 in areas surrounding neogenic hair follicles (Fig. 7c). By taking advantage of the fact that in mice the expression of Retn, which encodes for Resistin, is exclusive to mature fat cells, we generated Retn-lacZ HSCs BMT and Retn-lacZ CD45neg BMT models using the approach described above (Fig. 6a). In these mice, nuclear lacZ expression can be used to genetically identify functionally mature de novo adipocytes that may have derived from transplanted BM progenitors rather than from the host sources. We consistently observed formation of occasional lacZ-positive de novo adipocytes in the hair-bearing wounds of Retn-lacZ HSCs BMT mice (n = 72 animals) (Fig. 7f, g), but not in Retn-lacZ CD45neg BMT mice (n = 50 animals) (Fig. 7d). Importantly, in agreement with a previous report, we also observed integration of hematopoietic-derived lacZ-expressing adipocytes into normal fat of Retn-lacZ HSCs BMT mice (Fig. 7e).

20Next, we examined if the hematopoietic lineage in BMT mice contributes to de novo adipocytes via direct conversion rather than cell fusion. For that purpose, we generated BMT models in which we transplanted GFP+ HSCs into RFP+ hosts and R26R+/− HSCs into Fabp4-Cre hosts. In the first model, all fused cells can be identified on the basis of GFP/RFP double fluorescence. In the second model only fusion-derived de novo adipocytes contain Fabp4-Cre construct necessary for R26R reporter activation. We observed only very occasional GFP/RFP double positive cells in the hair-bearing wounds of GFP+ HSCs into RFP+ BMT mice (n = 8 animals) (Fig. 7h). LacZ-positive adipocytes were not observed in regenerated wounds and in normal fat of R26R+/− HSCs into Fabp4-Cre BMT mice (n = 13 animals) (Fig. 7i), in contrast to robust lacZ expression in fat of non-BMT Fabp4-Cre;R26R mice (Fig. 7j). These experiments confirm that hematopoietic-to-adipose lineage conversion, rather than cell fusion, takes place during de novo fat regeneration in the wound.

21We then studied wound healing in Cd45-Cre;R26R mice to verify that hematopoietic cells contribute to de novo fat physiologically and not only in the context of BMT models. In these mice, where Cre activity is restricted to the hematopoietic lineage, we observed consistent, albeit occasional formation of lacZ-positive de novo adipocytes (n = 9 animals) (Fig. 8a; Supplementary Figure 20a). Similarly, lacZ-positive adipocytes formed in the wounds of LysM-Cre;R26R mice (n = 12 animals) (Fig. 8b, c; Supplementary Figure 20b), and a portion of tdTomato+ cells, derived from 26 days PW hair-bearing wounds of LysM-Cre;tdTomato mice, differentiated into lipid-laden adipocytes when cultured under adipogenic conditions (n = 3 replicate cultures) (Fig. 8d). This suggests that hematopoietic lineage contribution to fat neogenesis is mediated via myeloid progenitors. Consistent with occasional distribution patterns of lacZ-positive adipocytes in wounds of LysM-Cre;R26R mice, LysM-Cre;Pparγ−/− mutants did not have a clear de novo fat defect (n = 20 animals) (Fig. 8e vs. 8f), characteristic of Sm22-Cre;Pparγ−/− mice described in Plikus, Guerrero-Juarez. Curiously, in both mouse models, we also observed occasional formation of neogenic hair follicles with lacZ-positive DPs and dermal sheath (Fig. 8a, b, red arrowheads; Supplementary Figure 20), suggesting that the lineage plasticity repertoire of myeloid cells during wound regeneration might extend beyond adipogenesis. Taken together, our lineage tracing studies help to establish the role of myeloid cells as one of the sources of adipogenic progenitors during wound regeneration.

22Traditionally, adult mammals are considered to have limited regenerative abilities and scarring is thought to be the default repair response. The notable exceptions to this rule are digit tip regeneration after amputation, and neogenesis of hair follicles,,, and fat in the center of large excisional wounds. Intriguingly, lineage studies reveal important differences in the regenerative strategies between these two systems. Epithelial and mesenchymal structures in the digit tip regenerate from several types of fate-restricted progenitors and no multipotent progenitors or lineage reprogramming events are observed,. In contrast, large skin wounds demonstrate broadened lineage plasticity. Although progeny of preexisting hair-fated bulge stem cells migrate into wound epidermis, they do not partake in hair follicle neogenesis. Instead, new hair follicles regenerate from non-bulge epithelial stem cells,, among other sources. Fat neogenesis is driven by lineage reprogramming of non-adipogenic wound myofibroblasts. DP neogenesis also likely relies on myofibroblast reprogramming strategy.

23Are all wound myofibroblasts identical or heterogeneous in terms of their origin, properties, and morphogenetic competence? While studying cell origin requires genetic fate mapping, and assessment of morphogenetic competence requires functional studies, scRNA-seq allows large-scale profiling of cell properties in complex tissues. Recently, this approach has been applied to studying cellular heterogeneity in the skin, including epithelial cells of mouse hair follicles,, and early wounds, and keratinocytes and fibroblasts, of human skin. Tabib, Morse identified two major populations of human skin fibroblasts, characterized by co-expression of SFRP2+/DPP4+ and FMO1+/LSP1+ markers. These further subdivide into several subpopulations. Philippeos, Telerman identified five fibroblast populations: upper and lower dermal fibroblasts, pericytes, and two uncharacterized populations. scRNA-seq has also been used to study heterogeneity of disease-associated fibroblasts in the synovial tissue upon rheumatoid arthritis, and in the lung upon bleomycin-induced fibrosis.

24Here, we studied fibroblast heterogeneity in the mouse model for wound-induced regeneration at 12 days PW when wound re-epithelialization is completed and preceding hair follicle neogenesis. At later time points 15 days and 21 days PW, we anticipated a decrease in SMA expression and change in fibroblast types. We show that wound fibroblasts can be broadly classified into two major populations on the basis of their transcription factor signatures and PDGF receptor expression patterns. Pdgfrahigh fibroblasts, which dominate in the wound, further subdivided into upper and lower dermal fractions, based on their localization and gene expression. Pdgfralow cells constituted nearly a quarter of all wound fibroblasts and localized throughout wound dermis. Lineage tracing will be required to determine which fibroblast population preferentially contributes toward newly regenerating hair follicle mesenchyme and adipocytes. Prominent additional heterogeneity exists within both populations, and pseudotime and RNA velocity analyses predict three distinct fibroblast-to-myofibroblast differentiation trajectories. Future work should aim to understand the functional significance of this heterogeneity in the context of regeneration and to lineage trace wound fibroblast origin to distinct skin fibroblasts in the unwounded skin.

25BM-derived progenitors, including myeloid cells, endothelial progenitors, and circulating mesenchymal stem cells can contribute new stromal cells toward injured tissues in various organs. Scar tissue in heart following myocardial infarction, cornea following keratectomy, and lung in pulmonary fibrosis contains BM-derived collagen-producing myofibroblasts. In skin, studies document BM giving rise to fibroblasts at the injury sites–, while some studies report this contribution to be minimal,. Such a discrepancy is likely attributed to several factors, including the type and extent of injury and experimental timing. The contribution from BM progenitors toward repairing tissues was shown to increase with the extent of injury,,, yet most of the previous studies were performed on small wounds. Also, many previous studies failed to evaluate long-term BM contribution, yet BM-derived fibroblasts in the scar were shown to peak after 22 days.

26Our data from large excisional wounds shows that the contribution from myeloid cells to wound fibroblasts is small yet significant, between 6% and 11.3%, depending on the assessment method. scRNA-seq analysis revealed an intermediate state of wound cells, displaying dual gene expression signature, enriched for contractile and mesenchymal extracellular matrix genes, as well as myeloid cell surface antigens and inflammatory cytokines. RNA velocity analysis showed that wound cells progress through this intermediate state in a myeloid-to-fibroblast direction. We also showed that at least a portion of these cells can convert into de novo adipocytes around neogenic hair follicles. Previously, adipogenic conversion of myeloid cells has been shown both in vitro, and in vivo in major adipose depots, such as in inguinal fat,; although some studies contested these observations,. Most recently, ~10% of adipocytes were shown to form from hematopoietic lineage in human subjects undergoing BMT treatment, while in another human BMT study, up to 35% of adipocytes were traced to transplanted BM sources.

27Overall, our findings illustrate the dynamic nature of fibroblast identities during wound healing, and the powerful wound-induced plasticity of myeloid lineage cells. scRNA-seq analysis infers several pathways fibroblasts follow during wound healing, including contractile and regenerative. The lineage relationship between myeloid cells and adipogenesis suggests further characterization of the factors influencing plasticity and fate switching could uncover potential therapeutic approaches to treating wounds and scars.

28The following mouse models were used in this study: C57BL/6J (JAX stock 000664), Retn-lacZ, Sm22-Cre (JAX stock 004746), Cd45-Cre, LysM-Cre (JAX stock 004781), Fabp4-Cre (JAX stock 005069), Pparγflox (JAX stock 004584), R26R (JAX stock 003474), tdTomato (JAX stock 007909), GFP (JAX stock 004353), RFP (JAX stock 006051), and Rag1−/− (JAX stock 002216).

29All animal experiments were carried out in accordance with the guidelines of the Institutional Animal Care and Use Committee of the University of Pennsylvania and University of California, Irvine. Animals were anesthetized with isoflurane, hairs were clipped, skin site was disinfected and a single full thickness excisional wound was created on their dorsum using scissors (squared s = 1.5 cm)–. Both male and female mice between 6 and 8 weeks of age were used. Following wounding, all animals were housed individually. Wounds were let to heal by secondary intention. No wound dressing was applied. Animals were used as biological replicates.

30Bone marrow cells were flushed from the long bones with Dulbecco’s Modified Eagle’s Medium supplemented with 5% heat-inactivated calf serum (Gibco). To obtain a single cell suspension, cells were filtered through a 45 µm nylon screen. Cells were blocked with anti-rat and anti-mouse IgG (Sigma) for 15 min. For isolation of HSCs, cells were then stained with CD45 (APC-Cy7, BD Pharmingen), CD150 (APC, Biolegend), CD48 (Pacific Blue, Biolegend), SCA1 (PE-Cy5.5, Invitrogen), and c-KIT (PE-Cy7, eBioscience), as well as PE-conjugated lineage markers including CD19, B220, CD3, CD4, CD8, GR-1, MAC-1, NK1.1, and TER119 (eBioscience). Cells were resuspended in DAPI to discriminate live from dead cells. HSCs were sorted by isolating Lineageneg, SCA1+, c-KIT+, CD150+, CD48neg cells using a BD FACSAria II flow cytometer (BD Biosciences). Approximately 150,000–2,000,000 CD45+ cells or 2300–4400 HSCs, or ~150,000–2,000,000 CD45neg cells were resuspended in phosphate buffered saline with 100,000–200,000 whole bone marrow cells for reconstitution and injected into lethally irradiated isogenic wild type host mice. For Sm22-Cre;R26R BMT experiments, cells were transplanted into immunodeficient Rag1−/− mice. Donor chimerism was assessed by measuring GFP or CD45.1 vs. CD45.2 levels by FACS in the peripheral blood 8 weeks after transplantation.

31To detect lacZ activity, freshly isolated wound tissue samples were incubated with X-gal reagent in lacZ staining buffer,. Samples were post-fixed in 4% PFA.

32PFA-fixed wound tissue samples were pre-incubated in propylene glycol and then stained with OilRedO buffer for 20 min. Samples were then washed in propylene glycol and stored in 0.05% aqueous solution of sodium azide.

33Tissues were fixed in 4% PFA, dehydrated, paraffin embedded, and sectioned at 7–12 µm thickness. Frozen tissues were sectioned at 12 µm thickness. Immunostaining was performed both on frozen and paraffin sections. Heat-based antigen retrieval was performed when necessary. Tissue sections were blocked in either 3% BSA, 3% Donkey serum, or 3% mouse IgG from Mouse on Mouse (MOM) kit (Vector Labs). The primary antibodies used were goat anti-FABP4 (1:200; R&D Systems), mouse anti-SMA (1:500; R&D Systems), rabbit anti-RGS5 (1:250; Proteintech), rabbit anti-LYZ (1:250; Abcam), rabbit anti-CRABP1 (1:800; Cell Signaling), goat anti-CYP26B1 (1:50; Novus Biologicals); rabbit anti-G0S2 (1:100; Proteintech), mouse anti-MEST (1:250; Proteintech), goat anti-PDGFRA (1:350; R&D Systems), chicken anti-β-Gal (1:500; Abcam), and rabbit anti-PLIN (1:800; Cell Signaling). Tissue sections were visualized with Nikon T inverted (Nikon) and Olympus FluoView confocal (Olympus) microscopes.

34Primary cells were isolated from 26 days PW pooled wound tissues of LysM-Cre;tdTomato mice. Tissues were micro-dissected, minced, and incubated in a Collagenase IV solution for 60 min with constant rotation. Post-incubation, cells were filtered through a 40 µm nylon screen, centrifuged, washed, resuspended in DMEM supplemented with 10% FBS, antibiotics, and antifungals and cultured at 37 °C in a water-jacketed incubator with 5% CO2 output. Once cells reached confluency, they were switched to adipocyte differentiation media (Cell Solutions) for 72 h, upon which they were switched to adipocyte maintenance media (Cell Solutions) for additional 7 days. Cells were incubated with BODIPY (Thermo) as per manufacturer’s directions and visualized using a Nikon T inverted (Nikon) microscope.

35Pooled wound tissues (n = 12 animals) were collected from Sm22-Cre;tdTomato mice on 12 days PW. Tissues were micro-dissected, minced, and incubated in a Dispase II/Collagenase IV/Liberase solution for 60 min with constant rotation. Post-incubation, cell aggregates were mechanically dissociated using GentleMACS (MACS). Single cell suspensions were treated with 1X RBC lysis buffer, washed, and re-suspended in 0.04% UltraPure BSA (Biolegend). Dead cells were removed using the MS columns of the Dead Cell Removal Kit (MACS) as per manufacturer’s directions. Live cells were resuspended in 0.04% UltraPure BSA and counted using the automated cell counter Countess (Thermo). GEM generation, barcoding, post GEM-RT cleanup, cDNA amplification, and cDNA library construction were performed using Single Cell 3′ v2 chemistry (10X Genomics). Library quality control metrics were as follows: 550.19 pg/µl with an average size ~454 bps. Libraries were sequenced on an Illumina HiSeq4000 platform (Illumina) (one lane, 100 PE). Cell counting, suspension, GEM generation, barcoding, post GEM-RT cleanup, cDNA amplification, library preparation, quality control, and sequencing were performed at the Genomics High Throughput Sequencing facility at the University of California, Irvine.

36Transcripts were mapped to the mm10 reference genome (GRCm38.91) using Cell Ranger Version 2.1.0. The web summary of Cell Ranger statistics on cells, sequencing and mapping metrics pre-quality control assessment for downstream bioinformatics analyses were as follows—cell metrics: ~22,322 estimated number of cells, ~84.2% fraction reads in cells, ~13,819 mean reads per cell, ~1101 median genes per cell, ~19,070 total genes detected, ~2448 median UMI counts per cell; sequencing metrics: ~308,471,010 total number of reads, ~98.5% valid barcodes; mapping metrics: ~90.4% reads mapped to genome, ~85.5% reads mapped confidently to genome, ~65.9% reads mapped confidently to transcriptome.

37For downstream analyses after initial Cell Ranger metric assessment, low-quality cells were removed to eliminate cell-specific biases. Quality control metrics included keeping cells displaying < 8000 UMI/cell and < 2500 genes/cell, and no more than 8% mitochondrial gene expression. Post-quality control, 21,819 cells remained for downstream bioinformatic analyses.

38Clustering of cells was performed using the Seurat R package. In brief, digital gene expression matrices were column-normalized and log-transformed. To identify cell clusters, principle component analysis (PCA) was first performed on the list of highly variable genes. A set of highly variable genes was then identified by binning the average expression of all genes into evenly sized groups and computing the median dispersion (variance divided by the mean) in each bin. Genes were selected for inclusion in PCA with an average expression > 0.01 and dispersion > 1.0. To identify significant PCs, we used the Jackstraw method (JackStraw function in Seurat). The top 40 PCs were used for clustering with the Louvain modularity-based community detection algorithm to generate cell clusters (FindClusters function, 13 clusters with resolution = 0.45). Marker genes were determined with p-value < 0.01 and log(fold-change) > 0.25 by performing differential gene expression analysis between the clusters using the likelihood-ratio test. To present high dimensional data in two-dimensional space, we performed t-SNE analysis using the results of PCA with significant PCs as input. For the subclustering analysis on all wound fibroblasts, we followed the same procedures as above. Briefly, the top 35 PCs were used for clustering and 12 subclusters were obtained with resolution = 0.6.

39We used cell cycle-related genes, including a previously defined core set of 43 G1/S and 54 G2/M genes. For each cell, a cell cycle phase (G1, S, G2/M) was assigned based on its expression of G1/S or G2/M phase genes using the scoring strategy described in CellCycleScoring function in Seurat. Cells in each cell cycle state were also quantified using the which.cells function in Seurat.

40We subclustered wound fibroblasts based on known fibroblast-specific markers. We used Monocle 2 to infer the pseudotime trajectories of these cells, followed by scEpath to identify pseudotime-dependent gene expression changes. To model gene expression changes in pseudotime, scEpath first divides the pseudotime into 10 equally spaced bins. Then the expression of each gene in each bin is estimated by the trimean of the expressions of this gene across all the cells located in this bin. Furthermore, scEpath smoothens the average expression of each gene using cubic regression splines. To determine the pseudotime-dependent genes that are significantly changed, we compared the standard deviation of the observed smoothed expressions with a set of similarly permuted expressions by randomly permuting the cell order (1000 permutations in our analyses). We considered all genes with a standard deviation > 0.5 and a Bonferroni-corrected p-value below a significance level α = 0.01 to be pseudotime dependent. To analyze pseudotime-dependent TFs, we used genes that are annotated in the Animal Transcription Factor Database (AnimalTFDB 2.0). For the pseudotime analysis of Lyz2+ wound myeloid cells and Lyz2+ fibroblasts, we followed the same procedures.

41RNA velocity was calculated on the basis of spliced and unspliced transcript reads as previously reported. Based on velocyto pipeline, annotation of spliced and unspliced reads was performed using the Python script velocyto.py on the Cell Ranger output folder. Only cells that were part of the pseudotime were considered in analyses. PCA analysis was performed with Pagoda2 with spliced expression matrix as input, and cell-to-cell distance matrix was calculated using Euclidean distance based on the top 40 principal components. RNA velocity was estimated using gene-relative model with k-nearest neighbor cell pooling (k = 100). Velocity fields were projected onto the pseudotime space produced by Monocle 2. Parameter n sight, which defines the size of the neighborhood used for projecting the velocity, was 3500. For other parameters, we used the standard R implementation of velocyto with default settings.

42Pooled wound tissues (n = 2–3 animals) were collected from Sm22-Cre;tdTomato mice on 12, 15, and 21 days PW. Wound tissue was micro-dissected and disaggregated into single cells with Dispase II (Sigma) and Collagenase IV (Sigma) at 37 °C with constant rotation. Single cell fraction was stained with Zombie Violet (1:1000; BioLegend) and FACS-sorted as Zombieneg;tdTomatohi with a BD FACSAria II flow cytometer (BD Biosciences). Pre-sorted, viable tdTomatohi single cells were re-suspended in DMEM supplemented with 10% FBS, antibiotics and antifungals, diluted with suspension reagent for attribution of optimal buoyancy, and loaded onto a large 17–25 µm 96-well microfluidic IFC (Fluidigm) for single cell capture in the automated C1 system for single-cell genomics (Fluidigm). Capture efficiency was assessed using bright field and fluorescent microscopy. Only cells captured singly (singlets) per micro-well were considered for downstream analyses. Double (doublets) and multiple (multiplets) cells captured per well were excluded. Lysis, RT, and cDNA pre-amplification were performed in loco (protocol 1.773×) with ultra-low input RNA reagents as suggested per manufacturer (Clontech). RNA spike-in controls were omitted. cDNA concentrations were estimated using Qubit 2.0 (Thermo) and cDNAs with concentration ≥ 1.0 ng/μl were used for downstream library preparation. Libraries were amplified using the Nextera XT v2 Index Kit (Illumina). Quality control on multiplexed libraries was performed using the Agilent Bioanalyzer and quantification was performed using KAPA for Illumina Sequencing Platforms (Illumina). Multiplexed libraries were sequenced as paired-end on an Illumina Next-Seq500 platform (Illumina).

43Demultiplexed, paired-end FASTQs were aligned to the mouse genome (mm10/gencode.Mv13) using Bowtie (version 1.0.0) with the following standard parameters: rsem-prepare-reference –bowtie –gtf and quantified using the RNA-seq by Expectation-Maximization algorithm (RSEM) (version 1.2.31) with the following standard parameters: rsem-calculate-expression -p $CORES –paired-end. Samples displaying ≥ 159,000 aligned reads were considered for downstream quality control filtering.

44Cells with a total number of expressed genes > 11,000, a proportion of counts in mitochondrial genes > 15%, and genes expressed in < 3 cells were excluded for downstream analyses. Fifty-six cells were removed from all replicates, leading to a total of 288 cells (116 cells at 12 days PW, 100 cells at 15 days PW, and 72 cells at 21 days PW) for downstream analyses. To adjust systematic variation in the relationship between transcript-specific expression and sequencing depth, we normalized our data using the quantile regression-based method SCnorm. Samples were then natural log-transformed and added a pseudo-count of 1 followed by batch correction for each time point using the Bayesian-based method ComBat using the ComBat function from the sva R package. The batch corrected data were used for downstream analyses.

45Single-cell western blot assays were performed using the ProteinSimple Milo platform with the standard scWest Kit as per manufacturer’s protocol. scWest chips were rehydrated and loaded with cells at a concentration of 1 × 105 cells in 1 mL suspension buffer. Doublet/multiplet capture rate in scWest chip microwells was determined with light microscopy (~1.3%, established from > 1000 microwells). scWest chips were lysed for 10 s and immediately followed by electrophoresis for 80 s at 240 V. Protein immobilization was achieved with UV light exposure for a total of 4 min. scWest chips were sequentially probed with primary and secondary antibodies for 1 h. Primary antibodies used were goat anti-mCherry (1:10, SicGen) and rabbit anti-LYZ (1:20, Abcam). Secondary antibodies used were donkey anti-goat Alexa Fluor 647 (1:40, Thermo) and donkey anti-rabbit Alexa Fluor 555 (1:40, Thermo). Slides were washed, centrifuge-dried, and imaged with the GenePix 4000B Microarray Scanner (Molecular Devices). Data was analyzed using Scout Software (ProteinSimple) and ImageJ (NIH). Debris, artifacts, and false positive signals were manually excluded during data analyses.

46Further information on experimental design is available in the Nature Research Reporting Summary linked to this Article.