Single-cell sequencing reveals the intermediate cell state and function of dermal papilla cells in the hair follicle cycle of cashmere goats


 Background:

Human hair loss and regeneration has stimulated interest in the natural hair cycle worldwide; however, such research is difficult because the periodicity of human or mouse hair is not visually obvious. Dermal papilla cells (DP cells) play an important role in the development of hair follicles, but knowledge of the differentiation and mechanisms of DP stem cells during transition through the hair follicle cycle are still limited, although some studies have reported that DP cells may have an intermediate cell state during differentiation, the classification and function of specific cell states are not clear.
Results:

Here, we used cashmere goats, that have obvious periodicity of hair follicles, as model animals and, based on unbiased single cell RNA sequencing, we identified and isolated DP cell data. Pseudotime ordering analysis was used to successfully construct a DP cell lineage differentiation trajectory and revealed the sequential activation of key genes, signaling pathways, and functions involved in cell fate decisions. At the same time, we analyzed the mechanisms of different cell fates 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 cells 1 revealed important functions in the process of apoptosis and cashmere shedding of secondary hair follicles; intermediate cells 0 initiated new follicular cycles and completed the migration of hair follicles and the occurrence of cashmere; and intermediate cells 15 are suggested to be DP progenitor cells.
Conclusions:

In development and apoptosis, inner bulge cells not only earlier than outer bulge cells, but occurred faster and was more thorough,this helps a deeper understanding of the role of bulge cells. Pseudogenes play another important role in function which promoted the competitive endogenous RNA (ceRNA) hybridization of pseudogenes.In different hair follicle cycles, DP cells will differentiate into different intermediate state cells and perform different functions, and the marker genes of the cells also changed. Intermediate cells 10 showed important functions in the growth of cashmere and maintenance of cashmere attachment to the skin; intermediate cells 1 revealed important functions in the process of apoptosis and cashmere shedding of secondary hair follicles; intermediate cells 0 initiated new follicular cycles and completed the migration of hair follicles and the occurrence of cashmere; and intermediate cells 15 are DP progenitor cells, this conclusion provides an unprecedented deeper understanding of the function of DP cells.


Introduction
Hair follicle development has been widely studied [1][2][3][4] and results show that they can have obvious periodicity [1,5,6]. In particular, the development of secondary hair follicles in cashmere goats can be divided into three stages: anagen, catagen, and telogen [7,8]. During the embryonic stage, the unspeci ed epidermis receives signals from the mesenchyme (" rst dermal signal") and subsequently forms layers of thickening epidermis known as placodes; these are the earliest morphological characteristics that mark the initiation of hair follicle morphogenesis [9,10]. Plaque sends out feedback signals and epithelial signals to promote the aggregation and coagulation of mesenchymal cells to form the initial structures of dermal papilla. During anagen, a dermal signal causes a series of gene and signal pathway changes, leading to hair follicle morphology [11]. Following the "second skin signal" [12], a series of signal transmissions [13] nally promote epithelial proliferation and down growth to form the hair follicle structure and to generate secondary hairs.
Dermal papilla (DP) cells are mesenchymal cells in the hair follicle, which regulate the development of the hair follicle and the growth of secondary hairs [14]. Many reports have revealed that DP cells are pluripotent stem cells [15,16]. Previous studies have shown that the number of DP cells determines the number of hair follicles and the density of secondary hairs, and have revealed the heterogeneity of DP cells in different types of hair follicles [17,18]. The hair of cashmere goats is divided into primary hair (guard hair) and secondary hair (ground hair or cashmere) [19,20]. The cashmere quality of cashmere goats has a signi cant relationship with the secondary hair follicles [21], but the speci c regulatory mechanism has not been fully explained. Therefore, the mechanism of action of DP cells on cashmere growth needs to be elucidated.
There are many complex cell types in the hair follicle structure. At present, DP cells, bulge cells, Hari germ (Hg) cells and matrix cells have been studied and described [20,[22][23][24]. Mok et al. show that during early hair follicle development, the molecular markers of speci c cell groups may change during cell state change [25], some of which are intermediate cells; however, these intermediate cells have not been clari ed. DP cells have many types of cell interference and different intermediate cell states, so research into their mechanisms has proved to be di cult.
Using single cell RNA(scRNA) seq to solve these problems, a recent study describes the status of intermediate cells formed during hair follicle development, which helps to analyze the fate of DP cells and determine the unrecognized intermediate cell status during development [26]. In the current study, we generated 4000 single cell scripts of cashmere goat dorsal skin from telogen and anagen. By using tdistribution random neighbor embedding(tSNE), we identi ed six major cell populations and isolated DP cells. Based on Monocle pseudotime ordering analysis, we successfully constructed the DP cell lineage differentiation trajectory and revealed the sequential activation of key genes, signaling pathways, and functions involved during cell fate decisions. Along with this, we analyzed the mechanism of different cell fates, which provided a molecular landscape for the action of DP cells during the initiation of the hair follicle cycle. Taken together, our data provided new insights into DP cell fate decisions during the hair follicle cycle, revealed the function of intermediate cells differentiated by DP cells, and delineated molecular information regarding the underappreciated DP cytodifferentiation stages.

Results
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 pro les (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][34][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 nally 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 identi ed the major cell types in the telogen and anagen skin of cashmere goats and identi ed cellspeci c characterization genes. These characteristic genes enabled true cellular heterogeneity identi cation 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 prebranch 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 pro le (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 corni cation. 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 immuno uorescence 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.  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 su cient 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 rst. 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 repseudotime ordered all the cells in branch 4 (  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, corni cation, epithelial cell differentiation, formation of the corni ed 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 con rmed 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 corni cation 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 corni ed 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 noncanonical 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 signi cantly 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 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 pro les (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 ful lled 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 xation 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 bril organization, and elastic ber 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, corni cation, and formation of the corni ed 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 repseudotime 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) 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 veri ed the change rule of the number of DP cells. The results showed that the expression of LEF1 increased signi cantly 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. Single-cell transcriptomics describe many complex biological processes in a detailed and comprehensive way [36], because it offers strong analytical ability regarding cell heterogeneity, which is conducive to the revelation of functions of different cell types in complex tissues. Recently, by using scRNA seq, seminal works have delineated an underappreciated intermediate pre-DP fate transition stage occurring prior to DC formation [26]. Here, we used the same technique to perform scRNA seq on hair follicles at two time points and isolate DP cells to provide a comprehensive understanding of genes, proteins, and signaling pathways.  [37]. For the signaling pathways we enriched 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; in particular, the upstream gene of the P53 signaling pathway, Tp53, predicts apoptosis of hair follicles [38,39]. In addition, we inferred that inner bulge cells had earlier expression patterns of the marker genes Krt6a [40] and Krt10 [41], and developed earlier than outer bulge cells, about which there have been few reports [42]. On the whole, our study revealed the function of intermediate cells 10 in promoting the growth of cashmere and maintaining its attachment to the skin.

Discussion
For cell fate 2, we enriched the Foxo-mediated transcription, p53 signaling pathway, PID REG GR pathway, and maintenance of location in cell. In particular, the p53 signaling pathway revealed the apoptotic process of cell fate 2. In addition, krt6a and krt10 indicated that apoptosis in inner bulge cells was not only earlier than outer bulge cells, but also occurred faster and was more thorough. It is worth mentioning that we obtained a large number of pseudogenes in geneset2, and contrary to previous views, we believe that they are useful. The process of apoptosis involves both the regulation of normal gene expression and the competitive inhibition caused by the combination of pseudogene transcripts and the genome, thus inhibiting the expression of normal genes [43,44]; they are similar to lncRNA in function. Our results support the competitive endogenous RNA (ceRNA) hypothesis.
For cell fate 3, we enriched the 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 in geneset1, which further revealed the function of the middle cell 1 in hair follicle apoptosis in cell fate 2. In addition, we enriched the TNFR2 non-canonical NF-kB pathway in geneset2, Wnt signaling pathway, signaling by NOTCH, and MAPK6/MAPK4 signaling, to reveal the function of intermediate cells 0 in hair follicle development [45,46]. These two intermediate cells revealed a crosstalk process of apoptosis and development during hair follicle transition [47]. It is amazing that one type of cell controls two processes of apoptosis and development at the same time.

Experimental Animals
The study used 24-month-old female Inner Mongolian cashmere goats from Jinlai Animal Husbandry Technology Co., Ltd. All experimental procedures designed in this study have been approved by the experimental animal management committee of Inner Mongolia Agricultural University.

Histological Analysis And Immuno uorescence Staining
Tissue blocks, isolated from the dorsal skin were xed with 4% paraformaldehyde (Leagene, Beijing, China) at 4 °C overnight. The next morning, the xed tissues were dehydrated in an ethanol solution and further incubated with xylene for 30 min. After incubation with xylene the samples were embedded in para n blocks. The embedded para n blocks were cut with a Leica RM2235 microtome (Leica, Nussloch, Germany) at a thickness of 5-7 µm and the samples were transferred to APES (ZSGB-BIO, Bejing, China) treated slides to avoid detachment.
For hematoxylin and eosin (H&E) staining, the slides were depara nized in 100% xylene solutions for 30 min and further rehydrated in an ethanol series. After rehydration the slides were stained with hematoxylin solution for 7 min followed by washing twice with distilled water for 5 min. After rinsing with 1% HCl (v/v) ethanol solution for 3-5 s, the slides were immediately washed with 45 °C water for 5 min.
Slides were dehydrated and then stained with 1% eosin ethanol solution and further rinsed with absolute ethanol solution for 10 min. Finally, the slides were mounted with neutral resin mounting medium and pictures were taken under an optical microscope.
For immuno uorescence staining analysis, we used the bond-RXm fully automated immunohistochemistry system (Leica Microsystems GmbH, Wetzlar, Germany), pal ™ 7-Color Manual IHC Kit 50 slides reagent (PerkinElmer, USA). First each Opal Fluorophore was reconstituted in 75 uL of DMSO. Before each procedure, dilute Opal Fluorophore in 1X ampli cation diluent at 1:100 was used to make an Opal Fluorophore working solution. The primary antibody dilution ratios were: Lef1 (1:500), Lum (1:500) and CD34 (1:800). Samples were baked in an oven at 65° C for 1 h and transferred to the bond-RXm system (Supplementary Table 5, and antibody information can be seen in Table 1). Pictures were taken under a Leica TCS SP5 II confocal microscope (Leica Microsystems GmbH, Wetzlar, Germany).

Single Cell Suspension Preparation
For each group, skin tissues were obtained from at least three females. To prepare the dorsal skin single cell suspension for single-cell RNA sequencing, the dorsal skin tissues were isolated through microdissection and a 0.25% trypsin/EDTA solution was used for tissue digestion at 37 o C for 30 min.
Samples were mechanically dissociated once every 10 min. An aliquot of 1 mg/ml collagenase IV (Sigma, St Louis, MO, USA) was used to digest the tissue with trypsin at 37 o C for 30 min. The skin tissues were mechanically dissociated into a single cell suspension by pipetting. Cell suspensions were then ltered through a 30 µm nylon cell strainer, to remove secondary hair debris (BD Falcon, BD Biosciences, San Jose, CA, USA) prior to single cell library construction.

Single Cell Library Preparation And Sequencing
Single cell barcoding and library preparation were performed based on the 10x Genomics single-cell RNA sequencing platform (10x Genomics, Pleasanton, CA, USA). Brie y, the single cell suspensions prepared above were immediately counted using a hemocytometer (TC20, Bio-Rad, Hercules, CA, USA) and the cell concentrations were adjusted to 1000 cells/µl prior to barcoding. To barcode the single cells with 10x Barcoded gel beads, the 10x Genomics Chromium Single Cell 3' Library and Gel Bead Kit v2 (10x Genomics Inc., Pleasanton, CA, USA, 120237) and the 10x Genomics Chromium barcoding system was used to construct a 10x barcoded cDNA library following the manufacturer's instructions. An Illumina HiSeq X Ten sequencer (Illumina, San Diego, CA, USA) was used for sequencing and pair-ended 150 bp (PE150) reads were generated for downstream analysis.
10x Sequencing Data Preprocessing CellRanger (v2.2.0) software was used for analyzing raw sequencing data according to the 10x Genomics o cial pipeline (https://support.10xgenomics.com/single-cell-geneexpression/software/pipelines/latest/what-is-cell-ranger). Brie y, the sequencing raw base call (BCL) les were rstly transformed into FASTQ les using the "cellranger mkfastq" function. The generated FASTQ les were then processed with "cell ranger count" wrapped function with the "--force-cells = 7000" argument to adjust sample size. Cell ranger count function used wrapped STAR software to align sequence to the reference genome. The output les containing gene expression matrices and barcode information of the CellRanger pipeline were then used for downstream visualization analysis.

Characterization Of Cell Clusters
Following CellRanger pipeline, quality control (QC) and cell clustering were analyzed with single-cell RNA seq Seurat software (v2.3.0) based on R environment (R version: 3.5.3, https://www.r-project.org/) following the online guide (https://satijalab.org/seurat/). We used " ltered_gene_bc_matrices" les generated by CellRanger as input les for Seurat. For each dataset, we rstly ltered cells with less than 200 unique detected genes and genes detected with less than 3 cells, then we used the "FilterCells" function to remove cells with a total number of detected genes (nGenes) less than 1750. After normalization, the variable genes for each dataset were calculated for downstream clustering assay.
To compare transcriptome pro les along three different developmental timepoints, we then merged three different datasets using the "RunMultiCCA" function implemented in Seurat. RunMultiCCA used a canonical correlation analysis to remove variation caused by sample source. After dataset alignment, we then performed a clustering analysis on the integrated dataset based on the t-distributed Stochastic Neighbor Embedding (tSNE) algorithms implemented in Seurat. To identify the genes speci cally expressed in clusters, we used the Seurat implemented "FindAllMarkers" function to calculate cluster markers and tSNE identi ed cell clusters were annotated based on previously reported canonical marker gene expression.
To subcluster cell clusters of interest for in-depth analysis and/or downstream differentiation trajectory construction, we used the Seurat implemented "SubsetData" function to extract clusters of interest. The extracted subclusters were then re-run through the Seurat pipeline, which provides higher resolution for dissecting cellular heterogeneity among particular cell types.

Constructing Single Cell Pseudotime Differentiation Trajectory
To interpret cell fate differentiation decisions, we used Monocle (v 2.4.0) to order single cells along pseudotime according to the o cial tutorial (http://cole-trapnell-lab.github.io/monoclerelease/docs/#constructing-single-cell-trajectories). To perform pseudotime ordering on particular cell types, we rstly subclustered interested cell types from the Seurat object, then, the Monocle object was constructed using the "newCellDataSet" function in Monocle. To order single cells along pseudotime, we used Seurat identi ed variable genes as ordering genes to construct single cell differentiation trajectories. The root state was set according to cell Seurat identi ed cell cluster labels and the "BEAM" function was used to calculate branch-speci c expressed genes. To plot a branch-speci c expression heatmap, we used the Monocle implemented "plot_genes_branched_heatmap" function and genes with qval < 1e-4 were regarded as input genes. Gene clusters were further divided into four clusters according to k-means.
For Western blot analysis, we diluted the extracted protein 20 times and diluted the BSA standard to prepare the standard protein. After denaturation, the protein was separated by electrophoresis, and then β -actin (1:5000) and LEF1 (1:1000) were transferred to the membrane for 90 min at 200 mA. Finally, color development exposure was carried out by an ECL color development system, and the gray value of the lm was analyzed by image-pro-plus.

Declarations
Data availability The single cell RNA sequencing data used in this research is deposited in NCBI GEO databases under accession number: GSE141284.

List Of Abbreviations
The list of abbreviations  UMI relationship for each dataset is also shown. Generally, the more UMI that was captured, the more genes were detected; right, average expression of genes. (C) left, The relationship between the number of unique molecules and the number of genes captured; middle, sum of gene expression; right, cell distribution in different periods.