Single‐cell transcriptomic analysis of small and large wounds reveals the distinct spatial organization of regenerative fibroblasts

Abstract Wound‐induced hair follicle neogenesis (WIHN) has been an important model to study hair follicle regeneration during wound repair. However, the cellular and molecular components of the dermis that make large wounds more regenerative are not fully understood. Here, we compare and contrast recently published scRNA‐seq data of small scarring wounds to wounds that regenerate in hope to elucidate the role of fibroblasts lineages in WIHN. Our analysis revealed an over‐representation of the newly identified upper wound fibroblasts in regenerative wound conditions, which express the retinoic acid binding protein Crabp1. This regenerative cell type shares a similar gene signature to the murine papillary fibroblast lineage, which are necessary to support hair follicle morphogenesis and homeostasis. RNA velocity analysis comparing scarring and regenerating wounds revealed the divergent trajectories towards upper and lower wound fibroblasts and that the upper populations were closely associated with the specialized dermal papilla. We also provide analyses and explanation reconciling the inconsistency between the histological lineage tracing and the scRNA‐seq data from recent reports investigating large wounds. Finally, we performed a computational test to map the spatial location of upper wound fibroblasts in large wounds which revealed that upper peripheral fibroblasts might harbour equivalent regenerative competence as those in the centre. Overall, our scRNA‐seq reanalysis combining multiple samples suggests that upper wound fibroblasts are required for hair follicle regeneration and that papillary fibroblasts may migrate from the wound periphery to the centre during wound re‐epithelialization. Moreover, data from this publication are made available on our searchable web resource: https://skinregeneration.org/.


| INTRODUC TI ON
Fibroblast heterogeneity is an important and emerging area of research in skin biology. 1 These cells perform diverse functions in the dermis and hypodermis of the skin. 1,2 The fibroblast subtypes that are closest to the epidermis are called papillary fibroblasts and are easy to observe histologically in neonatal mouse skin and in human skin. 3,4 In mice and in humans, there is an age-related reduction of the papillary region which may play a role in the decreased wound-healing abilities associated with aging. 5,6 Furthermore, the papillary fibroblast lineage in murine skin has been shown to be required to regenerate hair follicles as they can differentiate into dermal papilla. 3 Dermal papillae are another fibroblast subtype that resides at the base of the hair follicle and are necessary for their formation and maintenance. 7 Consequently, the proportion of papillary fibroblasts within a wound can affect the regenerative outcome. Reticular fibroblast subtypes reside in the reticular dermis, hypodermis and fascia. 3,8 They produce the extra-cellular matrix (ECM) of the dermis and possess a pool of progenitors for adipocytes. 9,10 Reticular fibroblasts do not differentiate into dermal papilla and produce the ECM of the scar. 3,8 Harnessing the functions of different fibroblast lineages during wound healing is an important step to achieving functional skin regeneration.
The development of scRNA-seq has provided unprecedented insights into the cellular heterogeneity of different tissues, including the ability to compare fibroblasts between wounds that scar and those that regenerate with hair follicles. [11][12][13][14] These comparisons of regenerating and scarring wounds using scRNA-seq have mapped diverse populations of fibroblast within the wound beds but have not fully understood how they can differentiate into dermal papilla to support hair follicle reformation. However, two recent publications investigating the differences between scarring and regenerating wounds by scRNA-seq now provide a foundational data set to investigate the ability to regenerate. 15,16 The genetic markers for cellular clusters in scRNA-seq experiments have recently been translated into a spatial map by histological analysis. 13,17 These experiments have allowed for the use of fibroblast markers to identify the spatial location of a fibroblast subtype histologically or computationally. For example, Crabp1 expression in wound fibroblasts marks the cells closest to the upper region of the wound bed. 13 These are referred to as upper wound fibroblasts. Fibroblasts that express Mest and Plac8 are represented in the lower regions of the wound and are referred to as lower wound fibroblasts. 13,17 The association of scRNA-seq gene expression and spatial location will provide additional insights into the molecular and cellular mechanisms of fibroblasts heterogeneity during the wound-healing process.
It has been recently shown that Hypermethylated in cancer (Hic1) transgenic CreERT2 mice can label both papillary and reticular fibroblasts at varying efficiencies. 15 Histological analysis of lineage traced Hic1-labelled cells revealed high contribution to regenerating dermal papilla in WIHN. However, scRNA-seq analysis indicated low contribution of Hic1CreER-tdTomato cells to regeneration-competent fibroblasts.
Here, we perform a detailed analysis to address this inconsistency and provide an explanation for these results. We also reanalysed three recently published scRNA-seq data of scarring and regenerative wounds to compare different fibroblast subtypes. In both conditions, we observed the diverging differentiation trajectories between upper and lower wound fibroblasts. Moreover, our analysis highlighted the striking difference in the relative abundance of upper wound fibroblasts between regenerating and scarring wounds. Finally, we performed a computational test to investigate the regenerative competence of upper wound fibroblasts, which unexpectedly revealed that the wound periphery also possesses the equal regenerative potential after the re-epithelization of large wounds.

| Abbasi et al., 2020 data
Generation of scRNA-seq data (GSE108677) is described. 15 In short, cells were isolated from large full-thickness wounds of Hic1-tdTomato mice treated with tamoxifen at P4-5. In addition, collected cells were flow sorted for viability and tdTomato expression. In this manuscript, we utilized data from four different 10× Genomics libraries: (1) large wound centre at D14 and tdTomato positive, (2) large wound centre at D14 and tdTomato negative, (3) large wound periphery at D14 and tdTomato positive and (4) large wound periphery at D14 and tdTomato negative. scRNA-seq BAM files are downloaded from SRA database. BAM files are processed by Velocyto.py (0.17.15) to generate loom files for downstream RNA velocity analysis. A custom reference genome was built using the cellranger mkref by adding the tdTomato transgene sequence to FASTA and GTF files of a prebuilt mm10 reference.

| Lim et al., 2018 data
We obtained FASTQ files from the single-cell RNA-sequencing data (GSE112671) recently published in Lim et al, 11

| Phan et al., 2020 data
Generation of scRNA-seq data (GSE153596) was described in Phan et al., 2020. Cells were isolated from small (2 mm) wounds at 7 days postwound from wild-type P2 and P21 mice (n = 3). A total of 6 libraries were generated using 10× Genomics 3′ Single-Cell Gene expression V2 kit. Libraries were sequenced on Illumina HiSeq 4000 with one full lane per sample. Cellranger v3.0.2 pipeline was used to process fastq files and generated output files. The output files were run through Velocyto.py (0.17.15) to produce .loom files for downstream and RNA velocity analysis. In this manuscript, we used all three libraries from P2 wounds and one library (P21_3) from P21 wounds.

| scRNA-seq data analysis
Loom files were analysed using SCANPY (1.5.1) 18  For scVelo analysis, we performed both stochastic and dynamical models, which resulted in similar outcomes. In this manuscript, the Stochastic model was presented.

| Comparative scRNA-seq analysis of scarring and regenerating wounds
To investigate differences between scarring and regenerative wound repair, we performed a scRNA-seq analysis from fibroblasts of two scarring wounds and two types of wounds that regenerate hair follicles ( We were able to identify a distinct population that expressed high levels of Tgfbi. This Tgfbi high population has overlapping expression with other lower lineage markings. Our reanalysis further confirmed the heterogeneity of fibroblasts detected in both scarring and regenerating wounds.

| Upper and lower wound fibroblasts have distinct differentiation trajectories in RNA Velocity analysis
To further investigate the dynamic nature of fibroblast heterogeneity in wounds, we performed RNA velocity analysis. RNA velocity estimates the future state of individual cells by calculating the ratio of unspliced and spliced mRNAs to predict lineage trajectories within scRNA-seq data. 19,22 We utilized scVelo, an integrated package for SCANPY, to perform this analysis on all four conditions (Methods).
Upper wound fibroblasts, which express Crabp1, are located closer to the wound epithelium and are thought to be the source of regenerative potential of WIHN. Our velocity analysis revealed that in scarring and regenerative wounds, upper wound fibroblast trajectories project away from lower wound fibroblasts (Figure 2A Table S1). The upper Crabp1 fibroblasts were stimulated to differentiate into dermal papilla cells with Shh activation ( Figure S2A). We infer that upper and lower wound fibroblasts have distinct differentiation trajectories in the wound.

| scRNA-seq analysis indicates that Hic1 reticular lineage fibroblasts may not contribute to the dermal papilla of WIHN
These scRNA-seq data are the first to contain definitive neogenic dermal papilla from large regenerating wounds. 12 Step 2a-b). Single-cell preparations were generated of LWC and LWP, which were then flow sorted for live cells. This strategy generated four 10× Genomic libraries which were as follows: LWC td-Tomato negative (LWC 14dpw Neg), LWC tdTomato positive (LWC 14dpw Pos), LWP tdTomato negative (LWP 14dpw Neg) and LWP tdTomato positive (LWP 14dpw Pos) ( Figure 3A Step 4). This elegant scRNA-seq approach provides us with an avenue to computationally test the potentiality of Hic1 lineage contribution to WIHN dermal papilla, and to ascertain which location and fibroblast subtype is the most similar to the regenerative dermal papilla of WIHN.
We investigated tdTomato expression within the scRNA-seq data by realigning all libraries to a reference genome containing tdTomato sequence for reanalysis ( Figure 3B). The newly generated UMAP contained all upper and lower fibroblast populations identified by Crabp1, Mest, Plac8 and Tgfbi expression ( Figure 3B and data not shown). Next, we overlayed the four conditions defined by wound location and flow-sorted tdTomato detection onto the new UMAP ( Figure 3C). As in Figure 3B,C,E,F, the regenerating dermal papilla cluster consisted mainly (85%) of fibroblasts that were negative for tdTomato by flow cytometry ( Figure 3C). In fact, only 3% of the neogenic dermal papilla originated from LWC 14dpw Pos flow-sorted library. We also verified that regenerating dermal papilla from LWC 14dpw Neg was also negative for tdTomato transcripts ( Figure 3D).  with the Hic1CreERT2-tdTomato transgenic mouse model has the potential to label the papillary and reticular dermis, with a preference for reticular fibroblasts. Furthermore, the papillary and reticular dermis has different roles in supporting WIHN. Papillary fibroblasts in wounds can differentiate into dermal papillae, while reticular fibroblasts do not directly support hair follicle regeneration.
Consequently, different labelling efficiencies of these fibroblast subtypes will influence the lineage tracing outcome in large wounds. In As a consequence, reporter genes may succumb to dropout in DP with droplet-based single-cell barcoding platforms ( Figure 3H and Figure S3). Uninjured skin processed using 10xV2 confirms tdTomato is similarly downregulated and dropped out as only 17% of DP are tdTomato+ve (compared to 65% of tdTomato+ve DS and 76% of tdTomato+ve interfollicular fibroblasts; Figure S1N,O). By flow sorting-enrichment of tdTomato-high subset (to limit contamination from tdTomato-ve fraction), Abbasi et al. may have inadvertently excluded tdTomato-low DP. This provides an explanation for tdTo-mato+ upper LWP identified as the direct source of tdTomato-low/ negative neogenic DP ( Figure 3F). Furthermore, we present two additional experiments that independently show similar results from a reanalysis of flow-sorted alpha-SMACreERT2-RosaeYFP+ve hair follicle mesenchyme from uninjured skin ( Figure S3). 26  To further study this sub-population, we quantified the contribution of each library to each of the three Leiden clusters identified as upper wound fibroblasts ( Figure 4C). Cluster 1 mainly composed of LWC tdTomato-positive cells (79%) and was projected to become cluster 2, where 75% of the cells were LWP tdTomato positive cells.
However, the RNA velocity trajectory of upper wound fibroblasts was disconnected from the neogenic DP when cells from all four libraries were included ( Figure 4A-C). Given the over-representation of peripheral fibroblasts in cluster 2 (which was closest to the neogenic DP), we hypothesize that there might be distinct differences between upper wound fibroblasts from wound periphery and centre. We next performed a computational test using different conditions/libraries of large wound fibroblasts ( Figure 4D,E). Since the neogenic DP were almost exclusively derived from the LWC tdTomato-negative cells, we wanted to evaluate their relationship with cells from LWC tdTomato positive and LWP tdTomato positive. Our analysis suggested that the LWC tdTomato positive was not associated with the neogenic DP, as seen by the proximity of cells on UMAP ( Figure 4D). RNA velocity also did not predict a trajectory of LWC-positive upper wound fibroblasts towards the neogenic DP (Table S3). In contrast, the upper wound fibroblasts from the periphery clustered closely with the neogenic DP, and RNA velocity analysis predicted a continuous trajectory from peripheral upper wound fibroblasts towards the centre neogenic DP (Table S2). Additionally, we have also identified a population of LWC tdTomato-negative cells that were closely associated with the neogenic DP ( Figure 4D).
These fibroblasts are most likely the source of the regenerating DP at the wound centre. Interestingly, the LWP upper wound fibroblasts also directly overlapped with this population, indicating their similarities in transcriptomic profiles ( Figure 4E). It is important to note that inferences drawn from gene expression analysis may not capture the multi-tiered regulatory dynamics that were highlighted in Abbasi et al., which may ultimately determine fibroblast function.
Nevertheless, based on this analysis, we inferred that the population of peripheral upper wound fibroblasts exhibit competence to support regeneration of DP similar to the upper wound centre fibroblasts.

| D ISCUSS I ON
The rediscovery of hair follicle regeneration in large wounds in murine skin (a.k.a. WIHN) combined with the application of modern transgenic technology has provided for an essential model to dissect the cellular and molecular mechanisms required to support functional skin regeneration. [10][11][12][13]27 The WIHN model system has been utilized in studies to reveal the importance of epidermal stem cell differentiation during wound repair and regeneration. However, one of the critical questions regarding the origins of the regenerative dermal papilla and fibroblasts that support WIHN remains Convergence of signalling pathways leading to changes in regulatory network activity ultimately drives the acquisition of regenerative competence. The fact that regenerative propensity of large wounds (indicated by the number of neogenic HFs) remains amenable to pharmacological and genetic perturbations [28][29][30][31] supports the notion that fibroblasts in wound periphery harbour a latent degree of regenerative competence.
Our reanalysis of scRNA-seq data from small and large wounds adds to the evidence for the role of restricted fibroblast lineages in wound repair and regeneration. 3,10,32 Likewise, the reanalysed data support the hypothesis that the migration of papillary fibroblasts is a critical component for skin regeneration. 3,33 In conclusion, our in-depth reanalysis of large wound scRNA-seq data reveals that the field of fibroblast heterogeneity is an understudied area of research with great potential to uncover critical knowledge for understanding skin development, wound healing and regeneration.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest.

AUTH O R CO NTR I B UTI O N S
QMP processed, analysed, interpreted the data and co-wrote manuscript. SS and JB provided the data, assisted in analysing and interpreting the data and co-wrote the manuscript. RRD conceived of the project, assisted in analysing and interpreting the data, and co-wrote the manuscript.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.   Table S1. List of 362 velocity genes. Table S2. Periphery test velocity genes Figure 4E. Table S3. Center test velocity genes from Figure 4D.