In order to decipher the transcriptome regulatory network and cell fate decisions of DP cells during the hair follicle cycle, we isolated and dissociated dorsal skin tissue of cashmere goats from telogen (HFT) and anagen (HFA), into single cells and performed droplet-based single cell RNA seq (Fig. 1A). After deleting the low-quality cells, we detected 17 045 genes from 1000 cells in telogen and 19 371 genes from 3000 cells in anagen (Supplementary Fig. 1A). Overlapping information of cell distribution and gene expression were obtained (Supplementary Fig. 1B, C). In order to analyze the heterogeneity of cells during the hair follicle cycle, we performed tSNE clustering of all the single cells (Fig. 1B), we found that only a small fraction of cells overlapped between telogen and anagen (Fig. 1B, left), in addition, tSNE clustering analysis revealed 16 cell clusters according to their gene expression profiles (Fig. 1B, right). Subsequently, we evaluated a series of well-recognized marker genes which we clustered, revealing seven major cells: Lum was highly expressed in dermal cells [27], Pecam1 was highly expressed in endothelial cells [28], Gata3 was highly expressed in matrix cells [29], Fcer1g was highly expressed in immune cells [30], Sox9 was highly expressed in bulge cells [31], and Plp1 was highly expressed in melanocyte cells [32] (Supplementary Fig. 2A); the genes that were highly expressed in DP cells were Lef1, Msx1, Pcna, and Hoxc13 (Fig. 1C)[22, 33–35]. These analyses together completed the true characterization of different cell populations (Supplementary Fig. 2B, C, Supplementary Table 1). In order to isolate the DP cell population, we performed further expression analysis of DP marker genes (Fig. 1D, Supplementary Table 2), and found that LEF1, MSX1, and HOXC13 were highly expressed in cell clusters 0, 1, 10, and 15, and Pcna was highly expressed in cell clusters 0, 1, 3, 4, 5, 6, 7, 10, 12, 14, and 15 (Supplementary Fig. 3A). Msx1, Lef1, Hoxc13 revealed the true characterization of DP cells; using the Lef1 gene as a marker gene, we created interactive plots to further characterize DP cells (Supplementary Fig. 3B). We finally determined that cell clusters 0, 1, 10, and 15 were DP intermediate cell populations (Fig. 1E). In addition, we also revealed the expression pattern of all cell marker genes (Supplementary Fig. 3C). In summary, we successfully identified the major cell types in the telogen and anagen skin of cashmere goats and identified cell-specific characterization genes. These characteristic genes enabled true cellular heterogeneity identification using single-cell RNA sequencing.
DP Intermediate Cell Separation And Cell Fate Determination
We initially performed pseudotime ordering of all cells based on the Monocle algorithm (Fig. 2A). After analysis, we observed that the DP cell population was mainly concentrated in state2, 3, and 4 of branch point 1, and state5 of branch point 2 (Supplementary Fig. 4A). We then extracted all DP cells and performed re-pseudotime ordering analysis (Fig. 2B). DP cells had four branch points in the pseudotime ordering of hair follicle processes. Cell cluster 0 was mainly concentrated at branch points 1 and 3; cell cluster 1 was mainly concentrated at branch points 2 and 4; cell cluster 10 was concentrated at pre-branch of the branch point 4; and cell cluster 15 was concentrated at the pre-branch of branch 4 and branch point 3 (Supplementary Fig. 4B). In addition, based on pseudotime ordering, we analyzed the expression levels of four marker genes in DP cells (Fig. 2C). Pcna was highly expressed in state 1, 2, 8, and 9 of cell cluster 0; Msx1 and Lef1 were highly expressed in all states of cell clusters 0 and 1; and Hoxc13 was highly expressed in state 5 and 6 of cell clusters 1 and 10. Finally, we performed all gene expression analyses along each branch (Fig. 2D) and combined these with tissue observations (Fig. 2E) to reveal DP cell mechanisms and intermediate cell transition nodes. First, we visualized the gene expression profile (Fig. 2D), divided the highly expressed genes into four genesets, and then enriched the GO functions and the KEGG signaling pathways (Supplementary Fig. 5A, Supplementary Tables 3 and 4). The following KEGG signaling pathways were enriched: FOXO-mediated transcription, Wnt signaling pathway, NOTCH signaling pathway, MAPK signaling pathway, BMP signaling pathway and complexes such as LARC complex, and PCNA-MSH2-MSH6 complex. GO functions included: epithelial cell differentiation, epidermis development, skin development, keratinocyte differentiation, epidermal cell differentiation, and cornification. We revealed the main signaling pathways and the functions of DP cells. Epithelial cell differentiation and epidermal development revealed the source of DP cells, and immunofluorescence showed the development of embryonic hair follicles (Fig. 2F). Keratinocyte differentiation revealed the regulation of DP cells during cashmere growth. It is worth noting that we enriched the rhythmic process and melanogenesis. We speculated that hair follicle development is rhythmic and that melanocytes respond to daylight intensity as an inducer of hair follicle activation. In addition, we studied the relationship between genesets and found that there is an important interaction between genesets (Supplementary Fig. 5B), which promote the formation of DP cells, and the growth of cashmere. At the same time, we visualized the GO functions and the KEGG signal pathways (Supplementary Fig. 5C), further revealing the GO and KEGG distribution of DP cells, and visualizing the interaction network between proteins (Supplementary Fig. 6), depicting the molecular landscape of DP cell protein interactions.
The transition of intermediate cells from cell cluster 10 to cell cluster 1: A foreshadowing of hair follicle apoptosis
In order to reveal the changes in DP cell fate and cell status, we first analyzed branch point 4, named it cell fate 1 (Fig. 3C, top), and performed gene expression analysis along branch point 4 (Fig. 3A, Supplementary Table 3). For pre-branch (geneset2), we enriched the compound catabolic process, but did not enrich the functions or signaling pathways associated with hair follicle development. For geneset1, we enriched keratinization, keratinocyte differentiation, epidermal cell differentiation, epidermis development, skin development, aging, and hair cycle (Fig. 3B, top). These functions are more like the process of maintaining guard hair growth and adhering to the skin, although they were enriched with epithelial cell differentiation, this was not sufficient to demonstrate the initiation of hair follicle development. Therefore, we initially speculated that this process was a continuation of the maintenance of the old hair follicle cycle; that is, the maintenance of hair follicle state during the telogen period, but before guard hair fall out (Fig. 3E). It is worth noting that we enriched the negative regulation of cell proliferation, and cytokine-mediated signaling pathway in geneset2, which further proved our inference. In order to further reveal this process, we extracted genes associated with GO functions (Fig. 3B, bottom), and found that most of these genes were keratin gene families and keratin-related protein gene families. We then visualized these genes (Fig. 3D) and found that their expression showed a decreasing trend according to pseudo-time ordering, which revealed the cause of the gradual decrease of hair follicles in observed hair follicle tissue. It is worth noting that KRT10 was included in these genes (Fig. 3D, geneset1). In previous studies, KRT10 was found to be the marker gene for outer bulge cells, and Krt6a the marker gene for inner bulge cells (Fig. 4, geneset1). Krt6a began to decline at an earlier stage, and it can be speculated that during root sheath cell apoptosis, the inner root sheath cells undergo apoptosis first. Subsequently, we enriched all KEGG signaling pathways at branch point 4: TP53 regulates metabolic genes, aging, regulation of epidermis development, transcriptional regulation by TP53, interleukin-10 signaling, Naba secreted factors, cytokine-mediated signaling pathway, and negative regulation of cell proliferation (Supplementary Fig. 7A). In particular, the negative regulation of cell proliferation further revealed the process of hair follicle maintenance towards the end of the telogen period. Finally, we re-pseudotime ordered all the cells in branch 4 (Fig. 3C, bottom), and found two branch points, which were composed of intermediate cells 1 and 10; intermediate cells 10 evolved toward intermediate cells 1, indicating that the period was still in telogen, but had just started to evolve to anagen. We then compared the pseudotime ordering analysis of gene expression between DP and branch point 4 cells (Supplementary Fig. 7B). As we thought, these genes were highly expressed in all cells of cell fate 1, but were only highly expressed in branch point 4 of pseudotime ordering in all DP cells, the emergence of intermediate cell state 1 marks the end of the process of keeping cashmere from falling out. In general, intermediate cells 10 played an important regulatory role in the process of hair follicle maintenance, and our conclusions support the evidence that the inner bulge cells apoptosis earlier than the outer bulge cells.
Intermediate cell 1 played an important role in hair follicle apoptosis and intermediate cell cluster 0 was the key cell for hair follicle development
In order to reveal the cell fate of branch 2, we performed gene expression analysis on cell fate 2 (branch point 2; Fig. 4A) First, focusing on the KEGG signaling pathways (Supplementary Table 3), we enriched FOXO-mediated transcription, p53 signaling pathway, PID REG GR pathway, and maintenance of location in cell (Supplementary Fig. 8A), which indicated apoptotic processes in hair follicles. For geneset1, we enriched keratinization, keratinocyte differentiation, epidermal cell differentiation, epidermis development, skin development, cornification, epithelial cell differentiation, formation of the cornified envelope, aging, maintenance of location in cell, and hair cycle (Fig. 4B, top). This result is consistent with cell fate 1, indicating that the geneset1 process of cell fate was a continuation of cell fate 1; we extracted the major genes (Fig. 4B, bottom) and visualized them (Fig. 4D, geneset1). Compared with cell fate 1, cell fate 2 had 21 more genes. The expression trend of these genes was the same as that in cell fate 1; gene expression was gradually reduced according to pseudotime ordering, which further illustrated that the geneset1 process was an apoptotic process of hair follicles. It is worth noting that this process not only included the keratin gene family and the keratin-associated protein gene family, we also extracted many regulatory genes such as Csta, Dap, Cnn3, Penk, etc. These genes have similar expression patterns to those of the keratin gene family. Similarly, we visualized the Krt6a gene which is the marker gene of inner bulge cells, and the expression level showed a decreasing trend. Compared with the outer bulge cell marker gene Krt10, the decreasing trend was more obvious, except for State5; it was almost no longer expressed in other states, which indicated that the apoptosis of inner bulge cells was earlier, faster, and more thorough than that of outer bulge cells. Unfortunately, we had not enriched any GO functions or KEGG signaling pathways associated with hair follicle development in geneset2 and geneset3, but it is worth mentioning that a large number of pseudogenes were found in geneset2: Loc102169182, Loc102189548, Loc102182562, Loc102171808, Loc108634465, Loc102184991, Loc102179869, Loc108637496, Loc102175427, and Loc102180996; these pseudogenes accounted for 84% of the total genes of geneset2, and they were all ribosome-related pseudogenes. These pseudogenes attracted our interest, but we did not investigated whether they had a role at present. We speculated that in the process of apoptosis, the regulation of sRNA or genes can promote normal genes to be transcribed into pseudogenes and then inhibit their expression, thus cooperating with the process of hair follicle apoptosis. In addition, we randomly extracted several genes in geneset2 and geneset3 for visualization (Fig. 4D, geneset2 and geneset3), and found that the expression patterns of these genes were completely different from those of geneset1. The expression level was increased according to pseudotime ordering; this attracted our attention, so we extracted all the cells in cell fate 2 and performed re-pseudotime ordering (Fig. 4C). Unexpectedly, in cell fate 1, the intermediate cell cluster 10 began to evolve to 1, but they had almost disappeared in cell fate 2. Meanwhile cell cluster 1 occupied a major proportion, and had evolved to cell states 0 and 15 at the end of pseudotime ordering. Thus, we further compared the expression of these genes in the pseudotime ordering of DP cells and cell fate 2 cells (Supplementary Fig. 8B). In geneset2 and geneset3, genes were highly expressed in cell state changes (Fig. 4C, bottom), combined with signaling pathways changes; we speculated that DP cells began to receive upstream signals at this location, and that new hair follicle development occurred during the transition from telogen to anagen. This result is consistent with previous studies, that the secondary hair follicles of cashmere goats entered anagen at the end of March; our tissue observations (Fig. 2E, Mar. and Fig. 4E) also confirmed these results. In cell fate 1, we also observed cell cluster 15, and it doesn't seem to be directly related to hair follicle development. From this, we speculated that intermediate cell 1 played an important role in hair follicle apoptosis and intermediate cell 0 was the inducing factor in hair follicle development during transition from telogen to anagen and is a key cell in hair follicle development.
Cell fate 3 reveals the development of hair follicles in anagen
In order to verify the conclusion in cell fate 2, we analyzed cell fate 3. First we revealed the gene expression pattern in branch point 3 (Fig. 5A) and found that geneset1 was still dominant in the process. There were 307 genes involved in this process, accounting for 45% of the total genes of branch point 3 (Supplementary Fig. 9B). We enriched GO functions with tissue cell development, skin development, epidermal cell differentiation, keratinocyte differentiation, and cornification and keratinization (Fig. 5B, top and Supplementary Fig. 9A and Table 3), indicating that the process of old cycle cashmere growth and attachment to the skin continued. Furthermore, we enriched with formation of the cornified envelope, epidermis development, hair follicle development, myelination, gliogenesis, and hair cycle and epithelial cell differentiation; these results indicated that the new hair follicle cycle had begun and the hair follicles were beginning to develop. It is worth mentioning that we had previously enriched melanosome transport. In another study we have shown that melanocyte-receiving signals are the triggers for the initiation of the hair follicle cycle. Here, DP cells also promoted the development of the hair follicle cycle by receiving melanocyte signals. In addition, we analyzed the expression of genes in geneset1 (Fig. 5E, geneset1) and the protein interaction network (Supplementary Fig. 9F, left). The expression of these genes decreased along with pseudotime ordering, further indicating that these genes were associated with cashmere growth, the previous cashmere fall off with the decline of gene expression.. We also enriched the KEGG signaling pathways: apoptotic signaling pathway, Rab protein signal transduction, PID DELTA NP63 pathway, PID P73 pathway, PID BETA CATENIN NUC pathway, PID P53 downstream pathway, ERAD pathway, and PID MTOR 4 pathway (Supplementary Fig. 9A, geneset1). Among them, PID BETA CATENIN NUC pathway and apoptotic signaling pathway further illustrated that geneset1 was the process of hair follicle apoptosis and cashmere shedding in the last cycle. For geneset2, we enriched GO functions: DNA replication, cell cycle, synthesis of DNA, rhythmic process, connective tissue development, gland development, osteoblast differentiation, stem cell differentiation, glial cell development, tissue morphogenesis, and tissue migration; these functions marked the development of new cycle hair follicles. In particular, osteoblast differentiation and stem cell differentiation played a crucial role in the initiation of the hair follicle cycle. It is worth noting that we also enriched the rhythmic process, from which we speculated that the hair follicle cycle transition is indeed rhythmic, in addition to being controlled by genes and hormones. At the same time, we enriched the KEGG signaling pathways and major complexes of geneset2 (Supplementary Fig. 9A, geneset2). It is worth mentioning that we enriched the TNFR2 non-canonical NF-kB pathway, Wnt signaling pathway, signaling by NOTCH, and MAPK6/MAPK4 signaling. These signaling pathways further revealed the hair follicle development process, and the protein interaction network of geneset2 also supported this conclusion (Supplementary Fig. 9F, middle). We extracted 258 genes from geneset2 (Supplementary Fig. 9B) and visualized the important genes. (Fig. 5E, middle). It was found that the expression levels of these genes increased significantly along with pseudotime ordering. For geneset3, we enriched eukaryotic translation, translation, cytoplasmic translation, signaling by ROBO receptors, ncRNA processing, and TNF-alpha/NF-kappa B signaling, which provide the energy for hair follicle development, in which the highly expressed genes included Eif3e, Rpl gene family, and the Sox gene family. In order to reveal the function of cells during hair follicle cycle transition, we performed re-pseudotime ordering of cells in cell fate 3 (Fig. 5C, bottom), and found that the main cells involved in the process were cell cluster 0 and cell cluster 1 which involved crosstalk between the apoptosis of old ground hairs and the development of new hair follicles; this further explained the function of intermediate cell 1 and intermediate cell 0. The comparison of the pseudotime ordering of DP cells and cell fate 3 cells also revealed the effect of gene expression in cell cluster 0 and cell cluster 1 on old cashmere apoptosis and new hair follicle development (Supplementary Fig. 9C). In general, cell fate 3 revealed the apoptosis of old cashmere and the development of new hair follicles (Fig. 5F), especially in the new cycle, where genes, proteins, and environmental interactions contributed to the development of hair follicles.
DP cell ultimate fate: hair follicle development and new cashmere growth
Finally, we focused on cell fate 4 and performed gene expression analysis along branch 1 (Fig. 6A); we then enriched the GO functions and KEGG signaling pathways for these gene expression profiles (Fig. 6B and Supplementary Fig. 10A, B ). We visualized the relationships between these functions and signaling pathways (Fig. 6D), and examined the interaction of genes and proteins (Fig. 6E, upper-left and Supplementary Fig. 10C, D). These results revealed that cell fate 4 was the golden age of hair follicle development and was the stage of transition from hair follicle development to cashmere growth. For geneset1, we enriched the KEGG signaling pathways with cell cycle, PID FOXM1 pathway, regeneration, signaling by NOTCH3, and FoxO signaling pathway (Supplementary Fig. 10A). The genes in these signaling pathways fulfilled the following functions of the hair follicle cycle: gland development, developmental growth, epithelial cell differentiation, connective tissue development, tissue regeneration, cilium assembly, vasculature development, epithelial cell migration, and tissue migration (Fig. 6B, top and Supplementary Fig. 10A), These functions indicated the intricacy of the hair follicle structure, especially connective tissue development, tissue morphogenesis, cilium assembly, and tissue migration. Furthermore, they revealed the complexity of hair follicle development, describing the development of other tissues, including glands, connective tissue, blood vessels, etc. These tissues achieved functions such as the fixation of hair follicles, the supply of nutrients, and the function of cashmere in regulating body temperature, moisture, and metabolism. During this period, the hair follicle migrates down the dermis and connects to the capillaries to provide nutrients for cashmere growth. In addition, we visualized some important genes in geneset1 (Fig. 6E, upper right) and protein interaction networks (Supplementary Fig. 10E, geneset1 PPI). The results showed that Pcna, Plk1, Tk1, and Top2a were mainly expressed in state1 of cell cluster 0, this result revealed that the function of geneset1 was mainly achieved in cell cluster 0. We enriched the GO functions and KEGG signal pathways of geneset2, including positive regulation of cell migration, cell morphogenesis involved in differentiation, homotypic cell-cell adhesion, maintenance of cell number, cell-cell communication, regulation of cell morphogenesis involved in differentiation, stem cell population maintenance, stem cell division, and osteoblast differentiation (Supplementary Fig. 10A, geneset2). These functions illustrate the three directions of DP cells: Firstly stem cell morphology changes, cells differentiate and divide to maintain stem cell populations and thus control cashmere types. Secondly, the same type of cells adhere to each other to form cell masses and migrate upwards to be surrounded by matrix cells. Finally, DP cells communicate with other cell types to complete hair follicle development and cashmere growth. After that, we paid attention to the signaling pathway in geneset2, and the following were enriched: PI3K-Akt signaling pathway, PID AR TF pathway, signaling by receptor tyrosine kinases, response to cAMP, rhythmic process, PID BETA CATENIN NUC pathway, signaling by PDGF, and PID P38 MK2 pathway. These signaling pathways are important for hair follicle development, we also enriched the GO functions with extracellular structure organization, collagen formation, collagen fibril organization, and elastic fiber formation, which complement the hair follicle structure. In particular, appendage morphogenesis further illustrated the generation and growth of cashmere at this stage. Visualizing genes in geneset2 (Fig. 6E, bottom left) and protein interaction networks (Supplementary Fig. 10E, geneset2 PPI), we found that the expression patterns of these genes were between geneset1 and geneset3, indicating that geneset2 was the process of transition from hair follicle development to cashmere growth. So we further analyzed geneset3 and found that the following were enriched: ribosome, cytoplasmic, eukaryotic translation initiation, protein localization to endoplasmic reticulum, cytoplasmic translation, epidermis development, keratinization, keratinocyte differentiation, epidermal cell differentiation, skin development, cornification, and formation of the cornified envelope (Supplementary Fig. 10A) These GO function results were similar to cell fate 1, but the gene and protein interactions expressed in geneset3 (Supplementary Fig. 10E, geneset3 PPI) were different from cell fate 1. This result revealed the production of new cashmere and was consistent with the results from June to September in the organizational observations (Fig. 2E, June–Sept. and Fig. 6F). It is worth noting that the old cashmere falls out before the new cashmere are produced. In the re-pseudotime ordering of cell fate 4, cell cluster 1 was almost non-existent, which further predicted the function of cell cluster 0 in hair follicle development. In general, cell fate 4 illustrated the process of hair follicle development, in addition, we re-pseudotime ordered cell fate 1 cells (Fig. 6C, bottom), and revealed the process of intermediate cell 0 ending its mission and apoptosis.
In total, we have summarized and revealed the function of four different intermediate cells: Intermediate cells 10 showed important functions in the growth of cashmere and maintenance of cashmere attachment to the skin; intermediate cell 1 showed important functions in the processes of apoptosis and cashmere shedding; intermediate cell 0 initiated a new follicular cycle, completed the migration of hair follicles and the growth of cashmere. It is worth noting that cell cluster 15 existed at all stages of hair follicle development, and the number of cells was relatively stable; we therefore speculated that cell cluster 15 was a DP progenitor cell. From this, we summarized the cell fate of DP cells in the hair follicle cycle (Fig. 6G). Total DP cells can be divided into three fates in the hair follicle cycle: Telogen cell fate, anagen cell fate, and catagen cell fate. Telogen cell fate mainly describes the function of DP cells on cashmere maintenance and attachment to the epidermis; anagen cell fate reveales the function of DP cells on hair follicle development and cashmere growth; catagen cell fate describes the process by which DP cells return to their original state after completing their function (inferred only by cell state, no in-depth study has been performed).
Finally, Western blot analysis was used to detect the expression of LEF1 protein (Fig. 6H left), which indirectly verified the change rule of the number of DP cells. The results showed that the expression of LEF1 increased significantly in the transition from telogen to anagen, and continued to increase in the transition stage from anagen to catagen (Fig. 6H right), indicating that the number of DP cells increased during hair follicle development. This agreed with our single cell transcriptome sequencing results. But because the number of hair follicles increased at the same time, we could not determine the number of DP cells in each hair follicle. In general, Western blot results revealed the quantitative changes of DP cells during the development of hair follicles.