An Integrated Transcriptome Atlas of Embryonic Hair Follicle Progenitors, Their Niche, and the Developing Skin.

INTRODUCTION
Hair follicle (HF) regeneration in the adult hair cycle is driven by quiescent and activated stem cells (SC) in the bulge and hair germ (Greco et al., 2009; Hsu et al., 2011), and coordinated by signal exchange from the dermal papilla niche (Rezza et al., 2014; Sennett and Rendl, 2012). While these adult SC and niche cells have been isolated and characterized at a genome-wide level (Blanpain et al., 2004; Greco et al., 2009; Lien et al., 2011; Morris et al., 2004; Rendl et al., 2005; Tumbar et al., 2004), a comprehensive and integrated analysis of the unique molecular features that define their embryonic precursors in nascent HFs, within the context of the entire developing skin, has been lacking. Identifying the internal and external signaling mechanisms during HF morphogenesis is key to understanding the dynamic epithelial-mesenchymal interactions at work during complex tissue development; it also is the first instrumental step in advancing clinical attempts to regulate hair growth and generate fully functional skin grafts including HFs.
HF morphogenesis is initiated during embryonic skin development after secreted epidermal Wnts activate broad dermal Wnt signaling activity (Chen et al., 2012), which in turn - through unknown downstream signaling - leads to hair placode induction in the epidermis and dermal condensate formation below (Figure 1A) (Millar, 2002; Schneider et al., 2009; Sennett and Rendl, 2012). Genetic fate mapping studies established placode cells as the earliest progenitors of all epithelial HF cells including adult SCs in the bulge (Levy et al., 2005), and dermal condensate cells as the precursors of dermal papilla/dermal sheath niche cells within the mature follicle (Grisanti et al., 2013a). Continued signaling between these progenitor cells is required for downgrowth and formation of mature HFs, but the precise identity and order of all molecular messengers driving this complex process are not yet fully known (Ahn, 2015; Biggs and Mikkola, 2014).
Prior investigations of these events have been narrow in their exclusive focus on only placode or dermal condensate cells. While placode progenitors are known to generate a handful of signals important for some aspects of hair morphogenesis – for example Eda for maintaining placodes (Laurikkala et al., 2002; Zhang et al., 2009), Fgf20 for inciting condensate formation (Huh et al., 2013), and Shh for promoting hair downgrowth (Chiang et al., 1999; St-Jacques et al., 1998) – much less is known regarding the dermal response and contribution to this crucial signaling exchange. Wnt signaling in dermal condensates is important for the progression of HF formation (Tsai et al., 2014), and a number of additional factors are distinctly upregulated in condensates compared to non-specialized dermal fibroblasts in embryo skin, but as of present few have proven required for HF formation (Grisanti et al., 2013a, 2013b; Rezza et al., 2015; Sennett et al., 2014). Importantly, the skin is incredibly heterogeneous by E14.5, when placodes and condensates first start to appear, and signaling from multiple sources in the micro- and macroenvironment could be important for directing hair growth and patterning through distinct mechanisms.
To systematically investigate the cellular complexity of developing embryonic skin and gain comprehensive insights into the molecular identity of HF progenitors and niche cells compared to non-hair inducing keratinocytes and fibroblasts, we conducted refined cell isolations and genome-wide transcriptome analyses by RNA-sequencing. Using double-transgenic reporter mice and specific antibodies, we isolated six distinct cell types from embryonic E14.5 mouse back skin, including placode progenitors and dermal condensate niche cells, as well as lineage-related epidermal keratinocytes and dermal fibroblasts, melanocytes and Schwann cells, and a mixed population comprised of all remaining skin cells. Therefore, any gene expressed in E14.5 skin can be attributed to a specific cell type and/or compartment using our inclusive gene expression atlas. We composed a molecular snapshot of an entire tissue with unprecedented cellular resolution, and mapped feasible modes of communication between specific cell types within the skin as HF formation begins. We further defined specialized signature expression profiles for each isolated cell type, composed of genes with the potential to control cell fates and in turn specific functionalities. Together with this work, we share our data in an integrative, searchable web database that enables the discovery and localization of genes of interest for further investigation. Our hope is that this publically available resource prompts the inception of additional studies so that the underlying molecular mechanisms of HF formation and skin development, including progenitor/niche fate acquisition and maintenance, will be further elucidated.

RESULTS

Isolation of HF Placode Progenitors, Dermal Condensate Niche Cells, and other Distinct Cell Types from Embryonic Skin
The first cellular constituents of new hair follicles (HFs) are epithelial placode cells that give rise to activated matrix progenitors and future bulge stem cells (SCs) of downgrowing HFs, and dermal condensate cells that form the future dermal papilla and dermal sheath niche. To gain comprehensive insights into the molecular makeup of these specialized cells we devised an innovative multicolor labeling and cell sorting strategy to purify placode (Pc) progenitors and dermal condensate (DC) niche cells during the first wave of HF morphogenesis at embryonic day (E)14.5 (Figure 1A). By simultaneously co-isolating epidermal keratinocytes (Epi), dermal fibroblasts (Fb), melanocytes (Mc), Schwann cells (Sch) and a population that contains all remaining skin cells (Neg) including an enrichment of endothelial and smooth muscle cells, we sought to define the unique molecular features of the progenitors and niche, along with other distinct cell types within the entire embryonic skin (Figure 1A). To this end, we employed a combinatorial approach of double-transgenic reporter mice with immunofluorescence staining of single cell preparations from E14.5 total back skin, followed by fluorescence-activated cell sorting (FACS) (Figure 1B). We first crossed Sox2GFP mice that express green fluorescent protein under the endogenous Sox2 promoter (Ellis et al., 2004) with Lef1-RFP mice that were engineered to express red fluorescent protein under a human Lef1 promoter fragment (Rendl et al., 2005). Sox2 is expressed in the mature dermal papilla (DP) and in embryonic DC niche precursors (Biernaskie et al., 2009; Driskell et al., 2009; Rendl et al., 2005; Tsai et al., 2010), and Lef1-RFP was previously used to isolate DP cells (Greco et al., 2009; Rendl et al., 2005). To simultaneously label all epithelial cells including Pc, we performed double-immunofluorescence staining for the cell surface markers E-cadherin (ECAD) and P-cadherin (PCAD) (Figure 1C). PCAD immunofluorescence was previously used at E17.5 to obtain a mixed population of Pc, downgrowing hair germ and hair peg cells from different HF formation waves (Rhee et al., 2006). At the earliest stage of first wave HF formation at E14.5, we found highest PCAD expression specifically in Pc, while ECAD was expressed by all epidermal cells. For the dermal compartment, GFP was strongly expressed in DC cells of E14.5 Sox2GFP;Lef1-RFP embryos (Figure 1D), as previously described (Clavel et al., 2012; Tsai et al., 2014). RFP was present in the DC and also broadly labeled Fb throughout the dermis (Figure 1D). In addition, RFP was expressed within the epidermal compartment in Mc, identified by co-expression of the marker DCT (Figure 1D). RFP+ Mc also expressed the cell surface marker E-cadherin (ECAD), allowing their distinction from RFP+Ecad− DC and Fb cells (Figure S1C). Additional combinations of immunofluorescence staining on GFP/RFP skin for ECAD, PCAD, or DCT together with the basement membrane marker Integrin-β4 (ITGB4) confirmed correct labeling and identification of each cell type (Figure S1A–C). In summary, our identified configuration of reporter and marker expression provided the necessary framework to proceed with cell isolations.
To purify all cell types by FACS, we first microdissected total back skin from E14.5 Sox2GFP;Lef1-RFP embryos and released all skin cells by combined dispase/collagenase digestion. Trypsin digestion was not required to obtain single cells and its omission allowed preservation of cell surface antigen detection (not shown). We then performed ECAD (allophycocyanin, APC)/PCAD (Brilliantviolet, BrV) double-immunofluorescence staining on a single skin preparation, followed by 4′,6-diamidino-2-phenylindole staining (DAPI) to label dead cells. Finally, all DAPI− live cells were subjected to the following nested FACS scheme to purify a total of eight separate cell populations encompassing the entire E14.5 embryonic skin (Figure 1E): ECAD shifted all epidermal cells with a further distinction by PCAD to isolate (1) ECAD+PCAD+ Pc (1% of live cells). The ECAD+PCAD− gate included a mix of Epi and Mc cells, which were differentiated on the basis of RFP expression to purify (2) ECAD+RFP− Epi keratinocytes (6.5%) and (3) ECAD+RFP+ Mc cells (0.4%). ECAD− cells were >80% of total back skin and included a mix of dermal subpopulations distinguishable by varying levels of GFP and RFP. Consistent with what we observed in sections, ~0.9% of cells in the dermis highly expressed GFP and RFP, representing the (4) GFPHighRFP+ DC compartment. (5) GFP−RFP+ Fb that exclusively expressed RFP made up ~46% of cells in the dermis. What had not been easily apparent from sections alone was the existence of cells expressing low levels of GFP that were either (6) GFPLowRFP+ (2.3%) or (7) GFPLowRFP− (0.8%). In addition, with flow cytometry it became clear that 14% of live cells represented a (8) GFP−RFP− “negative” (Neg) population. To account for all cells of the entire E14.5 skin and obtain a complete molecular snapshot at the genome-wide level, we collected these cells and defined their identity through marker gene explorations as described below.
We next performed qRT-PCR to test the expression of known marker genes in the correct corresponding sample (Figure 1F, Figure S1D, E). As expected, placode genes Dkk4, Pcad, Shh, Wnt10b and Edar were highly expressed by isolated Pc compared to all other samples, while the expression of epidermal keratinocyte markers Ovol1, Evpl and Krt1 was highly enriched in Epi cells. Mc cells uniquely expressed Dct, Mitf and Mc1r. From the dermal compartment, the DC expressed high levels of Tbx18, Sox2, Bmp4, Trps1 and Nrp2, while Fb cells were strongly enriched for known fibroblast markers Col1a1, Lox and Lrig1. Expression of neural crest and Schwann cell genes including Gfra1, Fabp7, Mpz and Gap43 by GFPLowRFP− cells indicated that these were Schwann (Sch) cells in varying stages of maturation. qRT-PCR analysis revealed that GFPLowRFP+ cells expressed both DC and Sch genes but consistently to a lesser extent than each cleanly isolated cell type, indicating that this sample was a mix of both (DC/Sch). Finally, GFP−RFP− Neg cells included both endothelial and smooth muscle cells, as evidenced by the exclusive expression of known marker genes Pecam1 (CD31), Flt4, Flk1 and Nos3, and likely encompassed other stray dermal cells that we did not specifically investigate further. Taken together, our qRT-PCR marker analyses verified our sorting strategy and confirmed that we were accurately able to isolate embryonic HF progenitor and niche cells, along with four other specific cell types and two mixed cell populations inclusive of the entire embryonic skin.

Genome-wide Transcriptome Analysis of Independently Isolated Skin Cell Types
In order to comprehensively examine the unique gene expression profiles of the isolated HF progenitors, niche cells and other cell types at the genome-wide level, we turned to RNA-sequencing (Figure 2A). Using the sorting strategy outlined above we isolated all eight cell populations in two separate FACS sorts to obtain biological replicates. Using only 6ng of total RNA we first generated amplified cDNA (NuGEN Ovation RNA-seq System) that faithfully maintained relative gene expression patterns, as determined by additional qRT-PCR comparisons to conventional cDNA (Figure S2A). We then produced cDNA libraries with unique barcoded adapters for each of the 16 samples (NuGEN Ovation Ultralow Library System) that were subsequently analyzed by next-generation multiplexed sequencing (Illumina HiSeq 2000). This approach resulted in a high-quality output, with a >30 Mean Quality Score (Q score) and >90% Perfect Index Reads for all samples (Figure S2B). On average, 55 million total reads and 48 million aligned reads were produced per sample (Figure S2C).
Next, we proceeded by mapping, aligning and quantifying these reads to compute differentially expressed genes between all cell populations. By ANOVA, we identified a total of >7500 genes that were significantly differentially expressed (FDR<0.05) in at least one out of the eight cell populations. To gain a first glimpse into the global similarities and differences of gene expression patterns between all embryonic HF progenitor, niche and other skin cell types, we interrogated their relatedness by calculating coefficients of determination (r2). As expected all replicates displayed the highest correlation among all comparisons (r2>0.97) (Figure 2B). While Epi and Pc cells were closely related, the comparison of either sample to any other resulted in plot-wide low r2 values. These data reflected the highly unique gene expression patterns that are shared by Epi and Pc cells but distinct from all others. Disregarding the mixed samples, DC cells and Fb were the next closely correlated cell types, with the comparison between Sch and Mc almost as tight. At the other end of the spectrum, DC and Mc samples were most dissimilar.
Next we performed hierarchical clustering of the differentially expressed genes allowing further grouping of genes and populations according to their expression patterns (Figure 2C). This analysis also established that all replicates clustered together (Figure 2D), and grouped epithelial Epi and Pc cells apart from dermal and neural crest cell types. Closely related Fb and DC cells clustered on one side with the DC/Sch population, whereas neural crest-derived Mc and Sch cells were clustered separately, and Neg samples localized to their own branch of the dendrogram. Principal component analysis provided an additional way to calculate and visualize the relationship between the eight samples (Figure 2E). Again and as expected, replicates grouped together, with Pc and Epi localized closely in their own space. Neural crest-derived Sch and Mc samples associated near the bottom of the plot, while dermal Fb cells accompanied DC in a separate corner. The DC/Sch samples appropriately appeared between the two purely isolated cell types, while the Neg samples were located separately.
Finally, gene ontology analysis of the differentially expressed genes across all cell types revealed significant enrichment for genes related to signaling and migration (Figure 2F), highlighting the central roles of intercellular cross-talk and dynamic cell rearrangement in promoting skin and hair development.

Genome-wide Gene Expression Analysis Corroborates Previously Established Gene Markers
Several prior studies have examined the expression and functional role of specific genes in Pc and DC cells during HF morphogenesis at E14.5 (Ahn, 2015; Biggs and Mikkola, 2014; Sennett and Rendl, 2012). Global gene expression in embryonic Epi and Mc cells has been explored by microarray analysis and expression of selected genes was examined in more detail (Bazzi et al., 2007a; Colombo et al., 2012). Further, several genes mark Sch cells at precise developmental stages (Jessen and Mirsky, 2005; Woodhoo and Sommer, 2008) and define heterogeneity in embryonic Fb (Driskell et al., 2013). To further validate the findings from our RNA-seq analysis we examined the fold expression changes of known marker genes from the established literature by comparing their FPKMs in a target sample with those of its most closely related cell type (Figure 3). All examined marker genes were highly enriched in the correct cell type supporting the accuracy of our isolation strategy and subsequent expression analysis methods. Additionally, we verified the marker gene expression pattern for select genes by qRT-PCR. Interestingly, the fold enrichment observed by qRT-PCR was frequently much greater than by FPKM comparisons, indicating that even minor but significant differences in gene expression found by sequencing could identify uniquely localized signature genes.

Defining Distinct Molecular Signatures of Embryonic HF Progenitors and their Niche Alongside Other Components of Developing Skin
Because unique expression of genes in a given cell type predicts functional relevance, we next defined cell type-specific molecular signatures for each of the six cleanly isolated skin cell populations (Pc, Epi, Mc, DC, Fb, Sch). For this we compared the FPKMs of all genes that were significantly differentially expressed across all samples or between any sample pairs, and genes expressed with FPKM>1 and at two-fold greater levels in one sample compared to all others were labeled signature genes (Figure 4A, non-overlapping sections). For the Pc versus Epi and the DC or Sch versus mixed DC/Sch comparisons we applied a 1.5-fold cutoff given the highly overlapping nature of these gene expression profiles. In total we identified 1728 signature genes for all six cell types, ranging from 102 genes for Pc to 395 genes in DC cells (Figure 4A, Table S1). Additionally, we classified highly enriched genes specific to cells categorized by epithelial, mesenchymal or neural crest origin, compared to all other six cell populations (Figure 4A, intersections). 336 epithelial genes were shared by Pc and Epi cells, while 146 distinct dermal genes were expressed by DC and Fb cells. Finally, 67 genes overlapped between neural crest-derived Mc and Sch cells (Table S2).
To gain global mechanistic insights into how each cell type might function during hair and skin formation we next computed the enriched gene ontology terms (GO) of all gene signatures (Figure 4B, Table S3) using Enrichr (Chen et al., 2013). Genes relevant to Wnt receptor signaling were highly associated with both Pc and DC signatures, consistent with the known central role of this pathway in promoting HF development (Andl et al., 2002; Chen et al., 2012; Reddy et al., 2001; Tsai et al., 2014; Zhang et al., 2009). Genes that promote epithelial cell differentiation were well-represented in the Epi signature, in addition to genes related to keratinization and lipid biosynthetic processing. In the Fb signature, genes related to adhesion and ECM organization prevailed. GO analysis of highly specialized Mc and Sch cells revealed equally distinct gene categories such as pigmentation or axonogenesis/myelination, respectively. Genes related to migration, adhesion and signaling were highly represented within the DC signature, and notably several top hits were related to neural development and axon guidance. This observation was previously made for mature dermal papilla (DP) cells of postnatal follicles as well (Rendl et al., 2005).
Finally, we explored how embryonic progenitor and niche gene signatures were related to their mature postnatal counterparts. Several prior studies established gene signatures of cells within the SC niche during HF growth and/or the regenerative hair cycle. Lien et al. defined gene signatures for HF SCs at multiple points during the hair cycle using microarray technology, and established an overarching adult SC signature consisting of 672 genes (Lien et al., 2011). From 102 embryonic Pc signature genes a total of 13 genes (13%, p<0.001) overlapped with the adult HF SC signature (Figure 4C, Table S4). These common genes largely have been unexplored so far in either placode or adult HF SC function, and are discussed in more detail below. A separate study examined genes unique to DP niche cells, the progeny of embryonic DC cells, in addition to establishing signatures of melanocytes and HF progenitor subtypes in growing HFs of postnatal day (P)5 mouse back skin (Rendl et al., 2005). Out of 395 signature DC genes, 31 (8%, p<0.001) overlapped with the postnatal DP signature (Figure 4C). Intriguingly, many common genes were related to cell-cell signaling, and while several have been explored regarding condensate function during HF formation (Sox2, Enpp2, Cxcr4), many are newly identified and discussed below. Interestingly, many genes were shared between embryonic and postnatal melanocytes, with 75 of 269 (28%, p<0.001) embryonic Mc signature genes also specifically expressed by postnatal melanocytes. The overlap included well-established melanocyte markers, but also identified new common genes. In summary, with genome-wide sequencing we defined comprehensive molecular signatures of embryonic HF progenitors, the niche, and other key cell types within the embryonic skin.

Many New Signature Genes Typify Unique Molecular Characteristics of Pc and DC Cells in Embryonic HFs
In addition to successfully detecting established Pc and DC genes (Figure 3), our analyses revealed many enriched genes not previously examined in embryonic mouse skin (Figure 5A, Table S1). In Pc, many such signature genes were related to cell communication and specifically modulation of canonical signaling pathways. Pc ligand Pthlh binds to receptor Pth1r highly expressed in DC and Fb cells, is known to be important for the formation of mammary glands and teeth, and likely plays a role in mediating epithelial-mesenchymal interactions during skin/hair development as well (Hiremath et al., 2012; Philbrick et al., 1998). Pthlh specifically modulates mesenchymal BMP signaling, a pathway that needs precise regulation for HF formation to occur (Botchkarev et al., 1999; Jamora et al., 2003). Crim1, a factor with potential to be expressed extracellularly or secreted to affect Vegf, Bmp, Pdgf, and Tgfβ signaling, was found within the Pc signature and also overlapped with the published HF SC signature. Dkk4 is a Wnt signaling inhibitor known to be highly expressed in the Pc (Bazzi et al., 2007b), but interestingly co-receptor Kremen2 was found highest expressed in the Pc. Scube3 was also found in the Pc signature, and is a secreted factor that activates Tgfβ signaling by interacting with Tgfbr2, highly expressed by DC, Sch, and Fb cells. Finally, because Wnt signaling within the epithelium is necessary to drive HF formation (Chen et al., 2012), it was interesting to note the presence of lone frizzled receptor Fzd10 in the Pc gene signature.
Aside from canonical signaling factors, genes representing other signaling modalities were found enriched in the Pc signature. Gap junction genes Gjb2 (Cx26) and Gjb6 (Cx30) were both enriched, and are known to form heterodimeric channels to facilitate communication between linked cells, representing one way Pc progenitors might coordinate their activity as a unit. Another signature gene Syt7, a calcium sensing member of the synaptogramin family, facilitates signaling by promoting vesicle fusion and exocytosis of secreted factors. Prostaglandin signaling could play a role in coordinating skin formation, as synthase Ptges was highly expressed in the Pc and Ptge receptors were present in both Epi and Fb cells. Intracellular components important for transducing signals were uniquely expressed in the Pc as well, including Nedd9 which was found highly expressed in both Pc and adult HF SC signatures. Another example, Art4 (CD297), regulates NAD metabolism to promote signaling through multiple cascades. Finally, Ifitm3 is a receptor known to be involved in interferon signaling, and was also identified as a highly expressed gene within the Pc signature in E14.5 skin.
A number of transcription factors were found differentially expressed in the Pc compared to interfollicular Epi. Ascl4 is a bHLH transcription factor, and although not well studied, has been implicated in skin development because of its restricted, high expression in fetal skin (Jonsson et al., 2004). Another uniquely expressed transcription factor was Sox21, previously characterized in cells of the growing hair cuticle (Kiso et al., 2009). Expression of Sox21 at E14.5 was not described in this prior study, but its ablation resulted in a severe, progressive alopecia, despite what appeared to be normal hair formation through the first week of life.
Multiple signature genes in the DC were similarly involved in cell-cell signaling, which could be particularly relevant in coordinating events during the earliest stages of HF formation. Canonical Shh and Wnt signaling pathway inhibitors Hhip and Wif1, respectively, were noticeably enriched in the DC signature. Overlap with the previously published postnatal DP signature suggests an essential role in condensate function during hair development. Curiously, indicators/potentiators of active Shh and Wnt signaling including Ptch1/2, Rspo3 and Wisp1 were also uniquely expressed in the DC signature. Similarly, while Fgf signaling inhibitors Spry1 and Fgfrl1 were highly expressed in the DC signature, so was activating ligand Fgf10 (Suzuki et al., 2000). Previous reports have highlighted the role of the DP in fine-tuning cell-cell signaling during hair growth (Clavel et al., 2012); it seems likely that the DC has a similar function during hair development.
More obscure signaling factors expressed exclusively in DC cells in embryonic skin included Ndnf (also known as Epidermacan), a relatively unexplored secreted factor; Dclk1, a kinase implicated in neural cell migration; Grin2a, an NMDA receptor subunit important for synaptic transmission; and Syt6, another member of the synaptogramin family instrumental in vesicular exocytosis and signaling execution. Pcsk5 is also important for the mechanics of signaling, as an enzyme that activates proproteins via site-specific cleavage, and its substrates include peptide hormones and adhesion molecules that could influence signaling through DC–cell contacts.
In addition to factors related to signaling, numerous transcription factors were unique to the DC gene signature. Interestingly, multiple Sox transcription factors were found enriched in the DC. The expression and function of Sox2 is well-established in the DP (Clavel et al., 2012; Driskell et al., 2009), but this analysis additionally highlighted the localized expression of both Sox9 and Sox18 in E14.5 DC. As these Sox genes belong to separate Sox family subgroups they likely act non-redundantly to influence condensate cell activity and maturation. Foxd1, a forkhead box transcription factor with a known role in kidney development, was especially enriched in DC compared to Fb cells (Fetting et al., 2013). Kruppel-like zinc finger Glis2, extensively studied for its role in kidney development as well, was identified in the DC signature while related protein Glis1 was also highly expressed in these cells. Both factors can act as transcriptional activators or repressors to influence a multitude of other intracellular signaling cascades.
The few genes described here represent only a small fraction of the most interesting new signature genes of Pc and DC. qRT-PCR and immunofluorescence for select genes verified expression at the RNA and protein level (Figures 5B, 5C, S3). In total we confirmed enrichment of over 200 genes within Pc, Epi, Mc, DC, Fb and Sch signatures by qRT-PCR (total of 82 genes shown in Figures 1, 3, 5, S1, S3). Our sampling approach to verify the sequencing data substantiated virtually all inspected signature genes, and the few factors highlighted in Figure 5 provide only a small glimpse of the hundreds of new signature genes defined for multiple cell types in embryonic skin by our inclusive and in-depth experiments (Table S1).

Distinct Expression Patterns of Genes Related to Canonical and Other Defined Signaling Pathways, Including Axon Guidance Signaling, in Embryonic Skin
Many signature factors in multiple cell types of embryonic skin were related to canonical signaling pathways Wnt, Shh, TGFβ and Notch, among others (Figure 4B, 5A), and have been studied in the context of Pc and DC cross-talk in great detail (Lee and Tumbar, 2012; Wang and Wu, 2014). To understand how canonical signaling factors might influence skin development on a broader scale, we constructed heat maps to illuminate patterns of relevant gene expression for specific pathways between all six isolated cell types (Figure 6A and S4A). Genes were mined from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Kanehisa and Goto, 2000; Kanehisa et al., 2014) and included in the maps as long as they were expressed (FPKM>1) by one of the six cell types included in the analysis.
This illustrative strategy made it easily apparent that in the case of Wnt signaling, most ligands were expressed from the epithelial compartment, with especially high expression in the Pc. The exceptions were non-canonical Wnt11 and Wnt5a/5b expressed strongly by DC and Fb cells. Multiple Wnt signaling inhibitors were also expressed highly in the dermal compartment, including Dkk1, Wif1 and all Sfrps. In a separate map Fgf20 expression, known to be important for DC formation (Huh et al., 2013), was correspondingly strong from Pc, while Fgfr1 expression was notably high in DC. Ligands, receptors, and downstream factors for Bmp signaling were expressed in both epithelial and dermal compartments (Botchkarev and Sharov, 2004). Bmp signaling inhibitors, including Grem1, Grem2, Chrd and Nog however, were strongly localized to Fb cells. Analysis of Notch signaling factors revealed strong and specific expression of Notch binding partner Dll1 in DC cells. Although no single Notch was exclusively found in DC cells, downstream effector Rbpj and direct targets Hes5 and Hey1 were highly expressed. Finally, pathway analysis confirmed Shh ligand was specifically produced by Pc cells, as has been well-established (Chiang et al., 1999; St-Jacques et al., 1998), and that both Ptch and Smo receptors were strongly present in DC and Fb cells. Downstream effectors of Shh signaling including Gli1 and Gli2 were also expressed in DC cells. In summary, whole transcriptome analyses of embryonic HF progenitors, the niche, and other key skin cell populations enabled construction of detailed maps for multiple canonical signaling pathways and a comprehensive overview of each pathway’s factors within embryonic skin.
We next used Enrichr (Chen et al., 2013) for KEGG pathway analyses to ascertain overrepresented signaling pathways in signature gene lists. Factors related to both Shh and Wnt signaling were prevalent in the Pc signature (Figure 6B, Table S5). General factors related to cell communication were found in both Pc and Epi signatures, while genes involved in melanogenesis were most significantly overrepresented in the Mc gene list. This analysis also revealed genes important for Shh, TGFb, and P53 signaling within the DC signature. As described above, the roles of several canonical signaling pathways are well established in driving HF formation, although the importance of P53 signaling has yet to be addressed. In addition, KEGG analysis found that factors related to axon guidance were significantly overrepresented within the DC signature. Correspondingly, genes related to ECM interactions were highly and significantly present in the Fb signature, while ErbB signaling that is known to be important for Schwann cell maturation was highlighted by KEGG analysis of the Sch signature.
The significant enrichment of genes related to axon guidance within the DC was particularly intriguing. On closer inspection, our analyses uncovered that both ligands and receptors important for signaling through Chemokine, Semaphorin, Neurotrophin, Neuropilin, Netrin, Ephrin, Slit/Robo and Shh pathways were part of the DC signature (Figure 6C). Importantly, many corresponding ligands and/or receptors were highly expressed in neighboring Pc, Sch or Fb cells. Overall, the prevalence of strongly expressed genes related to axon guidance signaling in multiple cell types within embryonic skin was apparent through GO analysis, KEGG pathway analysis, and inspection of signature gene lists, suggesting a functional role in promoting hair and skin development.

DISCUSSION
The notion that epithelial-mesenchymal crosstalk drives HF development, growth and cycling has prompted investigations to establish gene signatures for adult SCs and the niche (Blanpain et al., 2004; Driskell et al., 2009; Greco et al., 2009; Lien et al., 2011; Morris et al., 2004; Rendl et al., 2005; Tumbar et al., 2004). However, only limited insights into the molecular make-up of their embryonic precursors have been generated by studies that isolated and profiled total epidermis (Bazzi et al., 2007a), or a compound population of progenitors from different HF formation waves (Rhee et al., 2006), or condensate niche cells without collecting any other population for comparison (Tsai et al., 2014). Further, these prior studies were performed using microarrays that, while cutting-edge at the time, have been surpassed by unbiased deep-sequencing technology. To date, a comprehensive and refined analysis of embryonic precursors to both adult HF SCs and niche has been lacking, hampering the systematic discovery of molecular players coordinating HF morphogenesis. Here we solve this problem with a high-resolution, in-depth exploration of gene expression in the earliest placode (Pc) progenitor and dermal condensate (DC) niche constituents of nascent HFs within the broader context of the entire embryonic skin. We co-isolated pure Pc and DC cells from the first HF development wave, in parallel with lineage-related keratinocytes (Epi) and fibroblasts (Fb), neural crest-derived melanocytes (Mc) and Schwann cells (Sch), and a population inclusive of all remaining cells within embryonic E14.5 skin (Neg). Rigorously comparing the genome-wide RNA-sequencing analysis of each isolated cell type or population captured a molecular snapshot of the whole tissue with unprecedented cellular resolution, and precipitated signature gene discovery. Analyzing these signatures to assess gene ontology and signaling pathway involvement revealed factors related to axon guidance signaling were significantly represented in the DC signature and highly expressed in others, implying functionality in promoting the large-scale and dynamic cellular rearrangements that occur during HF formation.
Comparing the new embryonic signature gene lists to those previously published for adult HF SC and niche was particularly intriguing. As DC cells are precursors of postnatal dermal papilla niche cells (Grisanti et al., 2013a), shared genes likely confer niche cell identity or promote specialized hair-inducing functions. Meanwhile, the link between embryonic Pc progenitors and adult HF SCs is less clear-cut. Fate mapping Shh-expressing cells resulted in labeled progeny that contributed to the entire postnatal HF, including bulge SCs (Levy et al., 2005), suggesting the placode cells isolated in our experiments are indeed SC precursors. If so, the distinct signature lists of Pc and HF SCs reflect the disparate activities of an expanding progenitor pool compared to a mature, reserve SC compartment.
A notable strength of our sequencing strategy was the inclusion of all cells from a single preparation of total embryonic skin, categorized in terms of defined, isolated cell types accounting for the bulk of the skin (>85%), and one population encompassing all remaining cells. While gene expression analysis revealed this fraction largely contained endothelial and smooth muscle cells, any as of yet uncharacterized dermal subpopulations were enveloped in the mix. The advantage of collecting and sequencing this population means any distinguishing gene for such a subpopulation should be exclusively confined to this fraction – and additionally ensuring that the signatures established for the purely isolated cell types are definitive markers for those cells within the total embryonic skin.
While the results of this study showcase an impressive number of new signature genes, our approach also allows a fresh and comprehensive look at previously published reports of specific genes expressed by Pc and DC in the context of the entire developing skin. Multiple studies, for example, have assessed the crucial role of the Wnt signaling pathway in promoting hair induction, formation and growth. One such investigation systematically localized Wnt ligand expression in E14.5 skin using tissue in situ hybridization, and our RNA-sequencing data on purely isolated cells are consistent with that report (Reddy et al., 2001). At the same time, our heat map pathway analysis allows us to appreciate the relative localization of other components of the Wnt signaling pathway across all skin cell types, including receptors, inhibitors and downstream effectors. This global approach enables a birds-eye view not only of where important signaling factors are produced but also where their message is received.
Pathway analyses additionally featured the enriched expression of factors related to axon guidance signaling in the DC signature gene list. Multiple receptors related to Slit-Robo signaling, including Robo1, Unc5b, and Dcc were all found highly expressed by DC cells. Considering that Slit-Robo signaling confers repulsive information, and that ligands Ntn1 and Slit2 were present in the Fb signature, this dialogue could be crucial for sorting specialized DC precursor cells from dermal fibroblasts. Neurotrophic factor Ntf3 was part of the DC signature, with receptor Ntrk3 enriched in Sch cells. A potential role for DC cells in guiding or maintaining Sch cell development in embryonic skin has not yet been explored. Other genes related to axon guidance signaling and highly expressed in DC cells were Sema6a, with Plxn ligands again enriched in Sch cells, and Ephb1 with ligand Efnbs expressed in multiple cell types. Interestingly, select axon guidance signals were also produced from the Pc, including Ntf5 and Cxcl14 that remarkably overlapped with the adult HF SC signature. Ntf5 signals through receptors Ngfr and Ntrk2, enriched in DC and Fb cells, whereas Cxcl14 binds DC signature gene Cxcr4.
Importantly, a main aim of this project is to share the gained information on a companion website to promote additional future discovery (Figure 7). This RNA-sequencing endeavor was conducted to provide a global tissue-level look at gene expression with cellular resolution encompassing several key cell types present in E14.5 skin, as hair follicle formation is just starting. Seemingly any gene identified here as specifically enriched in one cell type could be crucial to its unique functionality but the task of investigating all signature genes is too great for any one lab to undertake alone. To promote the exchange of this information, we created a website with an easily searchable interface to act as a companion to this resource paper: the Hair Gene Expression Library (http://hair-gel.net/). This website will be maintained and periodically updated as new information becomes available. Researchers can query any gene of interest to interrogate if it is present and/or where its expression is most enriched within embryonic skin. Our ultimate hope is to inspire and enable additional studies that will further elucidate the complex molecular controls driving skin and hair morphogenesis.

EXPERIMENTAL PROCEDURES

Mice
Sox2GFP and Lef1-RFP mice were described previously (Ellis et al., 2004; Rendl et al., 2005), and the animals were housed in facilities operated by the Center for Comparative Medicine and Surgery (CCMS) at Icahn School of Medicine. All animal experiments were conducted in accordance with the guidelines and approval of the Institutional Animal Care and Use Committee at Icahn School of Medicine at Mount Sinai.

Immunofluorescence staining
Whole-mounted embryo skins or sections were fixed with 4% PFA, washed with PBS, then incubated with primary and secondary antibodies. Antibodies used included E-CADHERIN (Rat, 1:100, Invitrogen), P-CADHERIN (Goat, 1:100, R&D Systems), TRP2 (DCT) (Goat, 1:100, Santa Cuz), ITGB4 (Rat, 1:100, BD Pharmigen), EDAR (goat, 1:100, R&D), KREMEN2 (Goat, 1:500, R&D Systems), SOX21 (Goat, 1:100, R&D Systems), HHIP (Goat, 1:100, R&D Systems), or FOXD1 (Goat, 1:100, Santa Cruz). After washes samples were incubated with the Rhodamine Red-X-, AF488-, or DyLight 649-conjugated donkey anti-goat, rat, or rabbit secondary antibodies (Jackson Immunoresearch). Nuclei were counterstained with DAPI. Slides were analyzed using a Leica SP5 DM confocal microscope driven by the Leica LASAF software.

Cell isolation by FACS
To isolate cells from embryonic skin, Sox2GFP/Lef1-RFP E14.5 embryos were processed as previously described (Tsai et al., 2014). Embryos were collected in ice-cold PBS and back skins harvested by microdissection. Skins were incubated in dispase (Invitrogen) with 0.2% collagenase (Sigma-Aldrich) and 20U/ul of DNAse (Roche) for 40 minutes at 37C, filtered through 40um cell strainers, and then centrifuged at 350g for 5 minutes. Resuspended cell pellets were stained with primary antibodies against E-cadherin (Rat, 1:250, Invitrogen) and P-cadherin (Goat, 1:50, R&D Systems), followed by staining for secondary antibodies Donkey anti-rat APC (1:400, Jackson Immunoresearch) and Streptavidin-Brilliant Violet 421 (1:200, Biolegend). DAPI was added for dead cell identification. Cell isolations were performed on a BD Influx cell sorter at the ISMMS flow core facility.

cDNA generation, library manufacture, and RNA-sequencing
Total RNA obtained from FACS-sorted cells was purified with the Absolutely RNA Nanoprep kit (Agilent), quantified with the NanoDrop spectrophotometer (Thermo Scientific) and measured for quality control by the Agilent Bioanalyzer. Samples with RIN (RNA integrity number) scores 8.8 and higher were further processed. 6ng starting material was reverse transcribed and amplified to 5–7μg cDNA with the RNA Ovation RNA-seq System V2 (NuGEN). From 100ng amplified cDNA sequencing libraries were generated with 16 unique barcoded adapters (8 samples × 2 replicates) using the Ovation Ultralow DR Library System (NuGEN). Library concentration and quality was quantified by Qbit (invitrogen) and Bioanalyzer (Agilent) and subsequently sequenced on the IlluminaHiSeq 2000 platform using a 100nt single-read setting. Details of RNA-sequencing analysis are provided in Supplemental Experimental Procedures.

Real-time qRT-PCR
Total RNA obtained from FACS-sorted cells was purified with the Absolutely RNA Nanoprep kit (Agilent), quantified with the NanoDrop spectrophotometer (Thermo Scientific) and reverse transcribed using oligo(dT) primers (Superscript III First-Strand Synthesis System, Invitrogen). Real-time qRT-PCR was performed with a LightCycler 480 (Roche) instrument with LightCycler DNA Master SYBR Green I reagents. Differences between samples and controls were calculated based on the 2−ΔΔCt method and normalized to Gapdh. Measurements were recorded in duplicate. Primers used are listed in Supplemental Experimental Procedures.

Supplementary Material
1234567