Dermal Condensate Niche Fate Specification Occurs Prior to Formation and Is Placode Progenitor Dependent.

INTRODUCTION
Many adult tissues are maintained during homeostasis or repaired after injury by resident stem cells as a resource of reserve cells with multilineage potential. Stem cell functions are controlled both intrinsically and by a large variety of long-range neuronal and hormonal inputs and short range signals from the microenvironment or niche (Ge and Fuchs, 2018; Heitman et al., 2018; Rezza et al., 2014). In many cases, signals are derived from specialized niche cells that control stem cell quiescence for long periods of time, regulate stem cell self-renewal, or guide their differentiation into one or more lineages of the respective tissue. Neighboring niche cell types and their signals have been identified in several tissues such as bone marrow (Wei and Frenette, 2018), intestine (Shoshkes-Carmel et al., 2018), brain (Falk and Götz, 2017), lung (Lee et al., 2017) and skin (Ge and Fuchs, 2018; Sennett and Rendl, 2012). But a detailed understanding of the molecular controls of niche fate specification during embryonic development remains limited. Equally unclear is the precise timing with which niche cells emerge relative to stem/progenitor cells and their mutual impact on each other.
In hair follicles, rapidly dividing progenitors in the proximal bulb region receive signals from the centrally-located dermal papilla (DP) niche, which controls hair follicle progenitor proliferation, migration and differentiation into several lineages of the hair shaft and its channel during continuous hair growth (Clavel et al., 2012; Sennett and Rendl, 2012). During the hair cycle after bulb degeneration and rest, the DP niche sends activating signals to the stem cells in the hair germ to regenerate a new hair bulb ( Ge and Fuchs, 2018). Many studies have investigated the developmental processes and morphogenetic events during embryonic hair follicle formation, including the emergence of stem cell and niche compartments (; Sennett and Rendl, 2012). Around embryonic day (E) 12.5, secreted epidermal Wnts activate broad dermal Wnt/β-catenin signaling (Chen et al., 2012). This is upstream of yet unidentified dermal signals that at ~E13.5 induce the fate specification of hair follicle progenitors in patterned pre-placodes (Figure 1A, hair follicle stage 0) (Paus et al., 1999) that can be molecularly identified by Dkk4 and Edar expression (Sennett and Rendl, 2012). Progenitor cell migration then forms the physically identifiable placode (Pc), (Figure 1A, stage 1) (Ahtiainen et al., 2014). Following this, Pc progenitors signal back to the dermis for formation of dermal condensates (DC) (Figure 1A, stage 1, DC1), specialized cell clusters derived from fibroblasts (Fb), which act as signaling niches for Pc progenitors to regulate continued hair follicle development (Sennett and Rendl, 2012; Sennett et al., 2015). Pc progenitors also give rise to suprabasal Sox9+ precursors of the future adult hair follicle stem cells in the bulge after completion of hair follicle morphogenesis (Ouspenskaia et al., 2016). Continued signal exchange between the DC and Pc leads to progenitor proliferation and downgrowth of the polarized hair germ (Figure 1A, stage 2, DC2), with the DC at the leading edge that, after an elongated hair peg stage, becomes engulfed by the progenitors to form the DP niche of the hair-shaft producing bulb (Grisanti et al., 2013; Sennett and Rendl, 2012). Several key signaling pathways have been identified for progenitor and niche formation (Sennett and Rendl, 2012). Wnt/β-catenin signaling is essential and most upstream for Pc (DasGupta and Fuchs, 1999; Huelsken et al., 2001) and condensate (Tsai et al., 2014) formation. Eda signaling is then required for maintaining Wnt signaling and Pc stabilization (Zhang et al., 2009). FGF20, another Wnt target, is a key Pc-derived signal required for DC niche formation (Huh et al., 2013) in a process involving cell aggregation (Biggs et al., 2018).
The morphological detection of the DC after Pc initiation (Chen et al., 2012), combined with discovery of DC markers (Chen et al., 2012; Sennett and Rendl, 2012; Huh et al., 2013) and the definition of a DC gene expression signature (Sennett et al., 2015) has led to the current model of niche fate specification in response to progenitor signal(s) that coincides with condensate formation. However, the identification of unclustered dermal cells with localized high Wnt activity beneath Pcs (Zhang et al., 2009), and low-level expression of DC genes in immunophenotypically non-DC dermal cells (Sennett et al., 2015) raise the possibility of DC fate acquisition prior to cluster formation. Contextualizing these early events in hair follicle formation is further complicated by the rapidity with which morphogenesis occurs; Wnt signaling activity appears virtually simultaneously in Pc and DC (DasGupta and Fuchs, 1999; Tsai et al., 2014; Zhang et al., 2009), rendering parsing out the precise order and timing difficult. Recent work also demonstrated direct mesenchymal self-organization ability in the absence of an epidermal pre-pattern (Glover et al., 2017) further calling into question the origin and timing of niche fate specification (Figure 1A): 1) is the niche fate specified before niche cluster formation and signaling function and 2) how is the timing of niche fate acquisition related to progenitor fate and signaling?
Here, by combining 3D and 4D time-lapse live imaging with fluorescence-activated cell sorting, bulk and single-cell transcriptomics, proliferation kinetics, and mutant mouse models, we identify fated unclustered precursors of the DC niche prior to signaling center formation. Based on a molecular time-lapse of dynamically changing niche gene signatures and ordering the cell fates in developmental pseudotime we define the cell fate trajectory from Fb through DC niche precursors that are fated to become the clustered niche. In this process, DC precursors (pre-DC) undergo a transitional molecular state that is consolidated towards the DC niche biological fate through waves of transcription factor and signaling molecule upregulation. Finally, with epidermal Wnt signaling ablation, which prevents Pc formation, and with deletion of Pc-derived FGF20, which precludes DC aggregation, we demonstrate that early niche fate and emergence of unclustered DC precursors depends on pre-existing Pc progenitors and their signals.

RESULTS

Identification of Unclustered Precursors of Dermal Condensates
Previous studies have demonstrated that cell aggregation is the cellular mechanism inciting DC formation (Biggs et al., 2018; Glover et al., 2017), but when the DC niche fate is acquired and how this timing relates to placode progenitor specification remains unknown. To identify potential DC-fated precursors before cluster formation and niche establishment, we followed up on our previous work where we isolated with fluorescence-activated cell sorting (FACS) the DC as GFP+ cells from embryonic Sox2GFP knock-in reporter skin (Figure S1A) (Sennett et al., 2015). Both Sox2 promoter-driven GFP (Figures 1B and S1A) and Sox2 mRNA (Figure S1B) are highly expressed in the established DC (Driskell et al., 2009; Clavel et al., 2012), where condensed cells are visible underneath flat and barely invaginated EDAR+ hair placodes of stage 1 hair follicles (Figure 1B, open arrow) (Paus et al., 1999; Sennett et al., 2015). GFP levels remain high after initial Pc downgrowth and polarization of stage 2 hair follicle germs (Figure 1B, arrow). We named these two DC stages with high Sox2 expression DC1 and DC2 accordingly. We also observed by flow cytometry a low-level GFP population that was a mix of Schwann and DC-like cells (Figure S1A, S1B and S1C) (Sennett et al., 2015), but the difficulty of detecting low Sox2GFP levels in sagittal skin sections and absence of know-how to separately sort DC-like cells precluded further resolving this population (Figure 1B).
We surmised that low Sox2GFP-expressing DC-like cells could represent a precursor state that ramps up Sox2 expression towards a DC fate. To detect GFP low cells with high sensitivity we performed multicolor immunofluorescence in whole mount skin of E15.0 Sox2GFP embryos. In the top view, the anti-GFP signal was strong in DC1 and DC2 (Figure 1C, arrows). Also, low GFP levels were now detectable in Schwann cells (Figure 1C, asterisks), highlighted by ITGA6, forming a GFP low/ITGA6+ network within the developing skin (Figure 1C, asterisks). Importantly, we also found few loosely arranged GFP low/ITGA6− non-Schwann cells in the dermis in a location fitting a hair follicle pattern (Figure 1C, arrowhead). Co-labeling with previously identified DC marker FOXD1 (Figure S1C) (Sennett et al., 2015) labeled DC1, DC2 and these GFP low cells, suggesting that they might be DC precursors (pre-DC) before DC niche formation. High magnification images suggested that potential pre-DC cells are still unclustered, unlike the clearly visible cell aggregations of DC1 stage hair follicles (Figure 1D).
To demonstrate that potential pre-DC cells are associated with hair follicles we next performed whole mount immunofluorescence for SOX2 and EDAR on E15.0 PdgfraH2BGFP reporter back skin. Pdgfra is highly expressed in Fb-type mesenchymal cells (Figure S1C) (Hamilton et al., 2003), which are marked by nuclear GFP. In optical slices of 3D confocal scans, SOX2+ pre-DC cells at the dermal level were underneath EDAR+ Pc at the epidermal level and appeared to be unclustered compared to aggregated DC1 and DC2 (Figures 1E and S1D). Overlap of SOX2 with nuclear H2BGFP demonstrated the Fb-type of potential pre-DC cells. 3D reconstructions further illustrated the unclustered state of pre-DC cells (Figure S1E, Supplemental Movies 1–3). Z-projections demonstrated the close proximity of pre-DC to the overlying Pc (Figure 1F). Finally, to confirm that pre-DC cells are truly unclustered we measured, in 3D, the distance of each SOX2 and H2BGFP double-labeled nucleus to its five nearest neighbors for pre-DC, DC1 and DC2 (Figure 1G). Cells in DC2 clusters were significantly closer to each other compared to all others. Pre-DC cells, however, were significantly further apart compared to those of DC1 and DC2, confirming their loose initial arrangement. As other H2BGFP+ Fb were interspersed between pre-DC cells (Figure 1E), it suggests that early fate acquisition does not have to be synchronized among direct neighbors. In summary, our data thus far identified the existence of uncondensed potential DC precursor cells suggesting that the DC fate is specified prior to DC formation and function.

DC Precursor Cells Are in a Transitional Transcriptional State Towards a DC Niche Fate
After identifying potential DC precursors in stage 0 hair follicles with SOX2 protein, GFP reporter, and the DC signature gene FOXD1, we next sought to isolate this population for exploring the molecular controls of DC fate acquisition and its niche functions. Furthermore, the simultaneous presence of pre-DC, DC1 and DC2 from stage 0, stage 1 and 2 hair follicles, respectively, in E15.0 back skin provided us with an excellent opportunity to compare multiple parallel developmental stages in isolated cells. To separately purify Sox2GFP high DC1 and DC2 cells we utilized the previously identified DC marker CXCR4 (Sennett et al., 2014) that is highly expressed in DC2, has lower levels in DC1, and is absent in pre-DC (Figure 2A and 2B). To separately isolate pre-DC and Schwann cells, which are both Sox2GFP low and CXCR4−, we labeled pre-DC for mesenchymal marker PDGFRA and Schwann cells for ITGA6 to isolate them as Sox2GFPlow/CXCR4−/PDGFRA+/ITGA6− and Sox2GFPlow/CXCR4−/PDGFRA−/ITGA6+ populations, respectively. Additionally, we sorted dermal Fb as PDGFRA+/Sox2GFP- cells and PDGFRA−/Sox2GFP- cells as a negative (Neg) population containing all other remaining skin cells (Figure 2B). To confirm Fb and DC enrichment we performed qRT-PCR of known marker genes (Figure 2C). As expected, Lum was highly enriched in Fb, while Cxcr4 was absent in pre-DC and detected at increasing levels in DC1 and DC2; Sox2 expression was low in both pre-DC and Schwann cells and high in DC1 and DC2. Importantly, Schwann cell marker Fabp7 showed strong enrichment in Schwann cells and not in pre-DC (Figure S2A), and all other Fb/DC markers were also absent (Figure 2C). Finally, Foxd1 was highly expressed at all three DC stages, and not in Fb or Schwann cells, confirming successful isolation of pre-DC cells in agreement with our 3D imaging analysis.
We then proceeded with population-based transcriptome analysis for gene signature discovery by bulk RNA-sequencing (RNA-seq), as previously described (Sennett et al., 2015). We first globally compared the transcriptomes by principal component analysis of the top variably expressed genes (Figure S2B) and hierarchical clustering (Figure S2C), confirming the similarities between DC1 and DC2 and establishing multiple co-regulated gene groups within all six populations. Interestingly, pre-DC was positioned between Fb and DC1 suggesting pre-DC dually share characteristics of both dermal Fb and clustered DC. We next calculated significantly differentially expressed genes (DEGs) in each population by ANOVA (false discovery rate<0.05) and compared the resulting total of 6,454 genes. We established gene signatures using significant DEGs and selecting for genes with FPKM>1 and an expression fold change>2, resulting in 1,890 signature genes for all six cell types. To focus on the potential pre-DC and their relationship to Fb and the DC niche, we analyzed gene expression relationships in a four-way Venn diagram analysis of all significantly increased genes in Fb, pre-DC, DC1, and DC2 compared to Sch and Neg populations (Figure 2D and Table S1). A total of 1,157 genes were enriched in the four mesenchymal populations, including 692 population-unique signature genes. However, we also found several interesting population overlaps. As the DP and its embryonic DC predecessor are Fb-derived mesenchymal cells (Driskell et al., 2013), it was revealing that the pre-DC overlap with Fb (28 genes) was larger than both the unique overlaps of Fb with DC1 (1 gene) and with DC2 (5 genes) (Figure 2D, blue), suggesting that pre-DC shares transcriptional similarity with Fbs, fitting of a Fb-to-DC trajectory. Indeed, exploring gene expression shared between Fb and pre-DC demonstrated a transition of high-to-low expression from Fb to DC with intermediate levels in pre-DC (Figure 2E). Several genes included extracellular matrix components such as Ccbe1 and Eln (Figure 2F). Conversely, differential expression analysis confirmed a high overlap of 119 genes between pre-DC, DC1, and DC2 (Figure 2D, red), which included 58 previously described DC signature genes (Sennett et al., 2015). Also here, pre-DC expressed many of these shared DC signature genes at intermediate levels between Fb and DC1/DC2, suggestive of its developmental precursor status to morphologically identifiable DC (Figures 2E and 2G).
We next complemented population-based bulk RNA-sequencing with single-cell transcriptome analyses to allow for exploration of minute cell-to-cell transcriptional changes that define cell heterogeneity and differentiation (Grün et al., 201 6). Importantly, it also enables a high-resolution definition of lineage relationships along a developmental trajectory (Herman et al., 2018). To corroborate the developmental link of pre-DC between Fb and DC1/DC2, we used SORT-Seq (Muraro et al., 2016) to index-sort DC cells with our established strategy and obtain single-cell transcriptomes that are linked to the isolated cell identities (Figure 2H). We sequenced 1,371 Fb, 149 pre-DC, 152 DC1, and 151 DC2. Of those, 1,076 Fb, 107 pre-DC, 118 DC1, and 82 DC2 were used for analysis based on a minimum requirement of 6,000 transcripts/cell and established quality controls to exclude poorly sequenced cells. Based on population frequencies observed by flow cytometry it would require ~15,000 total unenriched cells to obtain similar transcriptome information of our sorted pre-DC cells (150x enrichment). We then applied the Rare Cell Type Identification 3 (RaceID3) clustering algorithm that was previously used to identify rare intestinal cell types and lineage relationships between stem cells and differentiated cells (Grün et al., 2016; Herman et al., 2018). RaceID3 clustered gene expression similarities were then visualized using t-distributed stochastic neighbor embedding (tSNE). Superimposing the Fb and DC cell identities onto the tSNE plot revealed clear distinctions between all isolated populations confirming our cell isolation strategy (Figure 2I). At the same time, it highlighted similarities between all DC stages and placed pre-DC between Fbs and DC1/DC2, similarly suggesting a lineage trajectory from Fb to DC through an intermediate pre-DC fate transition. As expected, all sorted DC populations expressed established DC signature genes, including Fgf10, Sox2, Foxd1, and Tbx18 (Grisanti et al., 2013; Sennett et al., 2015) (Figure 2J), and sorted Fb expressed established markers Col1a1, Lrig1, Dcn, and Lum (Driskell et al., 2013) (Figure 2K). Also here, expression patterns revealed gradual downregulation of Fb genes and upregulation of DC genes (Figure 2L). Altogether, these bulk and single-cell transcriptome data from purified embryonic skin cells indicate that unclustered pre-DC precursors are at a transcriptionally transitional state from Fb towards a permanent DC niche fate.

DC Precursors Are Developmental Intermediates Between Fibroblasts and the Mature DC Niche
As the identified pre-DC cells express genes characteristic of both Fbs and DC, we hypothesized that they are precursors of the mature DC niche at an intermediate developmental state between Fb and DC. To explore this lineage trajectory, we utilized RaceID3 to identify distinct cell clusters with similar transcriptome as a basis for calculating lineage relationships between related cell types (Figure 3A). Clusters 6, 7 and 8 were primarily comprised of sorted pre-DC, DC1 and DC2 cells, respectively (Figures 2H, 3A and S3A), again validating our sorting strategy. We then applied the Stem Cell Identification (StemID2) algorithm (Grün et al., 2016) to order RaceID3-identified clusters based on degree of transcriptional similarity (Figure 3B) inferring differentiation trajectories from single-cell RNA-sequencing data (Grün et al., 2016). Amongst DC clusters (clusters 6, 7, 8), pre-DC-enriched cluster 6 was closest to the Fb clusters. Amongst sorted Fb, cluster 5 was the sole calculated link to pre-DC cluster 6, suggesting that pre-DC cells emerged from closely related Fb (Figure 3B). It also demonstrated that pre-DC are the link for transition into DC1/DC2. To independently confirm the developmental trajectory, we utilized the Scanpy toolkit for Louvain clustering and diffusion maps (Wolf et al., 2018). Louvain clustering with the same 1,383 cells used for RaceID3 clustering, showed again clear distinctions between all sorted populations, placing pre-DC between Fb and DC1 (Figures 3C and S3B). Importantly, applying diffusion pseudotime, a method for ordering cells along lineage trajectories (Haghverdi et al., 2016), clearly demonstrated the lineage path from nearest Fb (Louvain cluster 5) through pre-DC towards DC1 and DC2 (Figures 3D, S3C).
Having established by single-cell RNA-seq analyses the differentiation trajectory from Fb to DC2, we next re-examined our bulk RNA-sequencing for gene expression trends. For this we performed hierarchical clustering noting four major trends: genes with highest expression in Fb (down A and B), and genes with highest expression in pre-DC (up C), DC1 (up D) and DC2 (up E) (Figure 3E and Table S2). 2,687 genes were downregulated from Fbs through DC2 (Figure 3F). Closer examination of these genes revealed major downregulation of Fb transcription factors, adhesion/ECM molecules, and signaling factors (Figures 3G and 3H) suggesting their involvement in rapid and sustained loss of Fb function as the DC fate is acquired and throughout DC differentiation.
Finally, we utilized our single-cell transcriptome data to define the pseudotemporal ordering of Fb and DC populations along the differentiation trajectory and explore patterns of gene regulation during niche fate acquisition. For this, we devised a self-organizing map using the Fb clusters closest to DC and the DC-enriched clusters (Figure 3I). We excluded Dlk1-expressing Fb clusters 1–3 that likely represented lower dermal Fb unrelated to DC (Figure S3D) (Driskell et al., 2013), and included Fb clusters 4 and 5 as both were closest to pre-DC-enriched cluster 6 (Figure S3E). Pseudotime gene expression analysis revealed refined trends of genes with highest expression in Fb (down clusters 1, 2), pre-DC (up clusters 3, 4), DC1 (up cluster 6), and DC2 (up clusters 5, 7) (Figure 3I and Table S3). Amongst 2,687 genes downregulated from Fb through DC2 in population-based RNA-seq (Figure 3F) and the 1,342 genes identified in single-cell transcriptomes (Figure 3J), 951 genes were in common (35% of bulk-identified targets, 70% of single cell-identified targets) (Figure S3F). Also here, these included several transcription factors (Lambert et al., 2018), adhesion and ECM components, and signaling molecules (Ramilowski et al., 2015) necessary for Fb function (Figure 3K). Creb3l1 is a secreted factor essential for collagen metabolism (Keller et al., 2018), while Irx1 and Twist1 are transcription factors that control expression of matrix-stiffening components (García-Palmero et al., 2016), and Zeb1 is known for promoting expression of vimentin (Liu et al., 2008). Additionally, expression of numerous collagens such as Col1a1, Col1a2, and Col3a1and other ECM components such as Dcn and Bgn, and P4hb were decreased between Fb and DC populations. Finally, we identified signaling molecules that are uniquely expressed in Fb, such as Efna5, Cntn4, Fbln1, and Wnt4 that are all downregulated on the path to DC lineage differentiation. Altogether these data demonstrate that Fb quickly halt transcription of Fb-associated genes as they transition through the pre-DC intermediate into the definitive DC niche fate.

Loss of Proliferation is Concomitant with Acquisition of DC Fate
It has been well established that niche cells in mature dermal papillae are quiescent and rarely divide (Tobin et al., 2003). Also the aggregated DC, as the embryonic DP precursor (Grisanti et al., 2013), was recently shown to be non-proliferative (Biggs et al., 2018). Therefore, we next explored in the unclustered DC precursors the precise timing of cell-cycle exit with respect to niche fate acquisition. As predicted, many cell cycle genes involved in molecular control of cell proliferation were downregulated from Fb to DC2 (Figure 4A). Interestingly, gene set enrichment analysis revealed a significant correlation of proliferation-associated genes in Fb compared to pre-DC, suggesting that exit of proliferation could be as early as during pre-DC fate acquisition (Figure 4B). Indeed, in pseudotime-ordered cells rapid declining expression of pro-proliferative genes such as Rbbp7, Mki67, and Pcna coincided with sharp upregulation of cell-cycle inhibitors, such as Cdkn1a and Btg1 (Zhu et al., 2013). Importantly, both events occurred prior to emergence of pre-DC and remained sustained through DC2 (Figures 4C and 4D), suggesting that cell cycle exit is a feature of pre-DC fate acquisition.
To independently test when proliferation exit occurs relative to DC fate acquisition, we performed EdU proliferation assays to mark proliferating cells during S-phase of the cell cycle (Figures 4E and S4A). EdU injection 6 hours prior to E15.0 analysis showed that 33% of Fb were proliferating (Figure 4F). By contrast, < 3% of DC1 and DC2 were going through S-phase suggesting that the vast majority of maturing DC cells had exited the cell cycle. Notably, only 3% of pre-DC incorporated EdU suggesting that already during early fate acquisition precursors shut down proliferation. Analyzing 12- and 24-hour chases showed increased EdU incorporation in pre-DC and DC1 fated cells, likely reflecting their developmental history from proliferating Fb at the time of EdU injection. Interestingly, with the 12-hour EdU chase DC2 had significantly lower proliferating cells than in pre-DC and DC1, and as low as 6h levels. Given that DC1 develops into DC2 during the 12-hour chase, these results suggest that DC1 progression to DC2 does not require cell proliferation.
Finally, to confirm that exit of proliferation is already occurring during DC fate acquisition before pre-DC detection, we tracked with 4D time-lapse live imaging the cellular processes up to pre-DC fate specification (Figure 4G). We live imaged embryonic ex vivo back skins from Tbx18H2BGFP mice that express high nuclear GFP levels in the DC (Grisanti et al., 2013) (Figure S4B). Importantly, pre-DC cells also already expressed high nuclear GFP, while Fb exhibited low GFP levels (Figure 4G). Live-tracking cell divisions during a 6-hour period prior to pre-DC fate acquisition detected actively proliferating interfollicular Fb (Figure 4G and 4H). By contrast, no pre-DC cells were dividing during the same time span, indicating that proliferation exit is occurring during DC fate acquisition in the Fb-to-pre-DC transition.

Three Distinct Waves of Transcription Factor and Signaling Molecule Upregulation Occur Along the DC Niche Developmental Trajectory
Besides proliferation shutdown and gradual loss of Fb gene expression signature, gradual upregulation of the DC niche signature occurs during DC fate acquisition (Figures 2E, 2G, 2L, 3E and 3I). We next explored the precise timing of molecular fate acquisition and consolidation from Fb to pre-DC, through DC1 and DC2. We discovered three distinct waves of gene upregulation for pre-DC fate acquisition, DC1 clustering, and hair germ downgrowth stages at DC2, in both bulk (Figure 5A) and single-cell (Figure 5C and Table S3) transcriptional analyses. The first wave was evident in gene clusters whose transient expression peaked early at pre-DC and then dropped, which likely includes factors necessary for fate acquisition and other early functions. Two additional waves of co-regulated gene clusters had sustained expression through DC2, but differed only in when their expression peaked: mid genes were enriched in genes achieving maximal expression at DC1, and late genes reached maximal expression in DC2, at the end of the pseudotime. Of the three upregulated clusters, the mid cluster contained the most heterogeneity in expression dynamics, and included genes with gradually ramped up expression beginning in pre-DC that is sustained throughout, and genes with more delayed upregulation patterns.
All three waves contained many transcription factors that could represent potential master regulators for DC specification and development (Figure 5B–5G and Table S4). Transcription factors expressed in the early wave included known key fate regulators such as Twist2 (Dermo1) (Li et al., 1995), the Wnt signaling transcription factor Lef1, and the TGF-β effector Smad3 (Figures 5B, 5D, 5F and S5A). Twist2 knock-out animals have sparse hair, among many other abnormalities (Šošić et al., 2003). Similarly, Lef1 knock-out results in abrogation of whiskers and back skin hair follicles (van Genderen et al., 1994). We confirmed Twist2 mRNA expression in pre-DC by in-situ hybridization, which is later shut down at the DC1 stage (Figure 5D). Many more transcription factors were in the mid gene cluster (Figures 5E, 5F and S5B). Mid-expressing transcription factors include Trps1, Tshz1, Foxd1, Prdm1, and Sox18. Mutations in Trps1 are associated with Ambras Syndrome, a hypertrichosis disorder (Fantauzzo et al., 2008). Mutations in Sox18 have been shown to impact differentiation of adult dermal papillae (Villani et al., 2017). Finally, we resolved many transcription factors whose expression peaked in DC2, including Notch target genes Hey1 and Hes5, as well as Glis2, Foxp1, and Alx4 (Figures 5E, 5F and S5C). Foxp1 has noted functions in maintaining hair follicle stem cell quiescence (Leishman et al., 2013). In total, we discovered 11 early-wave, 16 mid-wave, and 24 late-wave transcription factors (Figure 5G) that followed these temporally-restricted upregulation patterns along pseudotime, suggesting that a coordinated array of molecular changes is required for pre-DC to mature through DC1 and DC2.
As the dermal papilla acts as a major signaling niche for epithelial progenitors during postnatal hair growth and for stem cells in cycling adult hair follicles, we finally explored the expression and upregulation of signaling molecules during DC fate specification and niche formation. Also here, we discovered a distinct pattern of signaling waves (Figures 5H–5K and Table S4). Signaling factors expressed in the short time window around pre-DC include the Wnt inhibitor Dkk1, Cyr61, and Fst (Figures 5I, 5J and S5D). Cyr61 is dynamically expressed in wound healing to reduce fibrosis (Jun and Lau, 2010). Fst knock-out is associated with abrogated hair follicle development (Nakamura et al., 2003). The mid-acting signaling molecules included well known regulators Fgf10, the TGF-β modulator Ltbp1, and Sema6a (Figures 5I, 5J and S5E); late-acting signaling molecules included key ligands such as Igfbp4, Inhba, and Rspo3 (Figure 5I, 5J and S5F) – many of which continue to be expressed in dermal papillae (Rezza et al., 2016). The global display, as cumulative percent of total, corroborates a resolution of 8 early, 11 mid, and 24 late upregulated signaling factors (Figure 5K). Overall these data suggest that DC development, from fate acquisition in a subset of Fb through aggregated DC1 and DC2, is accompanied by three distinct waves of transcription factor and signaling molecule expression that in concert are likely regulating the specification of DC fate and initiation of niche function.

DC Cell Fate Specification Requires Signals from Pre-existing Placodes
Our data suggest that the transition to pre-DC fate from Fb and to niche formation is coordinated by waves of transcription factor and signaling molecule upregulation. Resolving whether this is a pre-programmed cell-autonomous process or whether it relies on intercompartmental signals from pre-existing epithelial placode progenitors directly establishes the developmental hierarchy of progenitor and niche fate specification. Wnt/β-catenin signaling in the epidermis is essential for Pc formation and hair follicle morphogenesis, including the formation of the mature DC (Huelsken et al., 2001; Zhang et al., 2009).We next asked whether the early pre-DC fate, which developmentally precedes DC1 and DC2, is established in the absence of Pc progenitors. To this end, we blocked Pc formation by conditional Wnt/β-catenin signaling ablation in embryonic epidermis of K14-Cre;β-cateninfl/fl mice. In E14.5 control skin we confirmed the developmental stages of SOX2+ ITGA6− DCs underneath EDAR+ Pcs, including pre-DC (Figures 6A and S6A). Additional double immunofluorescence with FOXD1 confirmed pre-DC cell identity (Figure S6B). In β-catenin-ablated skin, as predicted, EDAR+ Pc and SOX2+ ITGA6− mature DCs were absent (Figures 6B, S6A). Intriguingly, in the absence of Pc we also failed to detect groups of unclustered pre-DC, by either SOX2+/ITGA6− (Figure 6B, S6A and S6C) or SOX2+/FOXD1+/ITGA6− (Figure S6B) immunofluorescence, indicating that fate acquisition of pre-DC cells was entirely abolished. These data demonstrate that DC cell fate specification from Fb requires the presence of placode progenitors that likely provide essential Pc-derived signals.

Placode Progenitor-Derived FGF20 Is Required for DC Cell Fate Specification
A placode-derived signal potentially important for specifying the DC precursor fate is FGF20, which is both a well-known downstream target of Wnt/β-catenin signaling and also required for formation of the mature DC (Huh et al., 2013). Recent studies have further demonstrated that FGF signaling is required for cell migration to form mature DC clusters, however it is unclear whether it regulates DC fate acquisition itself (Biggs et al., 2018; Glover et al., 2017). Having discovered unclustered DC precursors prior to formation of mature DC, we next asked whether Fgf20 is required for induction of the pre-DC fate. To this end, we generated Fgf20-ablated embryos by crossing viable and fertile Fgf20Cre-GFP/Cre-GFP knockout mice with Fgf20LacZ/+ mice (Huh et al., 2012, 2015). In E14.5 heterozygous Fgf20Cre-GFP/+ control skins, we observed normal hair follicle development with multiple DC stages underneath GFP+ Pc (Figures 6C and S6D). In Fgf20 knockout skins of Fgf20Cre-GFP/LacZ embryos aggregated mature DCs were absent underneath GFP+ Pc, as described previously (Figures 6D and S6D) (Huh et al., 2013). Intriguingly, here like in the Wnt/β-catenin signaling ablation above, SOX2+/ITGA6− or SOX2+/FOXD1+ pre-DC cells were also entirely absent (Figures 6D, S6D, S6E and S6F), suggesting that FGF20 is a key Pc-derived signal essential for DC fate specification, prior to cell aggregation.

DISCUSSION
The timing of Fb-to-DC niche transformation has been traditionally described by the emergence of a characteristic tightly clustered DC beneath a thickened epidermal placode (Millar, 2002). At the earliest identified stage of bona fide HF morphogenesis, distinct Pcs develop before morphological DC clusters form. Recently, two studies have demonstrated that the structured DC forms through directed migration of dermal Fb (Biggs et al., 2018; Glover et al., 2017), but the precise timing and specific molecular drivers of DC fate-acquisition in relation to the physical aggregation process and to the timing of Pc progenitor specification remained unexplored.
Here we demonstrate that de novo expression of the definitive (clustered) DC markers, Sox2 (Clavel et al., 2012; Driskell et al., 2009) and Foxd1 (Sennett et al., 2015), first appears before any appreciable DC clustering. Taking advantage of Sox2 and other DC molecular markers, both pre- and post-clustering, we were now able to define the molecular landscape of the earliest fate-acquired DC precursors and their most closely related Fb, further refining our previous molecular characterization of the established DC niche (Sennett et al., 2015). Reinforcing our novel sorting strategy for population-based RNA sequencing analyses of distinct unclustered DC precursors and clustered DC stages, unbiased lineage prediction from single-cell transcriptomes further reconstructed the DC differentiation path of the bulk approach. Importantly, it also provided a highly resolved account of the molecular transitional states and developmental trajectory from non-committed Fb to DC precursors through DC1 and DC2.
Complementing the model of DC niche formation through directed migration from ex vivo live imaging (Biggs et al., 2018; Glover et al., 2017), we found that while Fb surrounding the DC can have robust proliferative activity, the earliest fate-acquired, Sox2+ DC precursors have significantly diminished proliferation indicating that the mature DC is a product of continual Fb recruitment (Figure 7A). We interestingly observed an early steady decline in the expression of proliferation-related genes in Fb most closely-related to pre-DC in pseudotime. This preemptive suppression of proliferation occurs prior to upregulation of early DC-fate genes and is particularly intriguing since the establishment of asymmetrically positioned DC is required in maintaining cell-fate asymmetry in the epithelium giving rise to both hair follicle stem cells and progenitors during early down growth of developing follicles (Cetera et al., 2018). Therefore, by extension, it can be posited that tight regulation of DC size by proliferation suppression is vital to maintain its asymmetric positioning for proper morphogenesis and hair follicle stem cell establishment. The potential necessity of cell-cycle exit for acquisition of DC fate, the cell intrinsic and/or extrinsic mechanisms to achieve cell cycle arrest, and the precise molecular controls that prevent over-recruitment through migration remain to be studied in more detail.
From both our bulk and single-cell transcriptome analyses, we identified transcription factors enriched in pre-DC including Twist2, Lef1 and Smad3 indicating an immediate shift in the gene regulatory network during early niche fate acquisition, potentially acting as trailblazers for downstream upregulation of mid- and late-expressed transcription factors and genes for DC formation and function such as those encoding cell-adhesion and signaling proteins (Figure 7B). The reconstructed lineage order in pseudotime also revealed a subset of Fb closely related to pre-DC raising the possibility of a yet earlier precursor stage to the DC niche fate. While pinpointing their precise spatial location poses a considerable challenge because of the paucity of distinctly expressed marker genes, clustering analyses place these cells broadly among upper dermal Fb which harbor multipotent dermal Fb progenitors (Driskell et al., 2013).
Epithelial-mesenchymal interdependence during early HF morphogenesis has been shown by multiple compartment-specific disruptions affecting the development of its counterpart (Chen et al., 2012; Tsai et al., 2014; Zhang et al., 2009). By preventing Pc formation through epidermal β-catenin ablation (Huelsken et al., 2001), we determined from the absence of pre-DC, early DC fate acquisition is not an autonomous pre-program of the upper dermis, but rather requires local directive signals from spatially patterned Pc. Genetic ablation of Fgf20 has previously demonstrated its requirement for DC formation by hair follicle stage 1 through cell migration (Biggs et al., 2018; Huh et al., 2013). Here we provide evidence that FGF20 directly promotes the early DC precursor fate before DC cluster formation. Similarly, recent identification of an Fgf20-expressing olfactory epithelial stem cell niche revealed a Wnt-regulated FGF20 requirement for formation of the underlying mesenchymal condensations that form nasal turbinates (Yang et al., 2018). Taken together these data indicate that the progenitors in nascent hair follicle pre-Pc signal to underlying Fb at least in part through FGF20 to induce the early DC fate prior to cluster-forming migration, and unequivocally demonstrate that progenitor fate precedes establishment of its supportive niche.
We now propose an updated model of progenitor and niche specification during early hair follicle morphogenesis (Figure 7A and 7B): Wnt-dependent Pc progenitors signal to underlying dermal Fb via the FGF/FGFR signaling axis initiating a cascade of dynamic transcriptional waves that regulate DC niche fate specification. Through multifaceted signaling interactions between Pc, pre-DC and unspecified interfollicular dermal Fb, fate specified DC precursor cells centrally migrate to form the clustered DC (DC1). Continued pre-DC fate acquisition and migration incorporates additional cells to the DC cluster, while in the earliest established DC the new transcriptional regulatory network – set up by the first transcriptional wave – begins to upregulate expression of the molecular machinery crucial for signaling functions from the DC niche to the progenitors above.

STAR METHODS

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Mice
Mice were housed at the Center for Comparative Medicine and Surgery (CCMS) at Icahn School of Medicine at Mount Sinai (ISMMS) and all performed animal experiments were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) at ISMMS. Timed matings were set up to obtain embryos at the ages indicated in figure legends (embryonic day 14.0 to 15.0), sex not known. Embryo ages were determined by designating the morning of vaginal plug detection as E0.5. Sox2GFP (Ellis et al., 2004), Tbx18H2BGFP (Grisanti et al., 2013) and PdgfraH2BGFP (Hamilton et al., 2003) fluorescent reporter mice were generated as previously described. Fgf20Cre-GFP and Fgf20LacZ mice for Fgf20 knockout were generated as previously described (Huh et al., 2012; Huh et al., 2015). β-catenin epidermal-specific cKO were generated by mating β-cateninfl/fl (Brault et al., 2001) with K14-Cre (Dassule et al., 2000) mice.

EdU incorporation
For EdU-uptake proliferation assays, pregnant females mated with Sox2GFP males were injected intraperitoneally with a single dose of EdU/PBS (500 µg/g body weight) 6, 12 or 24 hours prior to sacrifice. EdU detection was performed with the Click-It EdU AlexaFluor 647 Imaging Kit (Life Technologies) according to manufacturer’s instructions.

Immunofluorescence staining
Tissue sections were embedded in OCT (Tissue Tek) and cut at a thickness of 10 µm with a cryostat (Leica). Whole-mounted embryo skins or sections from E15.0 embryos were fixed in 4% PFA and washed with PBS. Whole-mount samples were additionally permeabilized for 2 hr in 0.3% Triton X-100 in PBS. Tissues were blocked in PBS-Triton with BSA/NDS for 2 hour at room temperature. Primary antibody labelling against EDAR (goat, 1:100, R&D), SOX2 (Goat, 1:100, Stemgent), FOXD1 (Goat, 1:100, Santa Cruz), GFP (chicken, 1:800, Abcam), and CXCR4 (Rat, 1:100, BD Pharmagen) were carried out overnight at 4ºC. After PBS washes, samples were incubated with Rhodamine Red-X-, AlexaFluor 488-, AlexaFluor 555-, or AlexaFluor 647-conjugated donkey anti-goat, rabbit, chicken or rat secondary antibodies (Jackson Immunoresearch, Invitrogen). Samples were again PBS washed, counterstained with 4’6’-diamidino-2-phenilindole (DAPI) and mounted onto slides with antifade mounting medium.

Image acquisition and processing
Whole mount immunofluorescence and stained tissue sections were imaged using Leica SP5 DMI confocal and Leica DM5500 widefield microscopes, respectively, both equipped with Leica LASAF software. For 3D whole mount imaging, z-stacks of up to 60 planes with 1 µm vertical intervals were acquired with Leica 40X, 63X or 100X oil-immersion lenses. 3D stack images from confocal microscopy were processed and analyzed with ImageJ/FIJI (NIH) or Imaris 3/4D image visualization and analysis software (Bitplane).

Single-molecule mRNA fluorescence in situ hybridization
Single-molecule mRNA fluorescence in situ hybridization (FISH) was performed on whole mount skin from E15.0 embryos using RNAscope (Advanced Cell Diagnostics) according to the manufacturer’s instructions. Briefly, back skins were fixed in 4% PFA, dehydrated through a graded methanol (MeOH) series (50%, 70% and 100% MeOH/H2O), and were stored at −20°C overnight in 100% MeOH. The following day, skins were rehydrated and pretreated with hydrogen peroxide and protease III (2.5 HD Reagent Kit-RED), and hybridized with the Twist2 probe (ACD #489121) overnight at 40 C. The following day, skins were washed, postfixed in 4% PFA for 15 minutes at room temperature, followed by signal amplification by provided hybridizing horseradish peroxidase-labeled probes and fast red signal detection, as described in the manufacturer’s kit instructions. Following the RNAscope protocol, slides were processed for co-immunofluorescence and imaged as described above. The red chromogenic substrate was visualized as autofluorescence in the red fluorescent channel.

E15.0 cell isolation and sorting
To generate single-cell suspensions from E15.0 Sox2GFP dorsal skins, tissues were collected and processed as previously described (Sennett et al., 2015). Briefly, back skins were harvested by microdissection and were pooled for FACS for bulk RNA-sequencing while skin from a single embryo was taken for single cell sorting prior to single-cell RNA-sequencing. Skins were then digested in a dispase (Invitrogen)/collagenase (0.03%, Worthington) solution with 20 U/µl of DNase (Roche) at 4ºC overnight followed the next day by 30 minutes incubation at 37ºC. Cells were then filtered through 40µm-pore nylon cell strainers, followed by centrifugation at 350xg for 10 minutes. Resuspended cell pellets were stained with primary antibodies against EpCAM (BrV605-conjugated rat Ab, 1:100, Biolegend), CXCR4 (PE-conjugated rat Ab, 1:100, Biolegend), ITAG6 (BrV421-conjugated rat Ab, 1:100, Biolegend), PDGFRA (biotinylated rat Ab, 1:50, Affymetrix). This was followed by secondary labeling with Streptavidin-APC (1:200, Biolegend). Dead cells were identified by DAPI uptake. Cells were sorted on a FACSAria II fluorescence cytometry sorter (BD Biosciences) for both bulk and single cell RNA-sequencing. For single cell RNA-sequencing, live fibroblasts, pre-DC, DC1, and DC2 were index sorted into 384-well hard shell plates (Biorad) pre-loaded with 5 µl of vapor-lock (QIAGEN) containing a 100–200 nl mix of RT primers, dNTPs, and synthetic mRNA spike-ins. After sorting, plates were immediately quick spun and frozen to −80°C.

Real-Time qRT-PCR
Total RNA from the FACS-sorted cells was purified with the Absolutely RNA Nanoprep Kit (Agilent) and quantified with the NanoDrop spectrophotometer (Thermo Scientific). Reverse transcription was performed with the SuperScript III First-Strand Synthesis System using oligo(dT) primers (Invitrogen). Real-time qRT-PCR was performed with a LightCycler 480 (Roche) instrument using LightCycler DNA Master SYBR Green I reagents (Roche). Fold changes between samples were calculated based on the 2−∆∆Ct method and normalized to Gapdh.

cDNA library preparation for bulk RNA-sequencing
Total RNA from the FACS-sorted cells was purified with the Absolutely RNA Nanoprep Kit (Agilent). RNA concentration and quality were measured by the Agilent Bio-analyzer. Samples with RIN (RNA integrity number) scores of 9.6 or higher were further processed. 5 µl of RNA from each sample was used for reverse transcription, followed by amplification using the RNA Ovation RNA-Seq System V2 (NuGEN). Using the Ovation Ultralow DR Library System (NuGEN), cDNA libraries were generated from 100 ng amplified cDNA with a unique barcoded adaptor for each sample. Library concentration and quality were quantified by Qbit (Invitrogen) and the Agilent Bioanalyzer. Samples were then sequenced on the Illumina HiSeq 200 platform using a 50-nt single-read setting.

cDNA library preparation for single-cell RNA-sequencing
Single-cell RNA-sequencing was performed using the SORT-seq protocol (Muraro et al., 2016). Briefly, after sorting cells were lysed at 65°C for 5 minutes. RT and second strand mixes were dispensed by the Nanodrop II liquid handling platform for cDNA library preparation (GC Biotech). The aqueous phase was separated from the oil phase after pooling all cells from one plate into one library, followed by IVT transcription. For each 384 well plate, 384 primers (1 library of 384 cells) was used for mRNA reverse transcription, conversion to double-stranded cDNA, and in vitro transcription. Each primer consists of a 24 bp polyT stretch, a 4 bp random molecular barcode (UMI), a cell-specific 8 bp barcode, the 5’ Illumina TruSeq small RNA kit adapter, and a T7 promoter. Single-cell double-stranded cDNAs were pooled and in-vitro transcribed for linear amplification following the CEL-Seq 2 protocol. Illumina sequencing libraries were prepared using TruSeq small RNA primers (Illumina) and sequenced paired-end at a read length of 75 bp using the Illumina NextSeq.

4D time-lapse live imaging
Tbx18H2BGFP E14.0 embryonic dorsal skin was harvested by microdissection in ice-cold PBS, sandwiched between a 8 µm Nuclepore filter (Whatman) and the bottom of a 35 mm Lumox culture dish (Greiner) such that the epidermal side of the skin was in contact with the Lumox membrane, i.e. the bottom of the dish. The skin-membrane sandwich was stabilized with Matrigel (Corning). Explants were cultured in DMEM (phenol red free, Invitrogen) containing 10% fetal bovine serum, 1% HEPES, and 1% penicillin-streptomycin. After overnight culture, 3D Z-stack images of up to 20 planes with 2 µm vertical intervals were acquired at 10 or 25 minute intervals during a total of 16–20 hours using a Zeiss LSM880 equipped with a 20X Plan Apo 0.8 NA air objective or a Zeiss Axio Observer Z1 Yokogawa spinning disk equipped with a 20x Plan Apo 0.4 NA air objective. During imaging sessions, explants were maintained in a live-cell chamber at 37ºC with 5% CO2 and humidity control.

QUANTIFICATION AND STATISTICAL ANALYSIS

Quantification of DC clustering state
3D reconstructed images for nearest neighbor measurements of pre-DC, DC1, DC1 and interfollicular dermal fibroblasts were generated in Imaris. Automated segmentation of each cell was then performed by placing a point sphere at the geometric center of each nucleus from which each point’s discrete coordinates were used to quantify the absolute distances between nuclei. Multiple hypothesis testing of the data from the average distance between 5 nearest neighbors among cells from pre-DC, DC1, DC2 and interfollicular fibroblasts were performed using Anova: single factor.

Quantification of proliferation by EdU-uptake
Percent EdU+ among pre-DC, DC1, DC2 and interfollicular dermal cells was manually counted from acquired 3D confocal scans in ImageJ/FIJI. Sox2GFP+/ITGA6− cells were considered as DC cells whereas GFP−/ITGA6− cells were counted as interfollicular dermal cells. For each injection time point skins from 2 embryos were analyzed, from which 3 DCs/regions of each morphological stage were quantified (i.e. 3 follicles x 4 stages per embryo per time point). Statistical analysis was performed using one-tailed Student’s t-test.

Quantification of proliferation by cell division during time-lapse live imaging
To measure proliferation with time-lapse live imaging, 4D images acquired from explants were processed and analyzed with Imaris. 15,625um2 areas (350×350pixels) containing pre-DC or interfollicular dermis were obtained by cropping from larger 20X fields. Areas containing pre-DC were identified by locating DCs at the final time point and back tracing in time through the clustering process over a 6-hour period whereas the areas containing interfollicular dermis were defined as not containing any DC over the entire imaging session. Total cell numbers and mitotic events were quantified manually. Hypothesis testing was performed using Student’s t-test.

Quantification of number of DCs per field of view
SOX2+ pre-DC were quantified in whole mount stained back skins of β-cateninfl/fl wild type and K14-Cre;β-cateninfl/fl conditional knockout, as well as in Fgf20Cre-GFP/+ heterozygous and Fgf20Cre-GFP/LacZ knockout. For each genotype 2 embryos were analyzed. Pre-DC were counted in 4 field of views of 10x widefield live scans (601,267 μm2) in each embryonic skin. Statistical analysis was performed using one-tailed Student’s t-test.

Bulk RNA-sequencing data analysis
All raw RNA sequencing reads were mapped to the mouse genome (mm10) with TopHat v2.0.3 (Trapnell et al., 2009) coupled with the Bowtie2 (Langmead and Salzberg, 2012) and aligned with default parameters. Transcriptomes were assembled and fragments per kilobase per million reads (FPKM) for each gene were computed with Cufflinks v2.1.1 (Trapnell et al., 2010) with default parameters. Differentially expressed genes (DEGs) were identified using Cuffdiff (with default parameters except for the library normalization method was upper quartile normalization, where FPKMs were scaled via the ratio of the 75 quartile fragment counts to the average 75 quartile value across all libraries) and ANOVA with the Benjamini-Hochberg correction for multiple hypothesis testing with significance cut off FDR <0.05. The Fisher exact test was used for enrichment analysis with the same multiple hypotheses testing correction procedure and cut off. Principle component analyses (PCA) were performed for samples using BioJupies (Torre et al., 2018). Population signature genes were defined by DEGs with a FPKM≥1, and a fold enrichment≥2 compared to all other populations. Signature genes of multiple population overlaps featured in the Venn diagram analysis were defined similarly as individual-population signature genes by which each overlapping population must individually meet minimal expression level and fold-enrichment constraints.
Selection of genes used in developmental trajectory analyses were obtained by excluding genes with too low expression across all stages (FPKMany population < 1) and genes with similar expression levels across each developmental stage (FPKMmax < 2*FPKMmin). The remaining genes were log2 transformed and z-score standardized to obtain 0-mean and unit standard deviation. Gene clustering was performed with Morpheus (Broad Institute: https://software.broadinstitute.org/morpheus) using Spearman rank correlation with average linkage. Gene ontology and KEGG pathway enrichment analyses were carried out using Enrichr (Chen et al., 2013). Gene Set Enrichment Analysis (Subramanian et al., 2005) was performed using GenePattern (Reich et al., 2006). Curated lists were used to identify transcription factors (Lambert et al., 2018) and signaling molecules (Ramilowski et al., 2015).

Single-cell RNA-sequencing data analysis
Analysis of single-cell RNA-seq data was performed independently in an R environment (R 3.5.1, RStudio 1.1.453) and a Python environment (Python 3.6.5, Jupyter 1.0.0). Cells with more than 6000 transcripts/cell were included in the analysis. Additionally, mitochondrial genes, ERCC spike-ins, and non-expressed genes were excluded from downstream analysis. The remaining 1383 cells and 18483 genes were analyzed independently using the RaceID3 package (R-based) (Grün et al., 2016; Herman et al. , 2018) or the Scanpy package (Python-based) (Wolf et al., 2018).
For RaceID3 analyses, normalization was performed by simple rescaling, which is done by dividing the transcript counts per cell by the total number of transcripts, and multiplying this by the minimum total number of transcripts across cells. Principal components with overrepresentation of genes related to GO annotations “cell proliferation” and “cell cycle” were removed using the CCcorrect function. Evaluation of transcriptome similarities by 1-Pearson’s correlation coefficient followed by k-medoids clustering was used to cluster cells and t-distributed stochastic neighbor embedding (tSNE) (Van Der Maaten and Hinton, 2008) was used for visualization.
Using the Scanpy module, highly variable genes were extracted, data was log-normalized, and principal components were computed. Principal components showing significant overrepresentation of genes linked to the GO annotations “cell proliferation” and “cell cycle” in the top or bottom 1% quantile of loadings were removed in a fashion similar to the CCcorrect function in RaceID3. Cells were clustered using the Louvain algorithm (Blondel et al., 2008), and visualized using UMAP (McInnes and Healy, 2018).
Pseudotime analysis was performed on a subset of cells: clusters 4, 5, 6, 7, and 8 in RaceID3 analysis and clusters 1, 2, 3, and 5 in Scanpy analysis which comprise dermal condensate cells and fibroblast precursors of the latter. Cell order along pseudotime was inferred using StemID2 (Grün et al., 2016; Herman et al., 2018) or diffusi on pseudotime (Haghverdi et al., 2016). Curated lists were used to identify transcription factors (Lambert et al., 2018) and signaling molecules (Ramilowski et al., 2015).
For easier reproducibility, the complete code is available online (https://github.com/rendllab and https://github.com/kasperlab).
Cumulative percent was calculated as the running sum of FPKM values divided by the total FPKM sum. Cumulative percent plots for early, mid, and late expressed genes were generated by calculating the mean cumulative percent values for each position (cell) in the SOM order.

Supplementary Material
23456789