Single-Cell Analysis Reveals a Hair Follicle Dermal Niche Molecular Differentiation Trajectory that Begins Prior to Morphogenesis.

Introduction:
How complex tissue structures initially emerge is a question in development that has been examined across diverse appendages (Dhouailly, 1975; Hardy, 1992; Jiang et al., 1999; Olivera-Martinez et al., 2004). The molecular and cellular events that precede the emergence of a morphological structure remain poorly defined, as conventional tools lack the ability to discriminate molecular differences over a large number of transcripts between individual cells. Consequently, efforts to define progenitor populations and the signals that lead to cell fate specification have fundamentally relied on morphologic segregation.
The hair follicle (HF) is one of the most tractable models to study appendage development, because the morphological stages of HF development are well-defined, and many of the molecular and cellular events that chronicle these stages have been characterized. HFs first appear as patterned epithelial thickenings (placodes) that are associated with underlying dermal condensates (DCs, Figure 1A) (Hardy, 1992; Millar, 2002; Paus et al., 1999; Xin et al., 2016). The DC is a cluster of specialized dermal cells that forms coordinately with the HF epithelium. DC cells express Sox2, one of the earliest DC markers of mouse primary and secondary HF types (Driskell et al., 2009; Sennett et al., 2015). Later stages of HF morphogenesis are marked by growth and organized differentiation of the HF epithelium as it envelops the DC, which matures into the dermal papilla. The dermal papilla serves as the permanent signaling center required for the cyclical regeneration of the adult HF.
Many of the molecular events essential for epithelial placode initiation, including epithelial Wnt/β-catenin and Ectodysplasin A receptor (Edar) activation were previously shown (Andl et al., 2002; Huelsken et al., 2001; Zhang et al., 2009). By contrast, the events that lead to DC cell specification and DC formation remain sparsely delineated. Previous studies, including those using bulk RNA-sequencing and live imaging, have characterized some of the signals and cell behaviors of differentiated DC cells once placodes are histologically apparent (Biggs et al., 2018; Glover et al., 2017; Sennett et al., 2015; Tsai et al., 2014). While these studies uncovered mechanisms by which differentiated DC cells condense into a morphological DC, the events that precede DC cell differentiation are largely unknown. Two studies suggested that DC cells are specified from apparently equivalent upper dermal cells in response to nearby epithelial placode signals (Biggs et al., 2018; Glover et al., 2017). These differentiated DC cells exit the cell cycle, migrate, and undergo cell shape changes to form a dense cluster. However, the potential molecular and cellular changes that precede DC cell differentiation were difficult to examine a priori. For instance, studies have shown that differentiated DC cells are quiescent (Biggs et al., 2018; Wessells, 1965; Wessells and Roessner, 1965). A lesser known fact is that cells surrounding quiescent DCs were noted to show a high rate of proliferation by Wessells in 1965 who proposed that these peri-DC cells divide to give progeny to the expanding DC. Since then, this decades-old theory has remained largely untested. The role of cell proliferation by a latent pre-DC population during DC formation and/or expansion was not addressed in these studies. As the DC is the central organizer required for HF growth and differentiation, our understanding of the mechanisms that direct HF development has stalled at this initial stage.
Based on seminal studies, including heterospecific recombination experiments, a prevailing theory posits that HF induction is controlled by the dermis, which determines appendage size and distribution (Dhouailly, 1975; Hardy, 1992; Jiang et al., 1999; Kratochwil et al., 1996; Olivera-Martinez et al., 2004). The molecular identity of this “first dermal signal” is unclear but appears to be unpatterned and is thought to act on the epidermis through paracrine mechanisms to promote placode specification. The molecular identity of this dermal signal was suggested by chick and mouse embryonic studies, which showed that the upper dermis is composed of Wnt/β-catenin-activated cells, leading to the hypothesis that transduction of Wnt signal in the dermis may be required to transmit the “first dermal signal” essential for HF induction (Chang et al., 2004; Chen et al., 2012; Fu and Hsu, 2013; Noramly et al., 1999). In support of this, early ablation of the transducer of Wnt signaling, β-catenin, in mouse dermal cells or loss of epidermal Wnt ligand secretion required for dermal Wnt/β-catenin signaling, resulted in a complete absence of HF placodes and DCs (Chen et al., 2012; Fu and Hsu, 2013). These key studies demonstrated that dermal Wnt/β-catenin signal transduction is a candidate pathway preceding DC formation and that could also be instructive for DC cell fate specification. However, the mechanism by which dermal Wnt signaling affects DC initiation remains unknown, as knowledge of the molecular and cellular events that lead to DC cell fate are unknown.
Single-cell RNA sequencing (scRNA-seq) technology offers the ability to distinguish fine molecular differences between individual cells that gives rise to DC cells. With scRNA-seq data, a diffusion map technique can be employed to codify the transcriptomic similarities and differences between dermal cells and to determine the diffusion distance between cells (Coifman et al., 2005). The resulting diffusion map assesses how cells relate to each other transcriptionally and infers patterns of differentiation that are based on changes across the transcriptome and not limited to single genes or pathways. This allows us to make unbiased predictions about molecular and cellular changes that precede DC differentiation. Using scRNA-seq data from two time points, before and once HFs have initiated, we infer a pattern, or trajectory, of transcriptional states through which differentiated DC cells pass. Coupled with in vivo studies, we show that DC cells originate from Wnt-activated dermal progenitors and that transduction of Wnt signal in the dermis is necessary to enter an intermediate transcriptional stage prior to differentiation. Based on scRNA-seq predictions, our in vivo proliferation analysis indicates that DC cells are quiescent progeny of a selectively proliferative population present prior to morphogenesis.

Results:

Dermal condensate cells specify a molecular differentiation trajectory that initiates prior to morphogenesis
To examine molecular patterns in developing skin before (E13.5) and at the time of HF initiation (E14.5), we microdissected dorsolateral skin from mouse embryos at both time points and obtained single-cell suspensions (Figure 1B). Droplet-based scRNA-seq was used to sequence the transcriptomes of live single cells FACS sorted from E13.5 and E14.5 skin replicates (Figure S1A, S1B) (Zheng et al., 2017). We first defined different skin populations, using t-distributed stochastic neighbor embedding dimension reduction (t-SNE) and unsupervised clustering. Dermal and keratinocyte populations were identified using Col1a1 for dermal populations and Krt10, Krt14 for keratinocyte populations. These populations were then subjected to downstream analyses (Figures 1C, S1A, S1B).
To examine the molecular changes in the dermis associated with DC differentiation, we constructed diffusion maps in which cells that are more molecularly similar are positioned closer and cells with more divergent molecular signatures are spaced further apart (Coifman et al., 2005) to infer the transcriptional states through which DC cells pass. Diffusion maps were constructed using the dermal scRNA-seq data sets obtained from E13.5 and E14.5 skin replicates (Figure 1D). We used the early DC marker, Sox2, to identify differentiated DC cells (Figures 1E, S1C) and observed that Sox2+ DC cells lie at the terminus of a trajectory that aligns with eigenvector 2, indicating that this specific trajectory represents cells at different transcriptional states that are arranged by how molecularly similar they are to differentiated DC cells. Cells closer to the Sox2+ terminus were more similar to DC cells, while those nearer to the origin were largely undifferentiated. Notably, this DC-specified trajectory contained both E13.5 dermal cells, which are more closely concentrated toward the origin, and E14.5 dermal cells that exclusively made up the terminus of the trajectory (Figures 1D, 1E).
We next assessed transcriptional changes that occur along the DC-specific trajectory by ordering cells as they occur along this trajectory, which we defined as pseudo-order. We first removed other trajectories, ordered the cells along the DC trajectory and smoothed gene expression values using a generalized additive model (Figures S1E-S1H) (Hastie, 1990). In addition to Sox2, other genes known to be expressed by DC cells (e.g. Bmp4, Ptch1) were also upregulated in cells closer to the terminus of the trajectory (Figure 1F). In contrast to Sox2, which showed a near binary expression pattern at the end of the DC pseudo-order, expression of some other DC markers was detected at both E13.5, prior to DC morphogenesis and at E14.5 in cells lacking Sox2 transcripts (Figure S2A). Hierarchical clustering of genes across the DC pseudo-order showed that expression of DC-associated markers and cell cycle regulation genes correlated with this DC differentiation trajectory (Figure S2C and Table S1).
To identify candidate signaling pathways that correlated with progression on the trajectory and that could be traced to an early progenitor population, we used KEGG pathway analysis. Pathways were prioritized according to the most significant correlation (p-value ≤ 0.05) with Wnt/β-catenin signaling correlating the most highly with the DC trajectory (Figures S2D, S2E). Visualization of the canonical Wnt target genes, Lef1 and Twist2, on the diffusion map showed that both genes exhibited increasing expression along the trajectory (Chen et al., 2012). This pattern is apparent at E13.5 and becomes more established at E14.5 (Figures 1G, S2F). Analysis of Wnt pathway genes along the DC pseudo-order showed that expression of Wnt target genes increased with increasing pseudo-order (Figure 1H). The direct transcriptional target and negative regulator of Wnt/β-catenin signaling, Dkk1, also paralleled this trend until the acquisition of Sox2 transcripts when Dkk1 expression dropped acutely in a manner that inversely correlated with Sox2 expression (Figures 1H, S1H) (Chamorro et al., 2005; Niida et al., 2004). Notably Lox, a previously reported pan-fibroblast marker, negatively correlated with expression of most DC markers over the pseudo-order (Figures 1F, S2A) (Sennett et al., 2015).
In situ hybridization (ISH) was performed to correlate the inferred gene expression patterns with spatial location. The dermal marker, Col1a1, showed expression throughout the dermis with decreased expression in the DC, while Lef1 and Axin2 expression were confined to the upper dermis and epidermis at both time points and were upregulated in the DC and HF epithelium at E14.5 (Figures 1I, S2B). Bmp4 and Dkk1 transcripts were detected similarly in the upper dermis at E13.5. Consistent with previous work, while Bmp4 was expressed in the DC at E14.5, Dkk1 appeared to be excluded from the DC in a peri-DC distribution (Figures 1I, S2B) (Andl et al., 2002; Chen et al., 2012). Lox-expressing dermal cells were largely distributed in the dermis beneath the Wnt-activated zone at both time points and were excluded from the DC at E14.5 (Figures 1I, S2B). Together, these data suggest that DC cells specify a molecular trajectory that is characterized by a continuum of transcriptional states that initiates by E13.5 and appears to localize spatially to the upper dermis at E13.5 and the DC/peri-DC region at E14.5.

The number of Wnt/β-catenin activated dermal progenitors regulates dermal condensate size
Loss of Wnt/β-catenin signaling in dermal progenitors results in an absence of HF placodes and DCs (Chen et al., 2012; Fu and Hsu, 2013). Based on the diffusion map, DC cells pass through progressive transcriptional states characterized by increasing levels of Wnt/β-catenin activation (Figures 1G, S2F). Consistent with previous work, we found that Sox2+ DCs at E14.5 were Wnt-activated, assessed by Axin2LacZ and by Lef1 staining, while Wnt-activated dermal cells were limited to the upper dermis at E13.5 (Figures S2G, 1I, S2B) (DasGupta and Fuchs, 1999; Zhang et al., 2009). Genetic lineage tracing of Wnt-activated dermal cells was done by administering tamoxifen (Tam) at E12.5 to Axin2CreER;Rosa-LSL-mTom/mGFP (Axin2CreER;mTmG) embryos (Figure 2A), resulting in ~52% recombined membrane GFP expressing (mGFP+) dermal cells at E13.5 that were predominantly confined to the upper dermis (Figure S3A). Three-dimensional (3D) volumetric quantification of mGFP+ DC cells in skin whole mounts taken from E14.5, E15.5 and P0 Axin2CreER;mTmG embryos showed that all DCs of all mouse HF types (primary, secondary, and tertiary) were highly enriched for mGFP+ cells (Figures 2A, S3A), showing that Wnt-activated upper dermal progenitors present prior to DC morphogenesis selectively contribute to DCs at later time points. At the same time, mGFP+ cells also contributed to the upper interfollicular dermis (IFD), but the proportion of mGFP+ IFD cells decreased over time, suggesting that mGFP-negative cells (e.g. lower dermal cells) were able to contribute to the upper IFD. Consistent with this, we observed that the proliferation rate of mGFP+ IFD cells decreased over time after E13.5 (Figure S3B) and correlated with a decreased proportion of Wnt-active cells in the IFD (Figure S2G).
We next asked if the lack of DCs caused by ablation of dermal β-catenin could not only be due to a paracrine effect on the epidermis but could also reflect a requirement for transduction of Wnt signal by dermal cells to acquire DC cell fate. To begin to address this, we used an inducible Cre recombinase driven by the Pdgfrα promoter to temporally ablate β-catenin expression in a varying number of dermal cells prior to HF initiation (PdgfrαaCreER;βcatfl/fl). Administration of a single high dose of Tam (50 μg/gm) at E11.5 (PdgfrαCreER(Hi);βcatfl/fl) resulted in a marked loss of Wnt-activated cells in the dermis by E13.5 assessed by Axin2 ISH (21% of control; Figures 2B, 2C). Tissue-specific recombination in the dermis and not epidermis was confirmed using a Cre reporter (Figure S3G). By E14.5, patterned HF placodes associated with Sox2+ DCs were observed in control embryonic skin (Figures 2D, S3G). As expected, mutant PdgfrαCreER(Hi);βcatfl/fl embryos showed a complete absence of HFs at E14.5, with no patterned expression of the HF epithelial markers, Krt17 or P-cadherin. The HF stem cell transcription factor, Sox9, which is also known to be expressed by DC cells (Sennett et al., 2015), was also not detected in PdgfrαCreER(Hi);βcatfl/fl embryos (Figures 2D, S3C). Likewise, the DC marker, Sox2, was not detected in the dermis of mutant embryos. HFs were not detected in PdgfrαCreER(Hi);βcatfl/fl embryos examined later at E16.5, showing that HF induction was not simply delayed (Figure 2E). To attain a lower level of Cre recombination, we administered a single low dose of Tam (30 μg/gm weight) to pregnant mice at the same time point (PdgfrαCreER(Lo);βcatfl/fl), which resulted in a greater proportion of Wnt-activated (Axin2+) dermal cells than that seen with high doses of tam (52% of control) (Figures 2F, 2G). At E14.5, PdgfrαCreER(Lo);βcatfl/fl embryos formed HFs of normal distribution and density, but that were significantly reduced in size (Figures 2H, 2I, S3G). Volumetric quantification of the number of Krt17+ epithelial cells per HF and Sox2+ cells per DC showed that both populations were proportionately reduced in number and correlated with decreased HF epithelial volume (Figures 2J, S3E). Asymmetric cell divisions occur early in HF development resulting in basal committed (P-cadHi) and suprabasal stem cell (Sox9+) daughter cells (Ouspenskaia et al., 2016; Xu et al., 2015). Despite the marked reduction in the number of Krt17+ epithelial cells per mutant HF, the proportion of Sox9+ HF stem cells and P-cadHi committed basal epithelial cells within mutant HFs remained similar to control HFs (Figure 2J).
In corollary experiments, varying expression of an activated form of β-catenin in the dermis (PdgfrαCreER(Hi);βcatfl(EX3/+) and PdgfrαCreER(Lo);βcatfl(EX3/+)) showed a Tam dose-dependent increase in the proportion of Wnt-activated dermal cells, which correlated with an increase in DC and HF epithelial size (Figures S3H-S3O). Enlarged hair germs showed increased numbers of Krt17+ HF epithelial and Sox2+ DC cells per HF (Figures S3I, S3K, S3M, S3O) (Chen et al., 2012) and increased HF epithelial volume (Figure S3J, S3N). For both E14.5 mutant PdgfrαCreER(Hi);βcatfl/(EX3/+) enlarged and PdgfrαCreER(Lo);βcatfl/fl small HFs, epithelial cell proliferation was similar to controls (Figures S3K, S3D), indicating that differences in HF or DC size were not due to changes in the rate of HF epithelial proliferation. Previous work in mice showed that loss of dermal Wnt/β-catenin resulted in decreased dermal cell proliferation at E14.5. However, upper dermal cell density was normal despite an absence of HFs, suggesting that dermal cell density per se does not control HF induction (Chen et al., 2012). Consistent with this, we found that dermal proliferation and upper dermal cell density were both decreased at E14.5 in PdgfrαCreER(Hi);βcatfl/fl and PdgfrαCreER(Lo);βcatfl/fl skin (Figure S3F), despite observing HFs in the latter. Together, these data indicate that the number or density of specifically Wnt-activated dermal progenitors regulates the size of the DC and HF epithelium.

Loss of dermal β-catenin results in impaired progression along a dermal condensate differentiation trajectory
Previous work suggested that upper dermal Wnt signaling, in response to epithelial, but not dermal, Wnt ligands may serve as the “first dermal signal” that induces HF placode specification through unknown tissue-tissue interactions (Chang et al., 2004; Chen et al., 2012; Fu and Hsu, 2013). Here, the requirement for dermal Wnt activation in DC initiation may be secondary to signals from the placode required for DC cell fate (Chen et al., 2012). However, based on our combined results, it was possible that Wnt/β-catenin signaling may also be required for dermal cells to acquire DC cell fate. In this scenario, only cells capable of transducing Wnt/β-catenin signal could acquire DC cell fate. To address this, we examined DC cells within the small HFs of PdgfrαCreER(Lo);βcatfl/fl mosaically recombined mutant embryos. We obtained scRNA-seq data from E13.5 and E14.5 PdgfdrαCreER(Lo);βcatfl/fl mosaically recombined mutant and control skin cells. Control and mutant scRNA-seq datasets were merged, batch-corrected, and keratinocyte and dermal populations were analyzed (Figures S4A-S4E). As expected, t-SNE plots of the merged mutant and control datasets showed age-dependent expansion of differentiated (Krt10) populations in both mutant and control samples (Figures S6C, S6D). HF epithelial populations were also identified in both mutant and control keratinocytes at E14.5. Consistent with our in vivo results, there was a smaller proportion of HF epithelial cells observed in the mutant keratinocyte population (13.2% control vs. 7.8% mutant) (Figure S6E). Consistent with mosaic ablation of β-catenin, we found that the mutant dermal population had a diminished Wnt-activated population at E13.5 and E14.5 (Figures 3A, 3B, S4F) and fewer putative DC cells as defined by expression of Sox2 (1.1% in mutants vs. 3.1% in controls; p <0.00001 chi-square test) (Figures 3A, 3B). However, DC cells from control and mutant embryos clustered similarly instead of segregating into divergent populations, suggesting that DC cells from controls and mutants are similar molecularly and that mutant cells were excluded from the DC population (Figure 3A). We next analyzed the Wnt target genes, Axin2 and Lef1, within the Sox2+ DC population from mutant and control dermal cells and found that although the DC population from mutant embryos was diminished in size, they were composed of Wnt-activated dermal cells similar to the control DC population (Figures 3C, 3D). Analysis of skin sections from PdgfrαCreER(Lo);βcatfl/fl embryos showed that all DC cells from mutant small HFs expressed markers of Wnt activation, confirming that DCs from mutant PdgfrαCreER(Lo);βcatfl/fl skin were composed of a non-mosaic population of Wnt-activated (wildtype) dermal cells (Figures 3E, 3F). These data indicate that transduction of Wnt/β-signal is required for dermal cells to acquire DC cell fate.
We next generated diffusion maps from both control and mosaic mutant dermal scRNA-seq datasets and identified the Sox2+ DC differentiation trajectory (Figures S5B, S5C). We found that the diffusion map of the merged data corresponded with the pattern observed in the control dermis (Figures S5D, S5E). The dermal population from mutant embryos also appeared to follow the same molecular trajectory as that of control cells but was more sparsely distributed toward the end of the DC-specific trajectory compared to controls (Figures S5E). This suggested that mutant dermal cells were impaired in their progression along the DC trajectory. Accordingly, density plot analysis showed an accumulation of mutant dermal progenitors prior to an intermediate phase of the DC trajectory compared to control dermal cells, indicating a relative block in differentiation at this stage (Figure 4A). Mutant dermal cells were further behind in this transcriptional progression, demonstrating lower values on the DC transcriptional trajectory than those of control dermal cells at E13.5 and E14.5 (p < 0.00001, permutation test). Specifically, Wnt-activated cells appeared to represent an intermediate cell population on the DC molecular trajectory (Figure 4B), and those cells from a mosaic population of mutant and wildtype cells that could be seen within this intermediate stage and at the terminus were Wnt-activated (wildtype). Consistent with this, we observed a significantly reduced number of PdgfrαCreER(Lo);βcatfl/fl dermal (Col1a1+) cells that expressed markers found within this intermediate stage, Bmp4 and Dkk1 (Figure 4C). Analysis of genes that were differentially expressed across the point on the trajectory when mutant dermal cells exhibit stalled progression showed that genes associated with extracellular matrix organization and cell adhesion (Table S3) were upregulated across this point. These data show that dermal β-catenin is required to progress into an intermediate phase on the DC-associated trajectory that is marked by increasing levels of Wnt/β-catenin signaling, cell proliferation, and expression of other DC markers before DC cell differentiation and DC morphogenesis. Our data suggest that the number of Wnt-activated cells that can enter this phase directly correlates with the number of DC cells.
Based on the hypothesis that the “first dermal signal” may act through paracrine mechanisms to regulate epithelial placode specification, we next analyzed the keratinocyte scRNA-seq datasets to assess potential non-cell autonomous placode-associated epidermal changes that could be affected by loss of dermal Wnt/β-catenin signaling (Figure S7). Notably, we found that despite the presence of a dermal DC trajectory at E13.5 (Figures 1E, S5E), we could not detect an E13.5 keratinocyte subpopulation that was specifically associated with a putative HF placode (pre-placode) cell type (Figures 5, S7D). Based on the wildtype keratinocyte diffusion maps, E13.5 keratinocytes appeared to be molecularly distinct from E14.5 basal keratinocytes, and E14.5 basal keratinocytes appeared to represent an intermediate transcriptional state between E13.5 keratinocytes and E14.5 HF placode cells (Figures 5, S7). Analysis of the placode differentiation trajectory, which aligns with eigenvector 5, showed that canonical placode markers were upregulated at the end of the E14.5 pseudo-order (Figures 5, S7). However, we were unable to detect placode-specific genes, including Wnt ligand genes, that showed an expression pattern that significantly correlated with this eigenvector at E13.5 (Figures S7E, S7F). Although these data do not exclude the possibility that dermal Wnt/β-catenin signaling functions to induce epithelial competency for placode cell fate, we could not detect a placode-specific trajectory at E13.5 when we detect a stall in a DC-associated trajectory by β-catenin mutant dermal cells.

Differentiated dermal condensate cells are quiescent progeny of highly proliferative progenitors
Previous work showed that differentiated DC cells are quiescent (Biggs et al., 2018; Wessells and Roessner, 1965). Consistent with this, measurement of proliferation by EdU incorporation confirmed that DCs were composed of largely quiescent cells (Figure 6A). However, examination of cell cycle phase of scRNA-seq cell populations showed that while almost all DC cells were in G0/G1 phase, Wnt-activated non-DC cells proliferate more than other dermal cells (Figure 6B). Analysis across the wildtype E14.5 DC pseudo-order showed that the proportion of proliferating cells progressively increased until the vast majority of cells were proliferating and that cells further along the pseudo-order began to enter G0/G1 phase prior to Sox2 expression. Upon Sox2 expression, virtually all cells were quiescent and in G0/G1 phase (Figure 6C).
According to the DC trajectory, cells show a marked increase in their rate of proliferation prior to quiescence and DC differentiation. Prior to HF morphogenesis, E13.5 dermal cells along the DC pseudo-order also demonstrated increasing rates of proliferation that correlated with progression on the E13.5 DC trajectory. Those cells at the terminus of the E13.5 trajectory did not yet show Sox2 expression (Figure 6D) but expressed other DC-associated markers (e.g. Dkk1, Lef1), which are expressed by cells closer to the E14.5 Sox2+ terminus (Figure 1H). Notably, we found that Cdkn1a, a gene that promotes cell cycle arrest and known to be expressed by Sox2+ differentiated DC cells, was also upregulated at the end of the pseudo-order at both time points while expression of the cyclin-dependent kinase activator, CyclinA2, correlated with the fraction of dividing cells at E13.5 and E14.5 (Figures 6D, 6I) (Biggs et al., 2018; Huh et al., 2013; Sennett et al., 2015).
Based on these data, we hypothesized that DC cells may represent progeny of highly proliferative progenitors. To address this, we performed EdU nucleotide pulse-chase experiments to trace proliferating dermal cells until DCs formed. EdU was administered to pregnant mice at E13.5. Embryos were examined either 3 hours after EdU administration (E13.5), one day later at E14.5, or two days later at E15.5 when DCs had expanded (Figure 6E). In contrast to the upper IFD, which showed dilution of EdU over time, DCs instead showed an increased proportion of EdU+ cells at E14.5 (52.1%±4.7) compared to upper IFD cells at all time points, including at E13.5 (39.8%±0.9) when cells were incorporating EdU (Figure 6E). This indicated that although DC cells are quiescent, they are progeny preferentially of cells that were proliferating at ~E13.5 when EdU was administered. Only occasional Sox2+EdU+ double-positive cells co-expressed Ki-67 (7.30±4.5 %Ki-67+; Figure 6F), and using the in vivo Fucci2 cell cycle reporter, virtually all traced EdU+ cells in the morphologic DC were in G0/G1 phase (Figure 6G). Traced EdU+ cells were still present in the expanded DCs at E15.5 (48.28±4.23 %EdU+), suggesting that these cells remained quiescent DC cells. Notably, we did not observe a statistically significant difference between the proportion of EdU+ DC cells at E14.5 and E15.5 despite an increased number of Sox2+ DC cells seen at E15.5 (Figures 6E, 6H), suggesting that those cells that were incorporated into the DC between E14.5 and E15.5 were quiescent by E14.5.
Examination of genes expressed as cells showed increasing proliferation followed by cells that were entering G0/G1 phase indicated that Wnt target genes, including Dkk1, showed increasing expression (Figure 6J). As Dkk1 is largely excluded from the Sox2+ DC population, we used this marker to help spatially locate putative DC progenitors, which were distributed in the upper dermis at E13.5 and in the peri-DC zone at E14.5 (Figures 1I, 6J). Quantification of proliferation showed that a higher rate of proliferation is shown by Dkk1+ cells in the upper dermis than by other upper dermal cells at E13.5 (54.0%±14.1 Dkk1+ vs. 22.9%±3.7 Dkk1−) and at E14.5 by Dkk1+ cells within both the peri-DC zone (57.4%±5.7 Dkk1+) and upper IFD (29.5%±3.9 Dkk1+ vs. 11.1%±2.1 Dkk1−) (Figure 6L). Consistent with this, co-detection of Lef1 and EdU in the peri-DC region showed that Wnt-activated Lef1+ peri-DC cells showed a higher rate of proliferation (29.6±0.4 %EdU+) compared to Lef1+ (17.9%±2.0) or Lef1− (6.3%±1.3) upper IFD cells (Figures 6K, 6L). Additionally, based on the DC trajectory, expression of CyclinA2 is upregulated in a manner that correlates with the proliferating fraction, while expression of Cdkn1a is upregulated when the fraction of quiescent cells increases. Detection of these markers showed that Cdkn1a is localized to Sox2GFP+ DC cells whereas CyclinA2+ cells were found within the IFD at a higher frequency in cells surrounding the Sox2GFP+ DC (53.9%±0.6 CyclinA2+ peri-DC vs. 36.2%±1.8 IFD) (Figure 6M). Taken together with our EdU pulse-chase analysis, these results suggest that dermal cells show a high rate of proliferation before entry into G0/G1 phase and subsequent DC differentiation. Notably, events on the E14.5 DC-associated trajectory, including increasing rates of proliferation followed by Cdkn1a expression, are detected at E13.5 prior to Sox2 expression or DC morphogenesis. Our combined scRNA-seq and in situ data implicate Wnt-activated Dkk1+ upper dermal and peri-DC cells as a potential selective population of DC progenitors that divide to give quiescent progeny that become Sox2+ DC cells (Figure 6J).

Dermal condensate size is dependent on active dermal proliferation
Similar to previous work, we found a dose-dependent decrease in the proportion of proliferating dermal cells in E14.5 PdgfrαCreER(Lo);βcatfl/fl and PdgfrαCreER(Hi);βcatfl/fl embryos (39.8%±3.3 control vs. 29.5%±3.7 PdgfrαCreER(Lo);βcatfl/fl and 21.3%±5.4 PdgfrαCreER(Hi);βcatfl/fl) (Figure 7A). As PdgfrαCreER(Lo);βcatfl/fl mutant embryos show smaller DCs and decreased proliferation, it was possible that the diminished DC size observed was due to a lack of progenitor cell divisions required for generating DC progeny. To address this hypothesis, we treated E13.5 PDGFRαH2BGFP skin explants with mitomycin C (MMC) to inhibit proliferation prior to apparent DC morphogenesis and then allowed explants to grow in culture without MMC for an additional 40-44 hours. Quantification of proliferation showed a dose-dependent decrease in dermal cell proliferation in MMC-treated explants (6.8%±3.2 EdU+ and 1.2%±0.5 Edu+) compared to vehicle-treated explants (22.0%±2.5 EdU+) (Figure 7B). Skin explants after 44-48 hours of culture showed that a normal pattern of HFs formed in MMC-treated explants. However, there was a dose-dependent decrease in DC and HF epithelial size in MMC-treated explants compared to controls (Figures 7B, 7C). We found that the proportion of dermal cells expressing Wnt target genes was similar between conditions (Figure 7B). Therefore, while Wnt/β-catenin signaling was maintained, inhibiting proliferation phenocopied the PdgfrαCreER(Lo);βcatfl/fl mutant, suggesting that Wnt/β-catenin signaling may function, in part, by promoting the proliferation of DC progenitors.

Discussion:

scRNA-seq analysis infers a transcriptional pattern of DC differentiation that is detected prior to morphogenesis
The dermo-epidermal interactions that govern appendage formation have been studied across different species, revealing several conserved principles of appendage induction and pattern formation (Dhouailly, 1975; Jiang et al., 2004; Olivera-Martinez et al., 2004). Despite the availability of tools to interrogate specific HF populations (Myung and Ito, 2012), investigating steps that precede morphologic changes has proved challenging, as we lack tools to resolve differences between histologically indistinguishable cells. This study evolved from recognition of this major limitation and offers an initial unbiased inference into this process.
Here, we used scRNA-seq to distinguish different transcriptional states present in embryonic skin, employing a diffusion map method to codify the relationship between these different states. This method infers molecular pathways of differentiation based on the expectation that during any differentiation process, not all transcriptional states are stable nor do they necessarily represent committed states (Coifman et al., 2005). By arranging these transcriptional states by similarity, inferred intermediate states through which a differentiated cell has passed are used to construct a trajectory. In an accurate diffusion map, cells cannot “jump” along the path, but necessarily have to pass through these intermediate transcriptional states, assuming that large transcriptional changes are not instantaneous.
This inferred transcriptional path is directionless, and an individual cell can potentially move in any direction: toward a differentiated state, toward a progenitor-like state, or even stall. Thus, while the interpretation is that mature differentiated cells have passed through “earlier” intermediate states, intermediates may not necessarily progress to reach a differentiated state. The diffusion map also does not in itself provide information about the identity or requirement for factors that characterize transcriptional changes over a differentiation path. However, using multiple time points, this method provides a unique opportunity to predict the molecular and cellular changes that occur during a differentiation process as well as the spatial location of these events to ultimately guide further in vivo investigations. Combining this method with in vivo analyses, we examined early transcriptional events within the epidermis and dermis to delineate an inferred pathway of DC differentiation.
We took advantage of previous work that showed that dermal Wnt/β-catenin signaling, in response to epithelial Wnt ligands, is essential for HF initiation and DC formation (Chen et al., 2012; Fu and Hsu, 2013). These studies provided one of the few molecular leads that could be leveraged to investigate mechanisms of DC induction. Here, we employed unbiased scRNA-seq analysis across two developmental time points, before and during DC morphogenesis, to infer a pattern of DC molecular differentiation, exploiting dermal Wnt signaling as a tool to cross-examine this pathway.
Using scRNA-seq data, we identified both differentiated DC and placode populations at E14.5. Using the diffusion map method, we could distinguish transcriptional differences between more undifferentiated dermal cells, including differences in expression levels of Wnt target genes. Our analysis correlates the DC pseudo-order with heterogeneous expression of Wnt target genes, and the varying expression levels of Wnt target genes also correlate with each other, suggesting that they represent part of a biological process rather than a random mix of biologic or experimental noise. We find that these changes are present at E13.5 and hypothesize that they represent early transcriptional events that precede DC cell differentiation. Consistent with this, we found that most cells at the terminus of the E13.5 DC-associated trajectory had adopted characteristic cell cycle phase changes known to occur during DC differentiation, including upregulation of Cdkn1a. In addition, the E14.5 DC-associated trajectory contained cells within the E13.5 trajectory, suggesting that intermediate transcriptional states that lie proximal to a differentiated DC state are also represented by those observed at E13.5.
By contrast, our scRNA-seq analysis was unable to detect robust molecular changes that correspond with a placode-specific differentiation axis at E13.5. Interestingly, E13.5 keratinocytes appear to represent a novel transcriptional state from that of E14.5 basal keratinocytes, and E14.5 keratinocytes appear to give rise to both suprabasal and placode populations. These inferred data suggest that dermal cells begin to molecularly resemble a DC cell prior to detection of a placode-specific differentiation pattern. Nevertheless, it is possible that molecular changes at E13.5 that are not specifically associated with a placode transcriptional state are required for placode specification.
As we do not detect Sox2+ cells until placode cells are evident at E14.5 and as a committed pre-DC state that precedes Sox2 expression has not yet been shown, our study does not address the hierarchy of placode and DC cell specification. Accordingly, placode-driven gene expression changes appear to correlate with events toward the end of the DC differentiation trajectory when Wnt target genes and other known placode-dependent genes (e.g. Ptch1) are highly upregulated. Further, recent work shows that Fgf20, which is expressed by placode cells, is necessary for DC morphogenesis and Sox2 expression, suggesting that placode specification precedes DC specification (Biggs et al., 2018; Huh et al., 2013). At the same time, it is plausible that a committed pre-DC population may exist prior to dermal Sox2 and epithelial Fgf20 expression, and placode signals are required to establish or maintain a differentiated Sox2+ DC population. Upon terminal DC cell differentiation, Sox2 is expressed and changes associated with DC morphogenesis such as cell migration and maintenance of quiescence are dependent upon placode signals. Nevertheless, although we detect upregulation of Cdkn1a at the end of the E13.5 DC trajectory and show that a selectively proliferative population, also reflected by the terminus of the E13.5 trajectory, gives rise to quiescent differentiated DC cells, further studies are needed to determine the existence of either a patterned intermediate population or a committed intermediate dermal subpopulation. To this end, our analysis can serve as a molecular framework to guide future investigations.

Dermal progenitors require transduction of Wnt/β-catenin signal to progress along the differentiation trajectory
Classical heterospecific recombination experiments demonstrated a critical role of the dermis in appendage induction. Studies in mouse and chick suggest that upper dermal Wnt/β-catenin signaling may serve as the “first dermal signal” that mediates appendage induction through paracrine mechanisms to induce placode specification. Consistent with this, previous work showed that upper dermal Wnt signaling, in response to secreted epithelial Wnt ligands, is required for placode and DC induction (Chen et al., 2012; Fu and Hsu, 2013). Our scRNA-seq and in vivo genetic data show that transduction of Wnt signal in dermal cells does not solely function to induce placode cell fate, but is also required for a dermal cell to acquire DC cell fate. Our scRNA-seq analysis further infers that dermal Wnt/β-catenin signaling is required to enter an intermediate phase along the DC trajectory in which cells undergo changes marked by increasing expression of DC markers and increasing proliferation followed by quiescence and differentiation. Notably, loss of dermal β-catenin did not induce a novel transcriptional state, but instead, resulted in a relative block in differentiation and stalled progression along the DC transcriptional axis. At the same time, forced activation of Wnt signaling in dermal cells was not sufficient to drive expression of DC markers throughout the dermis at E13.5 or E14.5, showing that dermal Wnt signaling is required, but not sufficient, for DC cell differentiation (Chen et al., 2012). It is possible that dermal Wnt signaling plays an additional downstream role that results in enlarged DCs, directly or indirectly, following expression of stabilized β-catenin. For instance, as dermal Wnt signaling is essential for dermal cell proliferation, and as enlarged DCs are quiescent, the enlarged DCs may reflect increased proliferation of pre-DC cells.

Dermal condensate cells are immediate quiescent progeny of a proliferative population
Our inferred pattern of DC differentiation shows that dermal cells undergo increasing rates of proliferation just before entering G0/G1 phase and acquiring Sox2 expression. These cells also express Dkk1 and other Wnt target genes. Given that the mean residence time of EdU is 26 minutes (Cheraghali et al., 1994), our experiments indicate that Sox2+ quiescent DCs at E14.5 preferentially incorporated cells that were replicating DNA at approximately E13.5. As virtually all Sox2+EdU+ label-retaining cells were observed to be in G0/G1 phase and remained within DCs at E15.5, these Sox2+ label-retaining cells are quiescent DC cells. This entry into quiescence may occur at least at or shortly after E13.5, as E14.5 DCs and expanding DCs from E15.5 embryos pulsed at E13.5 show the same frequency of EdU+ label-retaining cells. Our scRNA-seq and in situ studies suggest that these proliferative cells and their quiescent progeny that lack Sox2 expression are located in the peri-DC zone at E14.5 and revisits the decades-old postulation that highly proliferative peri-DC cells divide to give progeny to the expanding DC (Wessells and Roessner, 1965). Further studies are needed to determine the fate of those highly proliferative peri-DC cells.
According to both our scRNAseq and pulse-chase data, cell cycle exit occurs prior to expression of Sox2. The coordination of cell cycle phase and differentiation has been studied across different tissues, providing evidence that cells are more susceptible to external signals that induce differentiation during G1 phase (Ruijtenberg and van den Heuvel, 2016). We hypothesize that exit into G0/G1 phase may be necessary for DC cell differentiation. Recent work showed that Edar-induced epithelial Fgf20 expression is required but not sufficient for Cdkn1a expression by DC cells (Biggs et al., 2018). In contrast to our findings, this study also suggested that DC differentiation as defined by Sox2 expression precedes entry into G0/G1 phase, at least during the initial stages of DC morphogenesis. The signals that regulate the transition between proliferation and quiescence during DC cell differentiation remain unclear.

An inferred differentiation path correlates with spatial location
Our analysis suggests that the inferred pathway of DC differentiation that occurs at E14.5, such as Sox2 and Ptch1 expression, are both spatially and temporally associated with the presence of a placode and that factors from the placode determine definitive steps of DC differentiation. Accordingly, dermal cells that express these markers are spatially near to the placode epithelium. Similarly, at E13.5 when placodes are undetectable, the factors that influence the inferred dermal molecular changes could correlate with a morphogen gradient that emanates from the overlying epidermis. One candidate morphogen secreted by the epidermis is Wnt ligand. It was previously observed that the upper dermis is composed of cells that appear grossly to show low levels of Wnt signaling and low levels of expression of some DC markers. Although it was recognized that expression of these markers is likely heterogeneous in the upper dermis, quantitative resolution of this heterogeneity between cells was lacking (Chen et al., 2012; Noramly et al., 1999; Reddy et al., 2001; Zhang et al., 2009). Our scRNA-seq data demonstrates heterogeneous levels of Wnt signal activation in the dermis, which may reflect such a morphogen gradient. However, we were not able to detect a spatial gradient of Wnt signal response on a per cell basis in the upper dermis based on markers used in this study. Additionally, while our genetic data suggest that β-catenin-ablated upper dermal cells are transcriptionally similar to lower dermal cells at E13.5, forced activation of β-catenin is not sufficient to induce DC markers such as Bmp4 and Dkk1 (data not shown) (Chen et al., 2012). This suggests that a morphogen would be either acting through multiple signaling mechanisms or that there is at least one additional factor governing early transcriptional events prior to morphogenesis. Investigations into these potential signals may bring new insight into the molecular signals and cell populations that lead to DC cell commitment and DC formation.
Our study highlights the utility of unbiased single-cell approaches to overcome some of the challenges to delineating the mechanisms of appendage induction, owing to its range and sensitivity to detect differences between individual cells. More broadly, this innovative approach can help guide investigations into the development of other appendages through unbiased inference of progenitor populations and signals that are otherwise difficult to postulate and that were previously limited by histological distinction.

STAR Methods:

Contact for Reagent and Resource Sharing
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Peggy Myung (peggy.myung@yale.edu).

Experimental Model and Subject Details
No statistical methods were used to predetermine sample size. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

Mice
PDGFRαCreER mice (Collins et al., 2011) were bred to β-catfl/fl (Brault et al., 2001) and β-catfl(Ex3/+)) (Harada et al., 1999) mice to generate PDGFRαCreER;β-catfl/fl and PDGFRαCreER;β-catfl(Ex3/+) mice respectively. Axin2CreER (van Amerongen et al., 2012) mice were bred to Gt(ROSA)26Sortm4(ACTB-tdTomato–eGFP)Luo/J (mTmG) (Muzumdar et al., 2007). Reporter lines used to visualize Cre-recombination Gt(ROSA)26Sortm4(ACTB-tdTomato–eGFP)Luo/J (mTmG), Axin2LacZ (Lustig et al., 2002), PDGFRαH2BGFP (Hamilton et al., 2003) Sox2eGFP (Arnold et al., 2011) mice were obtained from Jackson Laboratory. R26Fucci2aR reporter mice were used to visualize the cell cycle phases (Mort et al., 2014). A random population of both male and female mice were used for all experiments, as it was not possible to sex the animals due to developmental age. All procedures involving animal subjects were performed under the approval of the Institutional Animal Care and Use Committee of the Yale School of Medicine.

Method Details

Tamoxifen induction of mice
Embryos were staged as days post coitum, with embryonic day E0.5 considered as noon of the day a vaginal plug was detected after overnight mating. For low-dose experiments, PDGFRαCreER(Lo);β-catfl/fl;mTmG and PDGFRαCreER(Lo);β-catfl(Ex3)/+ mice were given a single dose of Tamoxifen dissolved in corn oil (20mg/ml, Sigma) at E11.5 by oral gavage at 30 μg/gm and 15 μg/gm body weight respectively. For high dose experiments PDGFRαCreER(Hi);β-catfl/fl;mTmG and PDGFRαCreER(Hi);β-catfl(Ex3/+) mice were given a single dose at 50 μg/gm and 30 μg/gm body weight respectively. For Axin2CreER;mTmG lineage tracing experiments mice were administered a single dose of tamoxifen on E12.5 at 50 μg/gm.

EdU Incorporation Assay
To assess active proliferation, EdU was administered to pregnant mice intraperitoneally (25 μg/gm) and embryos were harvested after 1.5 hour. For pulse-chase experiments, 25 μg/gm EdU was administered to pregnant mice at E13.5 and embryos were harvested either after 3 hours, 24 hours (E14.5) or 48 hours (E15.5). For MMC treated culture experiments, 10μM EdU was added to the media for 1 hour before fixing the explants in PFA. EdU incorporation was assessed using the Click-it EdU Imaging kit Alexa 555 or Alexa 488 (Life technologies, c10338) according to manufacturer’s instructions. Briefly, skin explants were treated with a mixture of 1X Click-iT reaction buffer, CuSO4, Alexa Fluor azide dye and 1X reaction buffer additive, all provided in the Click-it EdU Imaging Kit, for 30 minutes at RT before washing in PBS.

Histology
Whole embryos were either cryoembedded in Tissue-Tek OCT Compound (Sakura Finetek) or 10% formalin-fixed paraffin embedded (FFPE) and used for histological analysis. FFPE skin was sectioned at 5 μm thickness. After deparaffinization, sections were subjected to 10 minutes antigen retrieval in citrate buffer. After washing in PBS, sections were blocked with 5% normal donkey serum, 1% bovine serum albumin and 0.2% Triton X-100 in PBS at room temperature (RT) for 1 hour. Sections were incubated overnight at 4 °C with the following antibodies: rabbit anti-Lef1(1:100, Cell signaling 2286S), goat anti-Pcadherin (1:400, R&D Systems, AF761), mouse anti-β-catenin (1:100, BD Sciences, 610153). After washing in 0.1% Triton-X 100, sections were incubated with the following secondary antibodies: Alexa Fluor 488-donkey anti-goat (1:200, Life technologies, A11055), Alexa Fluor 488-donkey anti-goat (1:200, Life technologies, A11057), Alexa Fluor 568-donkey anti-rabbit (1:200, Life technologies, A10042) for 1 hour at RT and DAPI mounted (Vector laboratories, H-1200). For β-catenin immunohistochemistry, biotinylated secondary antibody from M.O.M Kit (Vector Labs, BMK-2202) was used followed by detection using the ABC kit (Vector Labs, Pk-6100) and developed using DAB substrate kit (Vector Labs, SK-4100) according to the manufacturer’s instructions. Cryo-embedded blocks were cut at 10 μm thick sections and post-fixed with 4% paraformaldehyde (PFA) for 10 min at room temperature, followed by staining as described above. Combined images of DAPI, Alexa 488, Alexa 568 and Alexa 647 were taken using a Zeiss Axio Observer.

β-galactosidase detection
Whole embryos were fixed in 4% PFA for 45 minutes at 4 °C and then washed in PBS for 1 hour. Embryos were placed in 25 μl/ml of X-gal (VWR, 0428-100 mg) prepared in a staining solution of 5mM Potassium ferricyanide and 5mM Potassium ferrocyanide in 0.1% deoxycholate, 0.2% NP40 and 2 mM MgCl2 in PBS. Embryos were incubated in X-gal staining solution overnight at RT and washed with PBS before FFPE. Sections were subjected to fast Red (Sigma, N3020) nuclear counterstain.

Whole mount immunofluorescence
Dorsolateral skin from embryos were microdissected and placed on nucleopore filters (VWR, WHA 800281) and fixed overnight with 4%PFA at 4° C. Skin explants were blocked for 8 hours in 5% normal donkey serum, 1% bovine serum albumin and 0.5% Triton-X100 at RT. Explants were incubated overnight at 4°C in the following p rimary antibodies: rabbit anti-Sox2 (1:200, Abcam ab97959), rabbit anti-Sox9 (1:500, Millipore, ab5535), rabbit anti-Cytokeratin17 (1:1000, Abcam ab53707), goat anti-Pcadherin (1:400, R&D Systems, AF761), rabbit anti-Ccna2 (1:500, Abcam ab181591), chicken anti-GFP (1:500, Abcam ab13970). Explants were then washed in 0.2% Tween20/PBS for 6 hours on a rotator and then incubated with the following secondary antibodies: Alexa Fluor 488-donkey anti-goat (1:200, Life technologies, A11055), Alexa Fluor 568-donkey anti-rabbit (1:200, Life technologies, A10042), Alexa Fluor 488-donkey anti-rabbit (1:200, Life technologies, A21206) overnight at RT. Washes were carried out for 6 hours with 0.2% Tween20/PBS and explants were mounted on slides with Vectashield DAPI. Explants were imaged in 3 dimensions using the LaVision TriM Scope II (LaVision Biotec) microscope.

In-situ hybridization
Using the RNAscope 2.0 HD Detection-RED kit (ACDBio, 322350), single-molecule fluorescence in situ hybridization was performed according to the manufacturer’s protocol. Briefly, sections were deparaffinized, permeabilized with hydrogen peroxide followed by antigen retrieval and protease treatment before incubation with ISH probes for 2 hours. After probe hybridization, RNA signal was amplified and detected using the Amp 1-6 reagents provided within the kit for chromogenic development. Nuclear counter-stain was done using hematoxylin followed by 1.2% ammonia. Slides were then mounted with BioCare EcoMount (Biocare Medical, 320409). The following RNAscope probes were used: Mm-Bmp4 (401301), Mm-Ptch1 (402811), Mm-Lox (425311), Mm-Axin2 (400331), Mm-Dkk1 (402521) and Mm-Col1a1 (319371), Mm-Edar (423011), Mm-Cdkn1a (408551) and Mm-Lef1 (441861).
For dual staining following in situ hybridization, co-labelling with Sox2GFP reporter or any immunofluorescent stain was carried out according to the manufacturer’s protocol. Following ISH chromogen development as described above, sections were washed with PBS-T (0.1% Triton-X100), blocked with 5% normal donkey serum, 1% bovine serum albumin and 0.5% Triton-X100 in PBS at RT for 1 hour. Incubation with primary and secondary antibodies was carried out as described above for immunofluorescence stains and mounted using Vectashield DAPI mounting media.

Embryonic skin explant cultures
Embryonic skin explants harvested from E13.5 PDGFRαH2BGFP were cultured in DMEM/F12 (Gibco, 11039-021) supplemented with 10%FBS and penicillin-streptomycin. Skin explants were treated with Mitomycin C (MMC) (Sigma, 10107409001) in PBS for 3 hours at 1 μg/ml and 2.5 μg/ml. Explants were grown for additional 40 hours without MMC in culture; treated with 10μM EdU for 1 hour before fixing with 4% PFA for analysis.

Microscopy
IF stained paraffin and cryo-sections were imaged using the Zeiss Axio Observer Z1 equipped with a Plan APOCHROMAT 40x objective lens. Bright field images were taken using the Olympus BX61 microscope. Whole mount explants were imaged in 3 dimensions using the LaVision TriM Scope II (LaVision Biotec) microscope equipped with a Chameleon Vision II (Coherent) two-photon laser (940 nm for imaging of GFP, 880 nm for whole-mounts) to acquire z-stack images ranging from 50-120 μm (2 μm serial optical sections) using a 20X water immersion lens (numerical aperture 1.0; Olympus), scanned with a field of view of 0.3–0.5 mm2 at 800 Hz.

Single-cell dissociation
Embryonic dorsolateral/flank skin was micro-dissected and pooled (3-4 embryos per sample) and dissociated into a single-cell suspension using 0.25% trypsin (Gibco, Life Technologies) for 30 minutes at 37° C. Single-cell suspensions we re then stained with DAPI (Fisher Scientific, NBP2-31156) just prior to fluorescence-activated cell sorting.

Fluorescence-activated cell sorting
DAPI-excluded live skin cells were sorted on a BD Facs Aria II (Biosciences) sorter with a 100 μm nozzle. Cells were sorted in bulk for 10X Genomics platform.

Single-cell RNA sequencing and library preparation
Chromium Single cell 3’ Library and Gel Bead Kit v2 Chromium Single Cell 3′ Library & Gel Bead Kit v2 (PN- 120237), Chromium Single Cell 3′ Chip kit v2 (PN-120236) and Chromium i7 Multiplex Kit (PN-120262) were used according to the manufacturer’s instructions in the Chromium Single Cell 3′ Reagents Kits V2 User Guide. After cDNA libraries were created they were subjected to HiSeq 2500 (Illumina) sequencing.

10X scRNA-seq data pipeline to matrix
The transcriptomes of 7,067 and 6,096 live single cells from E13.5 skin replicates (n=2) and 6,394 and 6,152 single cells from E14.5 skin replicates (n=2) were sequenced. Raw 10X sequencing data was processed into a matrix employing the standard 10X CellRanger pipeline. Briefly, base call files were fastq format which were aligned to the mm10 reference genome followed by nUMI and barcode counting, constructing the nUMI count matrices. nUMI matrices were filtered, centered and normalized using Seurat (Butler et al., 2018; Macosko et al., 2015). Briefly, for each original run condition (i.e. E13.5 wildtype), cells were filtered to have > 1,000 genes, and nUMI > 2,500, but < 50,000. Data was then log scaled, centered and normalized to nUMI. A median of 10,882 and 8,469 nUMIs per cell for E13.5 replicates and 11,064 and 7,980 nUMIs per cell for E14.5 replicates were identified.

Principle Component Calculation
Principle components (PCs) were calculated using Seurat’s RunPCA function. Genes with dispersion > 0.8 were used for calculating PCs.

Cell type specification
For each original run, tSNE dimension reduction was performed on the normalized, centered, scaled nUMI count matrices employing the first 10 PCs. We then performed unsupervised clustering using the Seurat SNN clustering package, employing a resolution of 0.1. Clusters that were positive to Col1a1 were defined as dermal, while clusters positive for Krt10 or Krt14 were defined as keratinocytes. Any clusters co-positive for Col1a1 and Krt10 or Krt14 were considered doublets and removed from analysis.

Data merge and normalization
For construction of the dermal diffusion map of wildtype cells, we merged the wildtype E13.5 and wildtype E14.5 dermal cell matrices. We regressed out batch and nUMI using Seurat for data normalization, effectively normalizing by z-score. For the dermal diffusion map of all dermal cells, E13.5 wildtype, E13.5 mutant, E14.5 wildtype and E14.5 mutant datasets were merged and normalized. For comparisons between wildtype and mutant conditions in dermal cells and keratinocytes using t-SNE dimension reduction, wildtype and mutant data was merged and normalized (nUMI and batch) for E13.5 and E14.5 separately.

Diffusion Map Construction
To generate these maps, we computed PCA of the data and retained the first 40 principal components. Let {xi} for i=1,…,N be the 40-dimensional vectors representing the cells. We constructed an N × N sparse affinity matrix A by letting when xj is among the 50-nearest neighbors of xi, and σi is an adaptive bandwidth which is set to the distance from xi to its 50th neighbor. Next, we symmetrized the matrix, and normalized each row to sum to one, resulting in a Markov matrix. The eigenvectors of this sparse matrix corresponding to the largest eigenvalues are then computed using the ARPACK library interface to R implemented by the R package igraph (Csardi and Nepusz, 2006) and then visualized. We validate our results with a fixed bandwidth kernel, set to the mean distance of xi to its 50th neighbor (Figure S1D).

Gene Correlations of Differentiation Trajectory
Employing the diffusion map generated with WT E13.5 and E14.5 dermal cells, we calculated correlation statistics for the eigenvector 2 values and normalized gene expression for each gene (Table S1). The genes with the 100 highest positive correlations were submitted to DAVID for KEGG pathway analysis (Huang da et al., 2009).

Pseudo-ordering of cells
Diffusion maps were subset for cells lying along the DC differentiation trajectory, to allow us to find signals specific to the DC differentiation process and not associated with differentiation of other lineages. We removed ancillary trajectories; for wildtype analysis we removed cells with eigenvector 3 values > 0.01 and < −0.01 (Table S2). For the merged analysis we removed cells that had eigenvector 2 > 0.01 or < −0.05 and eigenvector 3 values > 0.04. To generate pseudo-order, cells were ordered by the DC differentiation vector, which we define as the eigenvector that corresponds with DC differentiation, as determined by gaining Sox2 and other putative DC markers. For the wildtype-only analysis the DC differentiation vector was the vector corresponding with the second largest eigenvalue, and for the merged analysis, this was the vector corresponding with the third largest eigenvalue. We employed a generalized additive model to smooth the data (Hastie, 1990), employing a Poisson distribution for gene expression and binomial distribution for cell cycle phase, using ggplot2 geom_smooth function.

Cell cycle estimation of scRNA-seq data
We employed Seurat’s cell cycle scoring using the following genes (Nestorowa et al., 2016). Briefly, averaged relative expression of these genes were used to calculate G2/M and S scores, which are used for binning cells into G2/M, S and G1/G0 bins.
G2/M genes: Hmgb2, Cdk1, Nusap1, Ube2c, Birc5, Tpx2, Top2a, Ndc80, Cks2, Nuf2, Cks1b, Mki67, Tmpo, Cenpf, Tacc3, Fam64a, Smc4, Ccnb2, Ckap2l, Ckap2, Aurkb, Bub1, Kif11, Anp32e, Tubb4b, Gtse1, Kif20b, Hjurp, Hjurp, Cdca3, Hn1, Cdc20, Ttk, Cdc25c, Kif2c, Rangap1, Ncapd2, Dlgap5, Cdca2, Cdca8, Ect2, Kif23, Hmmr, Aurka, Psrc1, Anln, Lbr, Ckap5, Cenpe, Ctcf, Nek2, G2e3, Gas2l3, Cbx5, and Cenpa,
S genes: Mcm5, Pcna, Tyms, Fen1, Mcm2, Mcm4, Rrm1, Ung, Gins2, Mcm6, Cdca7, Dtl, Prim1, Uhrf1, Mlf1ip, Hells, Rfc2, Rpa2, Nasp, Rad51ap1, Gmnn, Wdr76, Slbp, Ccne2, Ubr7, Pold3, Msh2, Atad2, Rad51, Rrm2, Cdc45, Cdc6, Exo1, Tipin, Dscc1, Blm, Casp8ap2, Usp1, Clspn, Pola1, Chaf1b, Brip1, and E2f8.

Quantification and Statistical Analysis

Image analysis
Raw image stacks acquired from whole-mount explants were imported into Fiji or IMARIS (BitPlane Scientific Software) for analysis. ggplot2, igraph and cowplot R libraries were used for graphical representation of the scRNA-seq data. For generation of one-dimensional heatmaps, scripts were adapted from https://github.com/KlugerLab/t-SNE-Heatmaps (Linderman et al., 2017). Diffusion maps were generated as described above.

Quantification
For quantification of hair follicle density, 3-dimensional whole-mount tiled mosaics of skin explants with a 1000 × 1000 μm field of view and a z depth of 60-100 μm were used for n=5 embryos. For volumetric quantification of DC cell number, HF epithelial cell number, percentage of Sox9+, EdU+ cells and mGFP+ cells, 3-dimensional whole mount mosaics stained with Sox2, cytokeratin17, P-cadherin, Sox9, anti-GFP and EdU were used to manually count positive cells, using ImageJ (Fiji) software. For volumetric quantification of cell cycle phases, 3-dimensional whole mount mosaics stained with markers such as Ki67, CyclinA2 or G0/G1 arrested cells (based on Fucci2 in vivo reporter) were done manually by counting nuclear positive cells using ImageJ (Fiji) software. For volumetric estimation of HF epithelium volume and DC volume, IMARIS (BitPlane Scientific software v9.1.2) was used. For quantification of ISH, sections with 4-5 dots per cell were considered positive (according to the RNAScope manufacturer’s instructions) and sections from a total of n=4 different embryos were examined. For quantification of sections stained for Lef1, Sox2, EdU, each cell marked by positive nuclear stain was considered positive and sections from n=5 different embryos were examined. For quantification of upper dermal cells, a distance of 30-40μm below the epidermis was considered as upper dermis. For quantification of cells in the peri-DC region an area in close vicinity to the dermal condensate was considered and cells that lie within that region were considered positive.

Statistical analysis
All statistical values are expressed as mean ± SD. An unpaired Student’s t-test was used to analyze data sets with two groups and *P <0.05 to, **P<0.01, ***p< 0.001 and ****P<0.00001 indicated a significant difference. When comparing more than two groups, P values were determined by one-way ANOVA with Tukey’s HSD test performed as the post hoc analysis. Statistical calculations were performed using Prism software package (GraphPad). For categorical data, a chi-square test was employed to estimate a P value. Permutation tests were used for estimating P values as indicated in the text.

Data and Software Availability
The data that support the findings of this study are available from the lead contact and on the NCBI GEOarchive (accession ID: GSE122043).

Supplementary Material
13Table S1: Related to Figure 1. Pearson correlation of genes along DC trajectory.

4Table S2: Related to STAR methods (pseudo-ordering of cells). Gene list of other dermal trajectories.

5Table S3: Related to Figure 4. Gene list of stalled point population along the DC trajectory.