Skip to main content
Download PDF

Gene network analysis reveals candidate genes related with the hair follicle development in sheep

BMC Genomics volume 23, Article number: 428 (2022) Cite this article

Abstract

Background

Merino sheep are the most famous fine wool sheep in the world. They have high wool production and excellent wool quality and have attracted worldwide attention. The fleece of the Merino sheep is composed predominantly of wool fibers grown from secondary wool follicles. Therefore, it is necessary to study the development of hair follicles to understand the mechanism of wool production. The hair follicle is a complex biological system involved in a dynamic process governed by gene regulation. The hair follicle development process is very complex and poorly understood. The purpose of our research is to identify candidate genes related to hair follicle development, provide a theoretical molecular breeding basis for the cultivation of fine wool sheep, and provide a reference for the problems of hair loss and alopecia areata that affect human beings.

Results

We analyzed mRNAs data in skin tissues of 18 Merino sheep at four embryonic days (E65, E85, E105 and E135) and two postnatal days (P7 and P30). G1 to G6 represent hair follicles developmental at six stages (i.e. E65 to P30). We identified 7879 differentially expressed genes (DEGs) and 12623 novel DEGs, revealed different expression patterns of these DEGs at six stages of hair follicle development, and demonstrated their complex interactions. DEGs with stage-specific expression were significantly enriched in epidermal differentiation and development, hair follicle development and hair follicle morphogenesis and were enriched in many pathways related to hair follicle development. The key genes (LAMA5, WNT10A, KRT25, SOSTDC1, ZDHHC21, FZD1, BMP7, LRP4, TGFβ2, TMEM79, SOX10, ITGB4, KRT14, ITGA6, and GLI2) affecting hair follicle morphogenesis were identified by network analysis.

Conclusion

This study provides a new reference for the molecular basis of hair follicle development and lays a foundation for further improving sheep hair follicle breeding. Candidate genes related to hair follicular development were found, which provided a theoretical basis for molecular breeding for the culture of fine wool sheep. These results are a valuable resource for biological investigations of fleece evolution in animals.

Peer Review reports

Background

Subo Merino (SBM) is a superfine wool-producing sheep breed in China. The average wool fiber diameter is 17–19 µm, which surpasses the standard textile count of 80 Nm [1] and has a far-reaching impact on the fine wool sheep industry. The growth and development of sheep wool is controlled by hair follicles (HFs), which are tiny organs attached to the skin that have a complex morphology, complex structure and periodic growth [2]. HFs are composed of multiple cells with very intricate interactions and are involved in the regulation of HF development, growth, regeneration and differentiation. The development of wool follicles has been described in detail for the Merino and it is well established that no new follicles are initiated after birth [3,4,5,6,7]. The first stage in follicle development is the proliferation of epidermal cells to form a placode beneath which an aggregation of dermal cells occurs and the two cell formations grow down together into the dermis. Progressively, the dermal cells move into the epithelial bud to form the pre-papilla and finally the epithelial bulb cells envelop the pre-papilla as the follicle lengthens and descends into the dermis. The stages as they occur in Merino sheep is first follicles formed are the primary follicles (PFs) followed by secondary follicles (SFs) and then secondary-derived follicles (SD) that branch from the SFs [4]. In Merino sheep, fibres from the SD constitute the bulk of the fleece. The first follicles to be initiated in the sheep fetus PFs are visible from 75 days of gestation and are producing a fibre by 90 days of gestation [3]. SFs do not appear until approximately 85 days of gestation. Some of these follicles will begin to branch (SD) at around 105 days [8]. HFs fully mature after birth; therefore, the number of HFs does not increase after birth. In du Cros et al. [9] description of the localization of epidermal growth factor immunoreactivity in sheep skin during wool follicle development, it was found that immunoreactivity was restricted to the periderm and intermediate layers of fetal epidermis at 55 d of gestation, when the first wave of wool follicles are initiated. This particular distribution persisted during subsequent development but never became associated with the basal cells of the epidermis. At approximately 105 d of gestation, however, reactions were detected in the outer root sheath as the follicles matured and in the differentiating cells of the sebaceous glands. Hutchison and Mellor [6] study the effects of maternal nutrition on the initiation of secondary wool follicles in foetal Scottish Blackface sheep. They found initiation of SFs usually takes place between about 95 and 135 days of gestation. Severe underfeeding during the first half of this period did not significantly inhibit the initiation of SFs, but severe underfeeding during the latter half of this period resulted in a significantly lower number of SFs and this number was not increased by refeeding ewes to a high level between 132 days and term. They concluded that SFs initiation is most affected by maternal undernutrition between about 115 and 135 days.

The differentiation of HFs is regulated by a variety of signaling pathways, including the bone morphogenetic protein (BMP), transforming growth factor beta (TGF-β) and Wnt signaling pathways [10, 11]. The specific expression patterns of these molecules in dermal papillae or stromal cells determine their functions during differentiation. However, knowledge regarding the corresponding cellular and molecular mechanisms is limited. After the primary HFs form in the sheep fetus, branches of the primary HFs form secondary HFs [12, 13]. Therefore, it is important to understand the associated molecular gene regulation mechanisms. The wool quality and commercial value of fine wool sheep are determined by the structure and characteristics of their HFs. This branching can be extensive and determines the final follicle population density in which about 80% of the follicles are SD follicles and several wool fibres emerge at the skin surface from the same orifice [14]. To improve the wool yield and wool quality of these sheep, it is necessary to study the factors affecting the formation of HFs and to deeply understand the molecular regulatory mechanisms of HF development. The processes of HF cell development and differentiation are regulated by a variety of genes and multiple signaling pathways; thus, identifying the major genes regulating the development and differentiation of HF cells has become the focus of research.

Wool fiber fineness, fiber length, wool bending, wool strength and hair flexibility determine not only the differences between wool products and other textile fibers but also the craft value of wool textile products [15]. Therefore, to rigorously improve all aspects of animal husbandry, it is important to study the HF development in sheep, and the results have notable practical production value. RNA sequencing (RNA-seq) can identify gene expression differences at the genome-wide level with high reproducibility, accuracy, and reliability and a wide detection range [16]. The practical potential of gene discovery in wool research is the provision of molecular markers for selective breeding and for altering wool growth and wool structure by other biological pathways such as sheep transgenesis that could lead to novel wool properties. In this study, RNA-seq was used to screen genes related to HFs development in order to analyze the expression regulation patterns of these genes at different stages of HFs development in SBM.

Results

Overview of sheep skin transcriptomic data

The 18 RNA-seq libraries from skin tissues at six developmental stages generated an average of 12 Gb paired-end clean reads with Q20 values > 95% (Additional file 1: Table S1). In the current study, an average of 84.98% of clean reads from the 18 samples were mapped in pairs. Furthermore, the reads were primarily aligned in gene regions, coding regions, and intergenic regions.

Analysis of genes associated with sheep skin

A total of 7879 DEGs and 12623 novel DEGs were identified in the six HF development stages of SBM sheep. A total of 1562 DEGs and 1762 novel DEGs were identified in G2/G1, of which 1082 DEGs were upregulated and 480 were downregulated, while 914 novel DEGs were upregulated and 848 were downregulated (Fig. 1a). When G3 was compared to G2, 2302 DEGs (including 1707 upregulated and 595 downregulated DEGs) and 2066 novel DEGs (including 1013 upregulated and 1053 downregulated novel DEGs) were detected. In total, 1520 DEGs (433 upregulated and 1087 downregulated) and 611 novel DEGs (251 upregulated and 360 downregulated) were found between G4 and G3. Between G5 and G4, 699 DEGs (433 upregulated and 266 downregulated) and 842 novel DEGs (595 upregulated and 247 downregulated) were identified. When G6 was compared to G5, 166 DEGs (36 upregulated and 130 downregulated) and 28 novel DEGs (7 upregulated and 360 downregulated) were identified. G6/G1 had the most DEGs, with 3507, and G6/G5 had the fewest DEGs, with only 10. This result indicated considerable differences between the beginning of HF development and postnatal development. The expression of DEGs is shown in heatmaps in Additional file 2: Fig. S1. Venn diagram analysis showed that TRPV6, MX2 and ENSOARG0000001910 were differentially expressed in all six groups (Fig. 1b).

Fig. 1
figure 1

a Histogram showing the numbers of DEGs and novel DEGs during HF morphogenesis. b Venn diagram of DEGs during HF morphogenesis 

Full size image

DEG enrichment analysis

To further explore the biological functions of DEGs related to HF development in SBM sheep during the fetal period, Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed on all 7879 DEGs. In the biological processes (BP) category, the DEGs were mainly associated with the term extracellular fibril organization, cardiac muscle fiber development, skeletal muscle tissue growth, regulation of microvillus assembly and positive regulation of the Wnt signaling pathway (Additional file 3: Fig. S2). In the cellular component category, the DEGs were mainly associated with the terms spectrin-associated cytoskeleton and elastic fiber. In the molecular function category, the DEGs were mainly associated with the term myosin heavy chain binding and monoamine transmembrane transporter activity (Additional file 3: Fig. S2a). It is presumed that the initiation of HF development is mainly driven by the above factors. The significant pathways were mainly concentrated in the cell adhesion molecule (CAM), complement and coagulation cascade, extracellular matrix (ECM) − receptor interaction, hypertrophic cardiomyopathy (HCM), hematopoietic cell lineage, and rheumatoid arthritis pathways (Additional file 3: Fig. S2b).

To further elucidate the functions of the DEGs, we also examined the top 10 GO functional enrichment BP terms and KEGG pathways of DEGs in adjacent comparison groups. For G2/G1, the GO terms were enriched mainly in the negative regulation of the canonical Wnt signaling pathway, HF development, and the negative regulation of the BMP signaling pathway, which are associated with HF or skin development (Fig. 2, Additional file 4: Table S2). The GO terms of G3/G2 were enriched in skeletal muscle contraction and muscle contraction. For G4/G3, the terms were mainly concentrated in heart and muscle development. For G5/G4, the terms were mainly enriched in BPs related to signal transduction and the immune response. For G6/G5, the terms were mainly enriched in collagen fiber organization and skin development. In addition, in the comparative analysis between 30 days after birth (G6) and the initial stage of HF development (G1), the DEGs were enriched in processes related to HF development, such as collagen fiber organization and establishment of the skin barrier. In summary, the DEGs in G2/G1, G6/G5 and G6/G1 were directly related to HF development, especially G2/G1. The DEGs in G3/G2 and G4/G3 were related to muscle development, and those in G5/G4 were related to immunity.

Fig. 2
figure 2

Top 10 GO terms (BP category) during HF morphogenesis

Full size image

KEGG enrichment analysis was also performed (Fig. 3, Additional file 5: Table S3). The phosphatidylinositol-3 kinase (PI3K)-Akt signaling pathway (e.g., G2/G1, G4/G3, G5/G4, G6/G5, and G6/G1), Hippo signaling pathway (e.g., G2/G1), peroxisome proliferator-activated receptor (PPAR) signaling pathway (e.g., G3/G2 and G4/G3), mitogen-activated protein kinase (MAPK) signaling pathway (e.g., G5/G4), and ECM–receptor interaction pathways (e.g., G2/G1, G5/G4, G6/G5, and G6/G1), which have been reported to play indispensable roles during HF development, were enriched. We also found that the Hippo signaling pathway was closely related to the development of primary and secondary HFs, and the PPAR signaling pathway and MAPK signaling pathways were related to the stimulation of acquired HFs and the HF cycle.

Fig. 3
figure 3

Top 10 enriched KEGG pathways during HF morphogenesis

Full size image

To explore the potential most key genes related to HF morphogenesis, we used Cytoscape to visualize the GO terms related to HF and skin development and the DEGs in KEGG pathways. We found that the enriched BP terms included HF development, the hair cycle process, the hair cycle, skin epidermis development, epidermis development, epidermal cell differentiation, skin development, epithelial cell differentiation, epithelial cell proliferation, epithelium development, and epithelium morphogenesis. The related genes included 41 genes (ACVR1B, APCDD1, BMP4, DKK1, DNASE1L2, DSG4, EGFR, FA2H, FGF10, FGF7, FOXN1, FZD3, FZD6, GLI1, HOXC13, HPSE, INHBA, KRT25, KRT27, KRT71, LAMA5, LGR4, LGR5, LHX2, LOC101115640, LRP4, NOTCH1, PTCH1, PTCH2, RBPJ, SFRP4, SHH, SMO, SNAIL, SOSTDC1, TGFβ2, TMEM79, TNF, TP63, WNT10A, and ZDHHC21) (Fig. 4a). These genes have certain relationships with the development of HFs, the skin, the epidermis and the epithelium. Similarly, in the pathway analysis, it was found that 17 of the above genes (ACVR1B, BMP4, DKK1, EGFR, FGF10, FGF7, FZD3, FZD6, GLI1, INHBA, NOTCH1, RBPJ, SFRP4, SMO, TGFβ2, TNF, and WNT10A) were enriched in the TGF-β signaling pathway, Wnt signaling pathway, Notch signaling pathway, Hippo signaling pathway, MAPK signaling pathway, VEGF signaling pathway and Hedgehog signaling pathway, and these pathways are known to be related to skin or HF development (Fig. 4b). Interestingly, in addition to the genes associated with the GO terms and KEGG pathways, many other genes were related to traits in which we were interested (such as WNT2, WNT3, TGFβ3, KRT14, BMP2, BMP5, BMP7, BMPER, WNT16, WNT5A, SFRP2, SFRP5, SAMD3, FGF19, FZD1, and FZD5). Similar to the above GO and KEGG analyses, we also carried out KEGG analysis of the pathways related to HF morphogenesis at different stages (Fig. 4c, Additional file 6: Fig. S3a-f) and visualized the genes in the pathways according to the different stages. We found that the results were similar to the above results and identified many genes related to HF morphogenesis.

Fig. 4
figure 4

a Gene network map of GO terms during HF morphogenesis. b Gene network map of KEGG pathways during HF morphogenesis. c DEGs in the PI3K-AKT, MAPK, Notch, TGF-β, Hippo, and Wnt signaling pathways during HF morphogenesis

Full size image

K-means analysis of DEGs

To identify the expression patterns of differentially expressed mRNAs during the six HF developmental stages, we performed k-means clustering of all DEGs and classified them into 10 clusters based on their expression changes (Fig. 5, Additional file 7: Fig. S4). GO analysis (see Additional file 8: Table S4) was carried out on the 10 clusters, and many of the identified BPs were related to the tissue types in which the genes were expressed. In cluster 4, we found that genes were enriched in epithelial cell differentiation (FZD1 and BMP7), HF morphogenesis (KRT71, WNT10A, KRT27, KRT25 and SOSTDC1) and sebaceous gland development (WNT10A and ZDHHC21), and their expression levels peaked in the G2 period. We also found that the genes in cluster 5 were related to BPs associated with hair fiber or capsule development in the negative regulation of the canonical Wnt signaling pathway (SOX2, SFRP2, LRP4, GLI1, SOX10, and GLI3), embryonic digit morphogenesis (SFRP2, LRP4, GLI2, and GLI3), the TGF-β receptor signaling pathway (TGFβ2, SMAD3 and TGFβ3), collagen fiber organization (TGFβ2 and SFRP2), cilium assembly and the fibroblast growth factor receptor signaling pathway (FGF19 and FGF12). Surprisingly, the genes in cluster 9 play key roles in skin morphogenesis (JUP, ITGB4, DHCR24, ITGA6, ASPRV1, and SLC27A4), establishment of skin barrier (CLDN4, TMEM79, LOC101110922, and ABCA12), and epithelial cell maturation (KCNE1, GJA1, and TMEM79).

Fig. 5
figure 5

Selected K-means clusters of DEGs corresponding to biological processes, with candidate genes shown next to each cluster. Each box is a differential mRNA cluster, each line in the box represents the expression trend of a differential mRNA, and the black bold line is the average trend of all differential mRNAs in the box 

Full size image

Weighted gene coexpression network analysis (WGCNA)

We constructed a gene expression matrix consisting of 12791 DEGs from the standardized data (Fig. 6a). The sample expression pattern heatmap showed that genes in the blue module were the most highly expressed at G1, G2 and G3. Those in the turquoise module were most highly expressed from G4 to G6. The similarity between the red and yellow modules was very high (Fig. 6b), and the expression was the highest at G1. Thus, the blue, turquoise, red and yellow modules were selected for further study. GO enrichment analysis was performed for genes in these modules. The 3 most significant terms in the BP category are shown in the figure. Different terms were enriched in each of the nine modules (see Additional file 9: Table S5). The turquoise module contained the most genes (4963), which were involved in negative regulation of cell proliferation, Ras protein signal transduction (GO:0,007,265), and positive regulation of cell migration (GO:0,030,335). The genes in the blue module were mainly associated with the term protein K48-linked ubiquitination (GO:0,070,936), cellular copper ion homeostasis (GO:0,006,878), and glycogen metabolic process (GO:0,005,977). The genes in the brown module were similar to those in the green module and were enriched in skeletal system morphogenesis (GO:0,048,705), face morphogenesis (GO:0,060,325), and heart development (GO:0,007,507). The yellow- and red-module genes were involved in chromosome segregation (GO:0,007,059) and translation (GO:0,006,412). KEGG analysis of the genes in these modules showed that the pathways with significant enrichment included gluconeogenesis, the insulin signaling pathway, the MAPK signaling pathway, the TGF-β signaling pathway, the PI3K-Akt signaling pathway, and other important pathways (see Additional file  10: Table S6). Based on the candidate genes associated with HF development, we found that in the yellow, green, red, brown and turquoise modules (Fig. 6c), the core genes (LAMA5, ACVR1B, EGFR, FZD1, ITFB4, PTCH1, WNT5A, GLI2, LGR4, PTCH2, TGFβ2, LRP4, LRP5, FGF10, TGFβ3, BMPER, SOX10, WNT16, DKK1, TMEM79, BMP4, DSG4, ZDHHC21, SOSTDC1, TP63, RBPJ, FZD3, KRT25, WNT10A, WNT3, BMP7, KRT14, FZD4, and ITGA6) may be involved in the control of the HF development process.

Fig. 6
figure 6

Results of WGCNA. a Division of gene modules and the correlations between the gene modules (y-axis) and the sample information (x-axis). The figure shows the clustering of genes, and the division of gene modules was based on this result. Branches of the same color were divided into the same gene modules. b Correlations between the gene modules and the module information and GO enrichment (top 3) analysis of the modules. In the panel, a darker color indicates a higher correlation, with red representing a positive correlation and green representing a negative correlation. c Visualization of the 150 hub genes related to HF morphogenesis in the modules. A darker color between genes indicates a higher correlation between the genes

Full size image

Protein–protein interaction (PPI) network

We performed PPI analysis on 200 hub DEGs (Fig. 7a). In the analysis, However, BMPER, FGFBP1, KRT75 and KRT14, which we know may be related to HFs or hair, were prominent. We next determined how the proteins interacted according to gene families related to HF development; thus, we identified KRT, WNT, FGF, IGF, NOTCH, BMP, TGF, and HOXC family genes for PPI analysis (Fig. 7b). We found that these genes constituted 3 interaction networks, among which KRT36 and KRT84 interacted independently, and 5 KRTAP members (KRTAP27-1, KRTAP24-1, KRTAP8-1, KRTAP15-1, and KRTAP11-1) interacted with each other. Most of the proteins had strong interactions. Interestingly, we found that KRT14 associates 17 KRT proteins with other proteins through NOTCH1 and FGF7 and that KRT14, NOTCH1 and FGF7 interact with each other, indicating that these three proteins play key roles in this interaction network. The above description KRT14 will become the focus of our attention.

Fig. 7
figure 7

a PPI network of the 200 hub DEGs during HF morphogenesis. b PPI network of KRT, WNT, FGF, IGF, NOTCH, BMP, TGF, and HOXC family genes

Full size image

Validation of the RNA-seq data

We compared the Real-time PCR (RT–PCR) results with the RNA-seq results (Fig. 8) and found that the expression levels of the BMP2, FGF5, HOXC13, IGF2, KRT14, MX2, SFRP2, KAP16 and TRPV6 genes in skin were consistent with the RNA-seq results. The results of KAP16 gene expression from G3 to G5 were not consistent with the sequencing results, but the expression trends were consistent in other periods. Therefore, we concluded that the sequencing results were accurate and reliable.

Fig. 8
figure 8

Results of RT–PCR and RNA-seq. The x-axis represents the time period, and the y-axis represents the relative expression

Full size image

Discussion

The characterization of the molecular mechanisms controlling differentiation and proliferation in mammalian HFs is central to our understanding of the regulation of normal hair growth, the foundations of hereditary hair loss diseases, and the origins of follicle-based tumors. Zhao et al. [17] studied the morphogenesis of sheep HFs at six developmental stages using the hematoxylin and eosin (H&E) staining method and showed the asynchronous development of sheep HFs. The Pc and DC started to form at E65, indicating the induction of HFs. At E85, the number of PFs increased and the SFs started to form. At E105, the SFs started to differentiate, and the number of secondary-derived follicles increased at E135. At E135, hair follicles matured with a complete structure, and majority hair shafts emerged through the epidermis. At P7 and P30, HFs entered into the anagen phase, during which the root of the hair divides rapidly, adding to the hair shaft (HS) [17]. Results are consistent with those reported by Rogers [14]. Liu et al. [18] studied the differential expression of genes during HF growth and development in Aohan fine wool sheep and concluded that 83 genes were differentially expressed. Nai [19] assembled and analyzed transcriptome data from cashmere goat skin and HFs and obtained 617 DEGs in primary and secondary HFs. Liu [20] found 4581 DEGs among HF development-related genes in inner Mongolian cashmere goats. Gao [21] f found a total of 2059 DEGs at the E60, E120, and newborn (NB) stages. In our study, a total of 7879 DEGs were identified in the six HF development stages of SBM sheep. Compared with previous studies, the current study identified more DEGs. The number of novel DEGs identified in this study was 12623. There were significantly more novel DEGs than known DEGs, probably because the test samples were newly obtained from the SBM variety and because this study is the first to investigate the development of HFs in SBM sheep. Although there have been many studies on HF development-related genes, it is necessary to conduct in-depth research on their functions [22].

In the K-means analysis and WGCNA, we found several DEGs, such as FZD1 [23, 24], GLI2, KRT25 [25, 26], LAMA5 [27, 28], LRP4 [29, 30], SOSTDC1 [31, 32], TGFβ2 [33, 34], TMEM79 [35], BMP7 [36, 37], WNT10A [38], ZDHHC21 [39], SOX10 [40, 41], ITGB4 [42], KRT14 [43], and ITGA6 [44], which are associated with the development of the epidermis and HF. They are also related to the development of the epithelium. Functional enrichment analysis showed that the DEGs were significantly enriched in negative regulation of the canonical Wnt signaling pathway, HF development, negative regulation of the BMP signaling pathway, establishment of the skin barrier, and positive regulation of epidermal cell differentiation and skin development, highlighting the central roles of these DEGs in hair morphogenesis. Notably, the fate of HFs is affected by typical Wnt/β-catenin signaling, BMP signaling, the TGF-β signaling pathway, the PI3K-Akt signaling pathway, and the Hippo signaling pathway.

Embryonic HF development and postnatal hair growth rely on intercellular communication within the epithelium and between epithelial and mesenchymal cells [19]. Rendl et al. [45] takes a novel multicolor labeling approach, using cell type–specific transgenic expression of red and green fluorescent proteins in combination with immunolabeling of specific antigens, to isolate pure populations of dermal papilla (DP) and four of its surrounding cell types: dermal fibroblast (DF), melanocytes (Mc), matrix (Mx) and outer root sheath (ORS). They found there was a remarkably high correlation between the DP and DF, and the Mx and ORS, respectively, highlighting the common mesenchymal origin of DP and DF and the close lineage relationship of Mx and ORS. The lowest correlation occurred between DP and Mx, revealing striking differences between the two populations whose signaling exchange orchestrates the dynamics of the hair growth. In recent years, with the development of single-cell sequencing technology, the classification of cells related to HF development and the analysis of differential genes have been studied [46,47,48]. Comparing the results of previous studies with the results of this study, we found the role of some genes in different HF cell populations (Fig. 9). These genes included the DP signature genes (i.e. WNT10A, SOSTDC, BMP, KRT14, FZD1, GLI2, LRP4, and SOX10), the ORS signature genes (i.e. LAMA5, WNT10A, SOSTDC1, TGFβ2, ITGB4, KRT14, FZD1, GLI2), the inner root sheath (IRS) signature genes (i.e. WNT10A, KRT14), the Mx signature genes (i.e. KRT25, SOSTDC1, ITGA6, BMP7, FZD1, WNT10A), the epidermal signature genes (i.e. TMEM79, KRT14, SOSTDC1, LRP4, TGFβ2, FZD1, ZDHHC21), the dermal signature genes (i.e. SOSTDC1, KRT14, KRT25, TGFβ2, BMP, WNT10A, ITGA6), the DC signature genes (i.e. WNT10A, SOSTDC1, BMP7, KRT14, FZD1, TGFβ2, LRP4), the Pc signature genes (i.e. WNT10A, SOSTDC1, FZD1), the keratinocyte signature genes (i.e. TGFβ2, ITGB4, KRT14). In addition, the hair disorder– associated genes (i.e. Gli2). The DP and ORS cells of the HF express Gli1 and Gli2 during anagen [49]. Närhi, et al. [32] found that Sostdc1 limits the number of developing HFs, and the size of primary hair placodes. Sostdc1 is essential for suppression of hair follicle fate in the normally hairless nipple epidermis. Suggest that functions of Sostdc1 can be largely attributed to its ability to attenuate Wnt/β-catenin signaling. Cross-talk of BMP and Wnt signaling pathways has been implicated in many aspects of biological events during embryogenesis and in adulthood. A secreted protein Wise and its orthologs (Sostdc1) have been shown to modulate Wnt signaling and also inhibit BMP signals and Wise also binds LRP4, a member of the LRP family functioning inhibitory to Wnt signals [50]. A homozygous missense mutation within KRT25 that causes autosomal recessive woolly hair in humans, which is consistent with findings in mutant mice. The identification of a KRT25 mutation as a cause of woolly hair in humans [51, 52]. All-trans-retinoic acid could inhibit hair follicle growth via inhibiting proliferation and inducing apoptosis of DPCs partially through the TGFβ2/Smad2/3 pathway [33]. Reddy et al. [24] found that expression of FZD1 in the placode correlates with expression of WNT10A and WNT10B in the placode. They suggest that canonical WNT signaling is likely activated in the placodes and DC of developing hair follicles by WNT10A and WNT10B, expressed in and secreted from epithelial cells, binding to FZD1, expressed in the placode epithelium and DC. Expression of FZD1 is detected in the DP and Mx of anagen follicles, and could potentially interact with DP and Mx cells, in particular WNT10A. In postnatal HFs in full anagen, FZD1 express in the ORS. At the germ and bulbous peg stages, WNT10A and WNT10B are expressed continuously in follicular epithelium during these later stages of morphogenesis. At the bulbous peg stage, WNT10A and WNT10B are strongly expressed in a cone of epithelial cells surrounding the DP. In rat whisker HFs WNT10A was expressed in the ORS, IRS, Mx and HS of anagen follicles [53].

Fig.9
figure 9

a Diagrams of HF development in Merino sheep. b Correspondence between different cell types and DEGs. Epi epidermal, Pc placode, DC dermal condensate, DP dermal papilla, Mx matrix, SG sebaceous glands, SwG Sweat gland, AMP Arrector pili muscle, Mc melanocytes, ORS outer root sheath, IRS inner root sheath, PF primary follicle, SF secondary follicle, SD secondary-derived follicle, HS hair shaft

Full size image

The characteristics of HF development and postpartum regeneration are significantly altered in microanatomy and cell viability experiments. HF development is controlled by a variety of signaling pathways, transcription factors and epigenetic regulatory factors (including miRNAs) [54]. Some studies based on HF gene regulatory networks have shown that the Wnt [10], TGF-β [55, 56], MAPK [57], Shh [58], BMP [59], PI3K-Akt [60, 61], Notch [62] and JAK-STAT [63] signaling pathways are widely involved in HF development, morphogenesis and circulation. The Wnt signaling pathway is capable of stimulating the conversion of HF stem cells from a resting state to a growing state. In most cases [64], the MAPK signaling pathway plays a major role, while the TGF-β signaling pathway [65] is the lead activator of other signaling pathways, and the Shh, Notch and JAK-STAT signaling pathways play fine-tuning roles [66].

In our study, we found that the PPAR signaling pathway was significantly enriched in G3/G2 and G4/G3, indicating that it plays an important role in the development stage of secondary HFs and secondary acquired HFs. Recent studies have shown that HF morphogenesis is not only related to HF cell proliferation and differentiation but also affected by other cells around the HF, such as sebaceous glands and sweat glands [67,68,69]. PPARs are members of the nuclear hormone receptor family and have become recognized as important mediators of lipid metabolism in adipocytes and sebaceous glands [70, 71]. A group of in vitro studies have demonstrated that PPARs play significant roles in cell differentiation, lipid synthesis and fatty acid uptake [65]. In addition, PPARs are engaged in the regulation of keratinocyte differentiation and the formation of functional skin barriers [72, 73].

Here, by analyzing DEGs at specific stages of development, we have provided insight into the genetic structure behind the complex biological process of HF development. However, functional verification experiments are needed to validate the candidate genes for traits related to HF development identified in this study and to evaluate their effectiveness in animal breeding. Previously, we performed sequencing to analyze noncoding RNAs [1, 17] and methylation [74]. In future research, the sequencing of skin should be performed using single-cell RNA-seq to determine the molecular drivers of HF development. In general, the discovery of candidate genes related to sheep HF development will help us improve and utilize sheep breeds. Finally, providing the necessary animal products for human beings is the ultimate goal of our animal husbandry workers.

Conclusion

In summary, we used RNA-seq to sequence the mRNAs in the six developmental stages of SBM HF and conducted bioinformatics analysis. We identified a total of 7879 DEGs and 12623 novel DEGs and determined 15 key candidate genes that are significantly related to HF morphogenesis. The DEGs selected in this study and the associated pathways (the TGF-β, Wnt, Hippo, PI3K-Akt, MAPK, PPAR and Hedgehog signaling pathways) are involved in HF development. Our study provides a new reference for the molecular basis of HF development and lays a foundation for further improving sheep HF breeding. Candidate genes related to HF development were found, which provided a theoretical basis for molecular breeding for the production of fine wool sheep. These results are a valuable resource for biological investigations of fleece evolution in animals and supply potential clues for understanding the molecular mechanisms of human hair development.

Materials

Animal selection and skin tissue collection

SBM is a subtype of the Merino breed of sheep in China and is famous for its excellent wool quality and yield and high survival rate. The tested individuals were selected from the sheep herd located at the Xinjiang Kechuang Animal Husbandry Breeding Center, Urumqi County (Latitude 43°01′08"–44°06′11"N, Longitude 86°37′56"–88°58′22"E), Xinjiang, China. Eighteen healthy SBM ewes (2–3 years old; mean fiber diameter (MFD), 18.1 ± 0.5 μm) were artificially inseminated with fresh sperm from an SBM ram (3 years old; MFD, 19.0 ± 0.4 μm) and then managed in the same flock. The day of insemination was designated embryonic day 0 (E0). Embryos were collected from 12 pregnant ewes (by electrocution followed by exsanguination) at four embryonic days (i.e., E65, E85, E105, and E135). Skin tissues were collected immediately after euthanasia. The skin tissues of postnatal lambs were collected in vivo at a depth of approximately 2 cm2 × 3 mm at P7 and P30. G1 to G6 represent hair follicles developmental at six stages (i.e. E65 to P30). All eighteen skin tissues were collected from the right mid-side region behind the shoulder blade of each individual and rinsed in 1 × phosphate-buffered saline (PBS). The samples were cut into small pieces, quickly placed into liquid nitrogen, and subsequently stored at − 80 °C for RNA-seq.

Total RNA isolation, library construction, and sequencing

Total RNA was isolated from the tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The RNA quality was verified using a 2100 Bioanalyzer RNA Nano Chip (Agilent Technologies, Santa Clara, CA, USA), and the total RNA content was measured using a NanoDrop ND-2000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). Only samples with RNA integrity number (RIN) scores > 8 were used for sequencing. Eighteen cDNA libraries were constructed using an NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer’s instructions, after which rRNA was removed using an Epicenter Ribo-Zero rRNA Removal Kit (Epicenter, USA). The purified libraries were quantified using a Qubit® 2.0 Fluorometer (Life Technologies, USA) and validated using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA) to confirm the insert size and calculate the molar concentration. Sequencing was performed using an Illumina HiSeq 2000 Genome Analyzer (Illumina Inc., San Diego, CA, USA).

Differential expression analysis and clustering of mRNAs

For known genes, we used Cuffdiff (v2.1.1) to calculate fragments per kilobyte per million reads (FPKM) for mRNAs in each sample based on TopHat BAM files and reference GTF files [75]. These counts were used to analyze differential gene expression using edgeR software [76]. For novel genes, we put together the genomic mapping results from all sequenced reads. Sequencing reads were assembled using Cufflinks (version 2.2.1) and compared with the known sheep genome (Ensemble Oar _v 3.1) using Cuffcompare to discover novel genes (relative to the original gene annotation file). Finally, HTSeq was used to quantify the expression of Novel genes. DEGs and Novel DEGs that were differentially expressed between any two comparison groups of samples representing specific developmental stages were identified using edgeR. The threshold value of P-value [77, 78] is determined by controlling FDR (false discovery rate), and the corrected P-value is Q-value. At the same time, we calculated the differential expression multiple (fold change) according to the FPKM value. The screening conditions of differential genes are as follows: Q-value ≤ 0.05 and fold change ≥ 2 were considered differentially expressed between the adjacent comparison groups (comparisons: G2/G1, G3/G2 G4/G3, G5/G4, G6/G5, and G6/G1) [79]. The expression patterns of DEGs were analyzed by systematic clustering to explore the similarities and relationships between the different libraries. Furthermore, the DEGs were subjected to K-means clustering using the Euclidean distance method associated with complete linkage on the BMK Cloud platform (https://www.biocloud.net/).

Functional enrichment analysis and gene annotation

Functional annotation was performed using the GO (http://geneontology.org) and KEGG (http://www.genome.ad.jp/kegg/) databases [80] based on GO and KEGG pathways. These analyses were conducted for the DEGs identified in each tissue at each stage of HF development, which were significantly enriched in dynamic expression patterns. The BPs and metabolic pathways significantly associated with the gene lists were determined based on their FDR [81]. Functional evidence was obtained on the basis of the relationships between the significant GO terms (FDR < 0.05) and the DEGs [82]. The Database for Annotation Visualization and Integrated Discovery (DAVID) 6.8 (https://david.ncifcrf.gov/tools.jsp) was used to perform functional annotation of the genes that were significantly enriched in the expression patterns [83].

Construction of a PPI network

Following the integration of the protein information in the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (https://string-db.org/) database with the DEGs, a PPI network of the identified DEGs was established. PPI network analysis was performed using STRING and Cytoscape software V3.8.0 (http://www.cytoscape.org/). First, STRING was utilized to analyze the correlation coefficients between genes. Interacting pairs with confidence scores > 0.5 were selected to build the PPI network, and the Matthews correlation coefficient (MCC) algorithm was used to calculate hub genes and the selected top 120 genes were visualized using Cytoscape software.

K-means analysis and WGCNA

DEGs were clustered with the R package K-means function, where K = 10 within the cluster package according to the Euclidean distance. WGCNA [84, 85] was applied to the FPKM expression data. A coexpression network was constructed with a beta value of 4. We calculated the coefficients of gene dissimilarity, performed hierarchical clustering of the genes and then determined the gene modules by the dynamic tree cut method. Through clustering analysis, modules close to each other were merged into new modules. Before WGCNA, we identified and filtered the selected gene set. We removed low-quality genes with an unstable impact on the results to improve the accuracy of network construction. The filtering criteria in this study were as follows: for each gene, the maximum count value in all samples was > 50, and the count was > 20 in at least 16 samples. The modules were functionally annotated using DAVID. Highly connected genes in each module, which are also known as hub genes, may play important roles in the module. Hub genes are conserved to a certain extent and are at the core of the gene coexpression network. These genes can act as genetic buffers to reduce the impacts of other gene mutations. We identified the top 150 hub genes in the modules that were most closely related to HF development differences, that is, the 150 genes with the highest connectivity in the modules, and used Cytoscape software to map the gene–gene interaction network to visualize the gene relationships.

Validation of RNA-seq data

Several differentially expressed mRNAs involved in HF development were selected and confirmed by RT–PCR with GAPDH as an internal reference. The primers used for RT–PCR are listed in Additional file 11: Table S7. Total RNA from the samples used for high-throughput RNA-seq was isolated and converted into cDNA using a PrimeScript™ RT Reagent Kit with gDNA Eraser (TaKaRa, Japan). RT–PCR was carried out on a CFX96™ Real-Time System (Bio–Rad, USA) using a TB Green Premix Ex Taq™ kit (TaKaRa, Japan) according to the manufacturer’s instructions. The thermal cycling conditions used in RT–PCR were 95 °C for 30 s followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. A reaction volume of 20 μL was used for RT–PCR according to the manufacturer’s protocol. The specificity of the SYBR Green PCR signal was confirmed by melting curve analysis. The RT–PCR experiments were performed in triplicate, and the average Ct value was used for further analysis. The 2−ΔΔCt method was used to determine the relative mRNA abundance.

Availability of data and materials

All RNA-seq data generated in this study were submitted to the NCBI SRA database under BioProject No. PRJNA705554.

Abbreviations

HF:

Hair follicle

DEGs:

Differentially expressed genes

SBM:

Subo merino

PF:

Primary follicle

SF:

Secondary follicle

SD:

Secondary-derived follicle

Pc:

Placode; DC: Dermal condensate

DP:

Dermal papilla

SG:

Sebaceous gland

ORS:

Outer root sheath

IRS:

Inner root sheath

DF:

Dermal fibroblast

Mx:

Matrix

Mc:

Melanocytes

Epi:

Epidermal

SwG:

Sweat gland

HS:

Hair shaft

AMP:

Arrector pili muscle

BMP:

Bone morphogenetic protein

TGF-β:

Transforming growth factor beta

GO:

Gene ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

BP:

Biological processes

CAM:

Cell adhesion molecule

HCM:

Hypertrophic cardiomyopathy

ECM:

Extracellular matrix

PPAR:

Peroxisome proliferator-activated receptor

MAPK:

Mitogen-activated protein kinase

mRNA:

Messenger RNA

NB:

Newborn

PBS:

Phosphate-buffered saline

RIN:

RNA integrity number

FPKM:

Fragments per kilobase of exon model per million mapped fragments

FC:

Fold change

FDR:

False discovery rates

DAVID:

Database for Annotation, Visualization, and Integrated Discovery

WGCNA:

Weighted gene co-expression network analysis

PPI:

Protein–protein interaction

RT–PCR:

Real-time PCR

STRING:

Search Tool for the Retrieval of Interacting Genes/Proteins

References

  1. Sulayman A, Tian K, Huang X, Tian Y, Xu X, Fu X, et al. Genome-Wide Identification and Characterization of Long Non-Coding RNAs Expressed During Sheep Fetal and Postnatal Hair Follicle Development. Sci Rep. 2019;9(1):8501. https://doi.org/10.1038/s41598-019-44600-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Galbraith H. Fundamental Hair Follicle Biology and Fine Fibre Production in Animals. Animal. 2010;4(9):1490–509. https://doi.org/10.1017/S175173111000025X.

    Article  CAS  PubMed  Google Scholar 

  3. Hardy H, Lyne A. The Pre-Natal Development of Wool Follicles in Merino Sheep. Aust J Biol Sci. 1956;9(3):421–41.

    Article  Google Scholar 

  4. Brook AH, Short BF, Lyne AG. Formation of New Wool Follicles in the Adult Sheep. Nature. 1960;185:51. https://doi.org/10.1038/185051a0.

    Article  CAS  PubMed  Google Scholar 

  5. Chapman RE, Hopkins PS, Thorburn GD. The Effects of Fetal Thyroidectomy and Thyroxine Administration on the Development of Skin and Wool Follicles of Sheep Fetuses. J Anat. 1974; 117(Pt 2):419–32. https://www.ncbi.nlm.nih.gov /pmc/articles/PMC1231415/pdf/janat00385–0188.

  6. Hutchison G, Mellor DJ. Effects of Maternal Nutrition on the Initiation of Secondary Wool Follicles in Foetal Sheep. J Comp Pathol. 1983;93(4):577–83. https://doi.org/10.1016/0021-9975(83)90064-6.

    Article  CAS  PubMed  Google Scholar 

  7. Chapman RE, Hardy MH. Effects of Intradermally Injected and Topically Applied Mouse Epidermal Growth Factor on Wool Growth, Skin and Wool Follicles of Merino Sheep. Aust J Biol Sci. 1988;41(2):261–8. https://doi.org/10.1071/bi9880261.

    Article  CAS  PubMed  Google Scholar 

  8. Hocking Edwards JE. Reduction in Wool Follicles Prior to Birth in Merino Sheep. Reprod Fert Develop. 1999;11(4–5):229–34. https://doi.org/10.1071/rd99049.

    Article  CAS  Google Scholar 

  9. du Cros DL, Isaacs K, Moore GP. Localization of Epidermal Growth Factor Immunoreactivity in Sheep Skin During Wool Follicle Development. J Invest Dermatol. 1992;98(1):109–15. https://doi.org/10.1111/1523-1747.ep12496010.

    Article  PubMed  Google Scholar 

  10. Andl T, Reddy ST, Gaddapara T, Millar SE. WNT Signals are Required for the Initiation of Hair Follicle Development. Dev Cell. 2002;2(5):643–53. https://doi.org/10.1016/S1534-5807(02)00167-3.

    Article  CAS  PubMed  Google Scholar 

  11. Foitzik K, Lindner G, Mueller Roever S, Maurer M, Botchkareva N, Botchkarev V, et al. Control of Murine Hair Follicle Regression (Catagen) by TGF-β1 in Vivo. FASEB J. 2000;14(5):752–60. https://doi.org/10.1096/fasebj.14.5.752.

    Article  CAS  PubMed  Google Scholar 

  12. Plikus MV, Mayer JA, de la Cruz D, Baker RE, Maini PK, Maxson R, Chuong C. Cyclic Dermal BMP Signalling Regulates Stem Cell Activation During Hair Regeneration. Nature. 2008;451(7176):340–4. https://doi.org/10.1038/nature06457.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Sander G, Simon Bawden C, Hynd PI, Nesci A, Rogers G, Powell BC. Expression of the Homeobox Gene, Barx2, Wool Follicle Development. J Invest Dermatol. 2000;115(4):753–6. https://doi.org/10.1046/j.1523-1747.2000.00122.x.

    Article  CAS  PubMed  Google Scholar 

  14. Rogers GE. Biology of the Wool Follicle: An Excursion into a Unique Tissue Interaction System Waiting to be Re-Discovered. Exp Dermatol. 2006;15(12):931–49. https://doi.org/10.1111/j.1600-0625.2006.00512.x.

    Article  PubMed  Google Scholar 

  15. Zhao B, Fu X, Tian K, Huang X, Di J, Bai Y, et al. Identification of SNPs and Expression Patterns of FZD3 Gene and its Effect on Wool Traits in Chinese Merino Sheep (Xinjiang Type). J Integr Agr. 2019;18(10):2351–60. https://doi.org/10.1016/S2095-3119(19)62735-8.

    Article  CAS  Google Scholar 

  16. Valerio C, Marianna A, Roberta E, Alfredo C. RNA-Seq and Human Complex Diseases: Recent Accomplishments and Future Perspectives. Eur J Hum Genet. 2013;21(2):134–42. https://doi.org/10.1038/ejhg.2012.129.

    Article  CAS  Google Scholar 

  17. Zhao B, Luo H, He J, Huang X, Chen S, Fu X, et al. Comprehensive Transcriptome and Methylome Analysis Delineates the Biological Basis of Hair Follicle Development and Wool-Related Traits in Merino Sheep. BMC Biol. 2021;19(1):197. https://doi.org/10.1186/s12915-021-01127-9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Liu N, Li H, Liu K, Yu J, Bu R, Cheng M, et al. Identification of Skin-Expressed Genes Possibly Associated with Wool Growth Regulation of Aohan Fine Wool Sheep. BMC Genet. 2014;15(1):144. https://doi.org/10.1186/s12863-014-0144-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Rile N, Liu Z, Gao L, Qi J, Zhao M, Xie Y, et al. Expression of Vimentin in Hair Follicle Growth Cycle of Inner Mongolian Cashmere Goats. BMC Genomics. 2018;19(1):38. https://doi.org/10.1186/s12864-017-4418-7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Liu Y, Wang L, Li X, Han W, Yang K, Wang H, et al. High-Throughput Sequencing of Hair Follicle Development-Related MicroRNAs in Cashmere Goat at Various Fetal Periods. Saudi J Biol Sci. 2018;25(7):1494–508. https://doi.org/10.1016/j.sjbs.2017.12.009.

    Article  PubMed  Google Scholar 

  21. Gao Y, Wang X, Yan H, Zeng J, Ma S, Niu Y, et al. Comparative Transcriptome Analysis of Fetal Skin Reveals Key Genes Related to Hair Follicle Morphogenesis in Cashmere Goats. PLoS ONE. 2016;11(3): e151118. https://doi.org/10.1371/journal.pone.0151118.

    Article  CAS  Google Scholar 

  22. Danilenko DM, Ring BD, Yanagihara D, Benson W, Wiemann B, Starnes CO, Pierce GF. Keratinocyte Growth Factor is an Important Endogenous Mediator of Hair Follicle Growth, Development, and Differentiation. Normalization of the Nu/Nu Follicular Differentiation Defect and Amelioration of Chemotherapy-Induced Alopecia. Am J Pathol. 1995;147(1):145–54.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Yang L, Yamasaki K, Shirakata Y, Dai X, Tokumaru S, Yahata Y, et al. Bone Morphogenetic Protein-2 Modulates Wnt and Frizzled Expression and Enhances the Canonical Pathway of Wnt Signaling in Normal Keratinocytes. J Dermatol Sci. 2006;42(2):111–9. https://doi.org/10.1016/j.jdermsci.2005.12.011.

    Article  CAS  PubMed  Google Scholar 

  24. Reddy ST, Andl T, Lu MM, Morrisey EE, Millar SE. Expression of Frizzled Genes in Developing and Postnatal Hair Follicles. J Invest Dermatol. 2004;123(2):275–82. https://doi.org/10.1111/j.0022-202X.2004.23215.x.

    Article  CAS  PubMed  Google Scholar 

  25. Zhang C, Li Y, Qin J, Yu C, Ma G, Chen H, Xu X. TMT-Based Quantitative Proteomic Analysis Reveals the Effect of Bone Marrow Derived Mesenchymal Stem Cell on Hair Follicle Regeneration. Front Pharmacol. 2021; 12. https://doi.org/10.3389/fphar.2021.658040.

  26. Morgan HJ, Benketah A, Olivero C, Rees E, Ziaj S, Mukhtar A, et al. Hair Follicle Differentiation-Specific Keratin Expression in Human Basal Cell Carcinoma. Clin Exp Dermatol. 2020;45(4):417–25. https://doi.org/10.1111/ced.14113.

    Article  CAS  PubMed  Google Scholar 

  27. Li J, Tzu J, Chen Y, Zhang YP, Nguyen NT, Gao J, et al. Laminin-10 is Crucial for Hair Morphogenesis. Embo J. 2003;22(10):2400–10. https://doi.org/10.1093/emboj/cdg239.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Wegner J, Loser K, Apsite G, Nischt R, Eckes B, Krieg T, et al. Laminin α5 in the Keratinocyte Basement Membrane is Required for Epidermal-Dermal Intercommunication. Matrix Biol. 2016;56:24–41. https://doi.org/10.1016/j.matbio.2016.05.001.

    Article  CAS  PubMed  Google Scholar 

  29. Ahn Y, Sims C, Logue JM, Weatherbee SD, Krumlauf R. Lrp4 and Wise Interplay Controls the Formation and Patterning of Mammary and Other Skin Appendage Placodes by Modulating Wnt Signaling. Development (Cambridge). 2013;140(3):583–93. https://doi.org/10.1242/dev.085118.

    Article  CAS  Google Scholar 

  30. Ohazama A, Johnson EB, Ota MS, Choi HJ, Porntaveetus T, Oommen S, et al. Lrp4 Modulates Extracellular Integration of Cell Signaling Pathways in Development. PLoS ONE. 2008;3(12): e4092. https://doi.org/10.1371/journal.pone.0004092.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Arikan V, Cumaogullari O, Ozgul BM, Oz FT. Investigation of SOSTDC1 Gene in Non-Syndromic Patients with Supernumerary Teeth. Med Oral Patol Oral Cir Bucal. 2018:0–0. https://doi.org/10.4317/medoral.22520.

  32. Närhi K, Tummers M, Ahtiainen L, Itoh N, Thesleff I, Mikkola ML. Sostdc1 Defines the Size and Number of Skin Appendage Placodes. Dev Biol. 2012;364(2):149–61. https://doi.org/10.1016/j.ydbio.2012.01.026.

    Article  CAS  PubMed  Google Scholar 

  33. Nan W, Li G, Si H, Lou Y, Wang D, Guo R, Zhang H. All-Trans-Retinoic Acid Inhibits Mink Hair Follicle Growth Via Inhibiting Proliferation and Inducing Apoptosis of Dermal Papilla Cells through TGF-β2/Smad2/3 Pathway. Acta Histochem. 2020;122(7): 151603. https://doi.org/10.1016/j.acthis.2020.151603.

    Article  CAS  PubMed  Google Scholar 

  34. Kim B, Yoon SK. Hairless Up-Regulates Tgf-β2 Expression via Down-Regulation of miR-31 in the Skin of “Hairpoor” (HrHp) Mice. J Cell Physiol. 2015;230(9):2075–85. https://doi.org/10.1002/jcp.24935.

    Article  CAS  PubMed  Google Scholar 

  35. Elias PM, Wakefield JS. Mechanisms of Abnormal Lamellar Body Secretion and the Dysfunctional Skin Barrier in Patients with Atopic Dermatitis. J Allergy Clin Immunol. 2014;134(4):781–91. https://doi.org/10.1016/j.jaci.2014.05.048.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Lv X, Gao W, Jin C, Wang L, Wang Y, Chen W, et al. Preliminary Study on microR-148a and microR-10a in Dermal Papilla Cells of Hu Sheep. BMC Genet. 2019;20(1):70. https://doi.org/10.1186/s12863-019-0770-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Kowtharapu B, Prakasam R, Murín R, Koczan D, Stahnke T, Wree A, et al. Role of Bone Morphogenetic Protein 7 (BMP7) in the Modulation of Corneal Stromal and Epithelial Cell Functions. Int J Mol Sci. 2018;19(5):1415. https://doi.org/10.3390/ijms19051415.

    Article  CAS  PubMed Central  Google Scholar 

  38. Wang K, Yamada S, Izumi H, Tsukamoto M, Nakashima T, Tasaki T, et al. Critical in Vivo Roles of WNT10A in Wound Healing by Regulating Collagen Expression/Synthesis in WNT10A-deficient Mice. PLoS ONE. 2018;13(3): e195156. https://doi.org/10.1371/journal.pone.0195156.

    Article  CAS  Google Scholar 

  39. Mill P, Lee AW, Fukata Y, Tsutsumi R, Fukata M, Keighren M, et al. Palmitoylation Regulates Epidermal Homeostasis and Hair Follicle Differentiation. PloS Genet. 2009;5(11): e1000748. https://doi.org/10.1371/journal.pgen.1000748.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Lezcano C, Ho J, Seethala RR. Sox10 and DOG1 Expression in Primary Adnexal Tumors of the Skin. Am J Dermatopathol. 2017;39(12):896–902. https://doi.org/10.1097/DAD.0000000000000872.

    Article  PubMed  Google Scholar 

  41. Slater NA, Googe PB. Sox10 Positive Breast Carcinoma Metastatic to the Skin. J Cutan Pathol. 2018;45(5):373–4. https://doi.org/10.1111/cup.13114.

    Article  PubMed  Google Scholar 

  42. Wee LWY, Tan EC, Bishnoi P, Ng YZ, Lunny DP, Lim HW, et al. Epidermolysis Bullosa with Pyloric Atresia Associated with Compound Heterozygous ITGB4 Pathogenic Variants: Minimal Skin Involvement but Severe Mucocutaneous Disease. Pediatr Dermatol. 2021;38(4):908–12. https://doi.org/10.1111/pde.14668.

    Article  PubMed  Google Scholar 

  43. Gunnarsson AP, Christensen R, Li J, Jensen UB. Dataset on Gene Expression Profiling of Multiple Murine Hair Follicle Populations. Data Brief. 2016;9:328–34. https://doi.org/10.1016/j.dib.2016.08.063.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Chen Z, Shen G, Tan X, Qu L, Zhang C, Ma L, et al. ID1/ID3 Mediate the Contribution of Skin Fibroblasts to Local Nerve Regeneration through Itga6 in Wound Repair. Stem Cell Transl Med. 2021;10(12):1637–49. https://doi.org/10.1002/sctm.21-0093.

    Article  CAS  Google Scholar 

  45. Rendl M, Lewis L, Fuchs E. Molecular Dissection of Mesenchymal-Epithelial Interactions in the Hair Follicle. PloS Biol. 2005;3(11): e331. https://doi.org/10.1371/journal.pbio.0030331.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Yang F, Li R, Zhao C, Che T, Guo J, Xie Y, et al. Single-Cell Sequencing Reveals the New Existence form of Dermal Papilla Cells in the Hair Follicle Regeneration of Cashmere Goats. Genomics. 2022;114(2): 110316. https://doi.org/10.1016/j.ygeno.2022.110316.

    Article  CAS  PubMed  Google Scholar 

  47. Ge W, Tan SJ, Wang SH, Li L, Sun XF, Shen W, Wang X. Single-Cell Transcriptome Profiling Reveals Dermal and Epithelial Cell Fate Decisions During Embryonic Hair Follicle Development. Theranostics. 2020;10(17):7581–98. https://doi.org/10.7150/thno.44306.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Wang S, Wu T, Sun J, Li Y, Yuan Z, Sun W. Single-Cell Transcriptomics Reveals the Molecular Anatomy of Sheep Hair Follicle Heterogeneity and Wool Curvature. Front Cell Dev Biol. 2021; 9. https://doi.org/10.3389/fcell.2021.800157.

  49. Amirfakhryan E, Davarnia B, Jeddi F, Najafzadeh N. Azelaic Acid Stimulates Catalase Activation and Promotes Hair Growth through Upregulation of Gli1 and Gli2 mRNA and Shh Protein. Avicenna J Phytomed. 2020; 10(5):460–71. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7508322/pdf/AJP-10-460.

  50. Lintern KB, Guidato S, Rowe A, Saldanha JW, Itasaki N. Characterization of Wise Protein and its Molecular Mechanism to Interact with Both Wnt and BMP Signals. J Biol Chem. 2009;284(34):23159–68. https://doi.org/10.1074/jbc.M109.025478.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Zernov NV, Skoblov MY, Marakhonov AV, Shimomura Y, Vasilyeva TA, Konovalov FA, et al. Autosomal Recessive Hypotrichosis with Woolly Hair Caused by a Mutation in the Keratin 25 Gene Expressed in Hair Follicles. J Invest Dermatol. 2016;136(6):1097–105. https://doi.org/10.1016/j.jid.2016.01.037.

    Article  CAS  PubMed  Google Scholar 

  52. Yu X, Chen F, Ni C, Zhang G, Zheng L, Zhang J, et al. A Missense Mutation within the Helix Termination Motif of KRT25 Causes Autosomal Dominant Woolly Hair/Hypotrichosis. J Invest Dermatol. 2018;138(1):230–3. https://doi.org/10.1016/j.jid.2017.08.035.

    Article  CAS  PubMed  Google Scholar 

  53. Lin C, Yuan Y, Chen X, Li H, Cai B, Liu Y, et al. Expression of Wnt/β-catenin Signaling, Stem-Cell Markers and Proliferating Cell Markers in Rat Whisker Hair Follicles. J Mol Histol. 2015;46(3):233–40. https://doi.org/10.1007/s10735-015-9616-5.

    Article  CAS  PubMed  Google Scholar 

  54. Andl T, Botchkareva NV. MicroRNAs (miRNAs) in the Control of HF Development and Cycling: The Next Frontiers in Hair Research. Exp Dermatol. 2016;24(11):821–6. https://doi.org/10.1111/exd.12785.

    Article  Google Scholar 

  55. Foitzik K, Paus R, Doetschman T, Dotto GP. The TGF-β2 Isoform is Both a Required and Sufficient Inducer of Murine Hair Follicle Morphogenesis. Dev Biol. 1999;212(2):278–89. https://doi.org/10.1006/dbio.1999.9325.

    Article  CAS  PubMed  Google Scholar 

  56. Massagué J. How Cells Read TGF-β Signals. Nat Rev Mol Cell Bio. 2000;1(3):169–78. https://doi.org/10.1038/35043051.

    Article  Google Scholar 

  57. Headon DJ, Overbeek PA. Involvement of a Novel Tnf Receptor Homologue in Hair Follicle Induction. Nat Genet. 1999;22(4):370–4. https://doi.org/10.1038/11943.

    Article  CAS  PubMed  Google Scholar 

  58. Kata B, Hamel PA. Alx4 Binding to LEF-1 Regulates N-CAM Promoter Activity. J Biol Chem. 2002;277(2):1120–7. https://doi.org/10.1074/jbc.M109912200.

    Article  CAS  Google Scholar 

  59. Botchkarev VA, Sharov AA. BMP Signaling in the Control of Skin Development and Hair Follicle Growth. Differentiation. 2010;72(9–10):512–26. https://doi.org/10.1111/j.1432-0436.2004.07209005.x.

    Article  Google Scholar 

  60. Cai B, Zheng Y, Ma S, Xing Q, Wang X, Yang B, et al. Long Non-Coding RNA Regulates Hair Follicle Stem Cell Proliferation and Differentiation through PI3K/AKT Signal Pathway. Mol Med Rep. 2018;17(4):5477–83. https://doi.org/10.3892/mmr.2018.8546.

    Article  CAS  PubMed  Google Scholar 

  61. Chen Y, Fan Z, Wang X, Mo M, Zeng SB, Xu RH, et al. PI3K/Akt Signaling Pathway is Essential for De Novo Hair Follicle Regeneration. Stem Cell Res Ther. 2020;11(1):144. https://doi.org/10.1186/s13287-020-01650-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Serrano CH, Ospina JP, Salazar PL, Cardona-Castro N. Notch Signaling Pathway Expression in the Skin of Leprosy Patients: Association with Skin and Neural Damage. Front Immunol. 2020; 11. https://doi.org/10.3389/fimmu.2020.00368.

  63. Efrat AK, Torres IL, Schejter ED, Daniel St J, Ben-Zion S. Drosophila Follicle Cells are Patterned by Multiple Levels of Notch Signaling and Antagonism Between the Notch and JAK/STAT Pathways. Development. 2007;134(6):1161–9. https://doi.org/10.1242/dev.02800.

    Article  CAS  Google Scholar 

  64. Gentile P, Garcovich S. Advances in Regenerative Stem Cell Therapy in Androgenic Alopecia and Hair Loss: Wnt pathway, Growth-Factor, and Mesenchymal Stem Cell Signaling Impact Analysis on Cell Growth and Hair Follicle Development. Cells-Basel. 2019;8(5):466. https://doi.org/10.3390/cells8050466.

    Article  CAS  Google Scholar 

  65. Paus R, Foitzik K, Welker P, Bulfone-Paus S, Eichmuller S. Transforming Growth Factor-Beta Receptor Type I and Type II Expression During Murine Hair Follicle Development and Cycling. J Invest Dermatol. 1997;109(4):518–26. https://doi.org/10.1111/1523-1747.ep12336635.

    Article  CAS  PubMed  Google Scholar 

  66. Luo K. Signaling Cross Talk Between TGF-β/Smad and Other Signaling Pathways. Csh Perspect Biol. 2017;9(1): a22137. https://doi.org/10.1101/cshperspect.a022137.

    Article  CAS  Google Scholar 

  67. Plikus MV, Baker RE, Chih-Chiang C, Clyde F, Damon DLC, Thomas A, et al. Self-Organizing and Stochastic Behaviors During the Regeneration of Hair Stem Cells. Science. 2011;332(6029):586. https://doi.org/10.1126/science.1201647.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Botchkarev VA, Yaar M, Peters EMJ, Raychaudhuri SP, Botchkareva NV, Marconi A, et al. Neurotrophins in Skin Biology and Pathology. J Invest Dermatol. 2006;126(8):1719–27. https://doi.org/10.1038/sj.jid.5700270.

    Article  CAS  PubMed  Google Scholar 

  69. Kamberov YG, Karlsson EK, Kamberova GL, Lieberman DE, Sabeti PC, Morgan BA, Tabin CJ. A Genetic Basis of Variation in Eccrine Sweat Gland and Hair Follicle Density. P Natl Acad Sci Usa. 2015;112(32):9932. https://doi.org/10.1073/pnas.1511680112.

    Article  CAS  Google Scholar 

  70. Thiboutot D. Regulation of Human Sebaceous Glands. J Invest Dermatol. 2004;123(1):1–12. https://doi.org/10.1111/j.1523-1747.2004.t01-2-.x.

    Article  CAS  PubMed  Google Scholar 

  71. Smith KR, Thiboutot DM. Thematic Review Series: Skin Lipids. Sebaceous Gland Lipids: Friend or Foe? J Lipid Res. 2008;49(2):271–81.

    Article  CAS  Google Scholar 

  72. Hanley K, Jiang Y, Crumrine D, Bass NM, Appel R, Elias PM, et al. Activators of the Nuclear Hormone Receptors PPAR alpha and FXR Accelerate the Development of the Fetal Epidermal Permeability Barrier. J Clin Invest. 1997;100(3):705–12. https://doi.org/10.1172/JCI119583.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Hanley K, Kömüves LG, Bass NM, He SS, Yan J, Crumrine D, et al. Fetal Epidermal Differentiation and Barrier Development in Vivo is Accelerated by Nuclear Hormone Receptor Activators 1. J Invest Dermatol. 1999;113(5):788–95. https://doi.org/10.1046/j.1523-1747.2000.00073.x.

    Article  CAS  PubMed  Google Scholar 

  74. Tian Y, Yang X, Du J, Zeng W, Wu W, Di J, et al. Differential Methylation and Transcriptome Integration Analysis Identified Differential Methylation Annotation Genes and Functional Research Related to Hair Follicle Development in Sheep. Front Genet. 2021; 12. https://doi.org/10.3389/fgene.2021.735827.

  75. Trapnell C, Pachter L, Salzberg SL. TopHat: Discovering Splice Junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11. https://doi.org/10.1093/bioinformatics/btp120.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Robinson MD, McCarthy DJ, Smyth GK. EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics. 2009;26(1):139–40. https://doi.org/10.1093/bioinformatics/btp616.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society. Series B (Methodological). 1995; 57(1). https://doi.org/10.2307/2346101.

  78. Benjamini Y, Yekutieli D. The Control of the False Discovery Rate in Multiple Testing Under Dependency. Ann. Statist. 2001; 29(4). https://doi.org/10.1214/aos/1013699998.

  79. Ren H, Wang G, Chen L, Jiang J, Liu L, Li N, et al. Genome-Wide Analysis of Long Non-Coding RNAs at Early Stage of Skin Pigmentation in Goats (Capra Hircus). BMC Genomics. 2016; 17(1). https://doi.org/10.1186/s12864-016-2365-3.

  80. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28:27–30. https://doi.org/10.1093/nar/28.1.27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: Tool for the Unification of Biology. The Gene Ontology Consortium Nat Genet. 2000;25(1):25–9. https://doi.org/10.1038/75556.

    Article  CAS  PubMed  Google Scholar 

  82. Cardoso TF, Cánovas A, Canela XO, González PR, Amills M, Quintanilla R. RNA-seq Based Detection of Differentially Expressed Genes in the Skeletal Muscle of Duroc Pigs with Distinct Lipid Profiles. Sci Rep-Uk. 2017; 7(1). https://doi.org/10.1038/srep40005.

  83. Cánovas A, Pena RN, Gallardo D, Ramírez O, Amills M, Quintanilla R, Moore S. Segregation of Regulatory Polymorphisms with Effects on the Gluteus Medius Transcriptome in a Purebred Pig Population. PLoS ONE. 2012;7(4): e35583. https://doi.org/10.1371/journal.pone.0035583.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Langfelder P, Horvath S. WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinformatics. 2008;9(1):559. https://doi.org/10.1186/1471-2105-9-559.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Wu Z, Hai E, Di Z, Ma R, Shang F, Wang Y, et al. Using WGCNA (Weighted Gene Co-Expression Network Analysis) to Identify the Hub Genes of Skin Hair Follicle Development in Fetus Stage of Inner Mongolia Cashmere Goat. PLoS ONE. 2020;15(12): e243507. https://doi.org/10.1371/journal.pone.0243507.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

The authors gratefully acknowledge the constructive comments from three reviewers

Funding

This work was supported by the China Agriculture Research System (Nos. CARS-39) and Agricultural science and technology innovation project of Shandong Academy of Agricultural Sciences (Nos. CXGC2021B20). The funders played no role in study design, collection, analysis, data interpretation, manuscript writing, or decision to submit the manuscript for publication.

Author information

Authors and Affiliations

  1. College of Animal Science and Technology, Gansu Agricultural University, Lanzhou, China

    Junmin He & Shuangbao Gun

  2. College of Animal Science and Technology, China Agricultural University, Beijing, China

    Bingru Zhao

  3. College of Animal Science, Xinjiang Agricultural University, Urumqi, China

    Xixia Huang, Cuiling Wu, Jingyi Mao & Jing Liu

  4. Key Laboratory of Genetics Breeding and Reproduction of the Fine Wool Sheep & Cashmere Goat in Xinjiang, Institute of Animal Science, Xinjiang Academy of Animal Sciences, Urumqi, China

    Xuefeng Fu & Yuezhen Tian

  5. Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences, Jinan, China

    Guifen Liu & Kechuan Tian

Authors
  1. Junmin He
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Bingru Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xixia Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xuefeng Fu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Guifen Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yuezhen Tian
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Cuiling Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jingyi Mao
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Jing Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Shuangbao Gun
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Kechuan Tian
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

JH and KT conceived and designed the study. JH, BZ, XH, and GL collected the samples and/or generated the data. JH, CW, JM, and JL conducted the statistical analyses. XF and YT managed the experimental sheep populations. JH, BZ, SG and KT drafted and revised the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Shuangbao Gun or Kechuan Tian.

Ethics declarations

Ethics approval and consent to participate

This study was conducted in accordance with the Basel Declaration and it is reported in accordance with ARRIVE guidelines. Collection of all tissue samples was conducted in accordance with the guidelines of the Ethics Committee of the College of Animal Science (approval number: 2006–398) of Gansu Agricultural University. All embryo samples were obtained by cesarean section after 12 ewes were electrocution followed by exsanguination. The skin tissues were collected immediately after euthanasia, used for scientific research. On postnatal days P7 and P30, lamb was observed to lose consciousness after using general anesthetic the skin tissue of postnatal lambs was collected in vivo, with a depth of about 2 cm2 × 3 mm. Then skin wounds were quickly sprinkled anti-inflammatory powder and carefully suture. We took good care of the lamb until the wound healed. All efforts were made to minimize animal suffering and discomfort.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

 Table S1. Result Statistics of mRNA Genome Comparison. Mapping ratio=Mapped reads/All reads, Mapped Unique reads: reads that have only one position in the genome.

Additional file 2:

 Fig. S1. Heatmap of DEGs during the hair follicle morphogenesis. The x-axis represents the sample number, where S1 to S3 is G1; S4 to S6 is G2; S7 to S9 is G3; S10 to S12 is G4; S13 to S15 is G5; S16 to S18 is G6.

Additional file 3:

 Fig. S2. Enrichment analysis of all DEGs. (a)Top 30 of GO enrichment; (b)Top 30 of pathway enrichment.

Additional file 4:

 Table S2. GO enrichment of DEGs.

Additional file 5:

 Table S3. KEGG pathway of DEGs.

Additional file 6: Fig. S3.

Go and KEGG network diagram during the hair follicle morphogenesis. (a-f) represents the comparison group of G1-G6. Circle is represented gene, square is represented KEGG; triangle is represented GO term.

Additional file 7: Fig. S4.

K-means clustering analysis of differentially expressed mRNAs among the six comparison groups.

Additional file 8:

 Table S4. GO analysis of K-means clusters.

Additional file 9:

 Table S5. GO analysis of module in WGCNA.

Additional file 10:

 Table S6. KEGG pathway of module in WGCNA.

Additional file 11:

 Table S7. Sequence of primers.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

He, J., Zhao, B., Huang, X. et al. Gene network analysis reveals candidate genes related with the hair follicle development in sheep. BMC Genomics 23, 428 (2022). https://doi.org/10.1186/s12864-022-08552-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08552-2

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us