Gene Network Analysis Reveals Candidate Genes Related with the Hair Follicle Development in Sheep

doi:10.21203/rs.3.rs-1335321/v1

Download PDF

Research Article

Gene Network Analysis Reveals Candidate Genes Related with the Hair Follicle Development in Sheep

https://doi.org/10.21203/rs.3.rs-1335321/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background: Merino sheep is the most famous fine wool sheep in the world. It has high wool production and excellent wool quality, which has attracted worldwide attention. The fleece of the Merino sheep is composed predominantly of wool fibres grown from secondary wool follicles. Therefore, it is very necessary to study the development of hair follicle in order to understand the mechanism of wool. The hair follicle is a complex biological system involved in a dynamic process governed by gene regulation. The development process of hair follicle is very complex and unclear. The purpose of our research is to find candidate genes related to hair follicle development, provide molecular breeding theoretical basis for the cultivation of fine wool sheep, and provide side reference for the problems of hair loss and alopecia areata that perplex human beings.

Results: We analyzed 18 mRNA data of skin tissues of 18 Merino sheep at 4 days (E65, E85, E105 and E135) and 2 days after birth (D7 and D30). 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 proved their complex interactions. DEGs with stage-specific expression were significantly enriched in epidermal differentiation and development, hair follicle development and hair follicle morphogenesis, and 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, ITGA6, GLI2) affecting hair follicle morphogenesis was identified by network analysis.

Conclusion: This study provides a new reference of 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 of molecular breeding for the culture of fine wool sheet. These results are a valuable resource for biological investigation of fleece evolvement in animals and also supply some potential clues for understanding the molecular mechanisms of human hair development.

Hair follicle

Merino sheep

RNA-seq

Differentially expressed genes

Network

The Subo Merino (SBM) is a superfine wool-producing sheep breed in China. The wool fiber diameter is average 17–19 µm, which surpasses the standard textile count of 80 (Nm) [1], has a far-reaching impact on fine wool sheep industry. The growth and development of sheep wool is controlled by hair follicle (HF), which are tiny organs attached to the skin with a complex morphology, complex structure and periodic growth [2]. HF are composed of more than 20 kinds of cells with very intricate interactions and are involved in the regulation of HF development, HF growth, regeneration and differentiation. HF begin to develop and fully mature after birth, so the number of HF does not increase after birth.

Differentiation of HF is regulated by a variety of signaling pathways, including the BMP (bone morphogenetic protein), transforming growth factor beta (TGF-β) and Wnt signaling pathways [3, 4]. The specific expression patterns of these molecules in dermal papillae or stromal cells determine their functions during differentiation. However, knowledge regarding the cellular and molecular mechanisms is limited. After the primary HF of the sheep fetus form, the branches of the primary HF form secondary HF [5, 6]. Therefore, it is important to understand the molecular regulation mechanism of the genes. The wool quality and commercial value of fine wool sheep are determined by the structure and characteristics of their HF. To improve the wool yield and wool quality of these sheep, it is necessary to study the factors affecting the formation of HF and to deeply understand the molecular regulation mechanism of HF development. The processes of development and differentiation of HF cells are regulated by a variety of genes and multiple signaling pathways, so finding 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 not only determine the difference between wool products and other textile fibers but also determine the craft value of wool textile products [7]. Therefore, in order to rigorously improve all aspects of animal husbandry, it is important to study the development of HF in sheep, which has notable practical production value. RNA-seq can identify gene expression differences at the genome-wide level with high reproducibility, accuracy, and reliability and a wide detection range [8]. In this study, RNA-seq was used to screen genes related to HF development in order to analyze the expression regulation patterns of these genes at different stages of HF development in SBM.

Overview of the sheep skin transcriptomic data

The 18 RNA-seq libraries from the skin tissues of 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% clean reads for 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 from the six HF development stages of SBM sheep. A total of 1562 DEGs and 1762 novel DEGs were identified out in G2/G1, of which 1082 DEGs were upregulated and 480 DEGs were downregulated, while 914 novel DEGs were upregulated and 848 were downregulated, respectively (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. A total of 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 least DEGs, with only 10; this result showed that there were considerable differences from the beginning of HF development to postnatal development. The expression of DEGs were showed in heatmaps in Additional file 2: Fig. S1.the Venn analysis showed the TRPV6, MX2 and ENSOARG0000001910 were differentially expressed in all six groups (Fig. 1b).

Enrichment analysis of DEGs

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 of 7879 DEGs. In the biological processes (BP) category, the DEGs were mainly associated with the extracellular fibril organization, cardiac muscle fiber development, skeletal muscle tissue growth, regulation of microvillus assembly and positive regulation of the Wnt signaling pathway terms (Additional file 3: Fig. S2). In the cellular component category, the DEGs were mainly associated with the spectrin-associated cytoskeleton, and elastic fiber terms. In the molecular function category, the DEGs were mainly associated with the myosin heavy chain binding and monoamine transmembrane transporter activity terms (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, 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 the DEGs in adjacent comparison groups. For G2/G1, the GO terms were enriched mainly in the negative regulation of the canonical Wnt signaling pathway, hair follicle development, and 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 BP related to signal transduction and the immune response. For G6/G5, the terms were mainly enriched in collagen fiber organization, 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, establishment of the skin barrier. To sum up, DEGs in G2/G1, G6/G5 and G6/G1 are directly related to hair follicle development, especially G2/G1. G3/G2 and G4/G3 are related to muscle development, and G5/G4 is related to immunity.

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, G6/G1), Hippo signaling pathway (e.g. G2/G1), peroxisome proliferator-activated receptor (PPAR) signaling pathway (e.g. G3/G2, G4/G3), mitogen-activated protein kinase (MAPK) signaling pathway (e.g. G5/G4), and extracellular matrix (ECM)–receptor interaction (e.g. G2/G1, G5/G4, G6/G5, 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 HF, and the PPAR signaling pathway and MAPK signaling pathway were related to the stimulation of acquired HF and the HF cycle.

To explore the most potentially 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, TGFB2, TMEM79, TNF, TP63, WNT10A, and ZDHHC21) (Fig. 4a). These genes have a certain relationship with the development of the HF, skin, epidermis and 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, TGFB2, 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 we were interested in (such as WNT2, WNT3, TGFB3, KRT14, BMP2, BMP5, BMP7, BMPER, WNT16, WNT5A, SFRP2, SFRP5, SAMD3, FGF19, FZD1, and FZD5). Similar to the other 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.

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 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 BP 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), hair follicle 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, the TGF-β receptor signaling pathway (TGFB2, SMAD3 and TGFB3), collagen fiber organization (TGFB2 and SFRP2), cilium assembly and the fibroblast growth factor receptor signaling pathway (FGF19, FGF12). Surprisingly, the genes in cluster 9 play key roles in skin development, skin morphogenesis (JUP, ITGB4, DHCR24, ITGA6, ASPRV1, and SLC27A4) and skin barrier establishment (CLDN4, TMEM79, LOC101110922, and ABCA12), epithelial cell maturation(KCNE1, GJA1, TMEM79), but their expression trends were opposite, perhaps because the development of skin is a continuous process from the embryonic stage to the postnatal stage, and different genes are expressed in opposite manners at different stages despite having very similar functions, highlighting the complexity and diversity of genes. In conclusion, the genes in clusters 4, 5 and 9 were related to HF development.

Weighted Gene Co-expression Network Analysis (WGCNA)

We constructed a gene expression matrix consisting of 12791 DEGs from the standardized data (Fig. 6a). Based on the standard of a scale-free network, the applicable power value for this test was four. Sample expression pattern heatmap showed that genes in the blue module were the most highly expressed at G1, G2 and G3. The turquoise module was most highly expressed at G4 to G6. The similarity between the red and yellow modules was very high (Fig. 6b), and the expression was highest at G1. Thus, the blue, turquoise, red and yellow modules were selected for further study. GO enrichment analysis was performed for genes in the modules. The 3 most significant terms in the BP category are shown in the figure. Different terms were enriched for each of the nine modules (see Additional file 9: Table S5). The turquoise module had the most genes (4963), which were involved in negative regulation of cell proliferation, Ras protein signal transduction (GO:0007265), and positive regulation of cell migration (GO:0030335). The genes in the blue module were mainly associated with the protein K48-linked ubiquitination (GO:0070936), cellular copper ion homeostasis (GO:0006878), and glycogen metabolic process (GO:0005977) terms. The genes in the brown module were similar to those in the green module and were enriched for skeletal system morphogenesis (GO:0048705), face morphogenesis (GO:0060325), and heart development (GO:0007507). The yellow and red module genes were involved in chromosome segregation (GO:0007059) and translation (GO:0006412). 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, TGFB2, LRP4, LRP5, FGF10, TGFB3, BMPER, SOX10, WNT16, DKK1, TMEM79, BMP4, DSG4, ZDHHC21, SOSTDC1, TP63, RBPJ, FZD3, KRT25, WNT10A, WNT3, BMP7, FZD4, ITGA6) may be involved in the control of the HF development process.

protein–protein interaction (PPI) network

We performed PPI analysis on 200 hub DEGs (Fig. 7a). However, BMPER, FGFBP1, KRT75 and KRT14 that we know may be related to hair follicles or hair are prominent. We next wanted to determine how the proteins interacted by 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 a key role in this interaction network. The above description KRT14 will become the focus of our attention.

Validation of the RNA-seq data

We compared RT–PCR results with 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 for 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.

Characterization of the molecular mechanisms controlling differentiation and proliferation in mammalian HF 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. [9] studied the morphogenesis of sheep HF at six developmental stages using the H&E dyeing method and showed the asynchronous development of sheep HF. Liu et al. [10] studied the differential expression of genes during HF growth and development in Aohan fine wool sheep. They concluded that 83 genes w were differentially expressed in Aohan fine wool sheep. Nai [11] assembled and analyzed transcriptome data for cashmere goat skin and HF and obtained 617 DEGs in primary and secondary HF. Liu [12] found 4581 DEGs among HF development-related genes in inner Mongolian cashmere goats. Gao [13] 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 HF 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 [14].

According to the K-means and WGCNA analysis, we found several DEGs, such as, FZD1[15, 16], GLI2[17–19], KRT25[20, 21], LAMA5[22, 23], LRP4[24, 25], SOSTDC1[26, 27], TGFβ2[28–31], TMEM79[32], BMP7[33–37], WNT10A[38], ZDHHC21[39], SOX10[40, 41], ITGB4[42], and ITGA6[43–45], which are associated with the development of the epidermis and HF [46, 47]. They are also related to the development of the epithelium. Functional enrichment analysis showed that the DEGs were significantly enriched for negative regulation of the canonical Wnt signaling pathway, HF development, negative regulation of the BMP signaling pathway, establishment of the skin barrier, positive regulation of epidermal cell differentiation and skin development, highlighting the central roles of these DEGs in hair morphogenesis. Notably, it has been shown that the fate of HF is affected by typical Wnt/β-catenin signaling, BMP signaling, the TGF-β signaling pathway, the PI3K-Akt signaling pathway, and the Hippo signaling pathway.

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) [48]. Some studies based on HF gene regulatory networks have shown that the Wnt [3], TGF-β [49, 50], MAPK [51], Shh [52], BMP[53], PI3K-Akt[54, 55], Notch[56] and JAK-STAT[57] 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 [58], the MAPK signaling pathway plays a major role, while the TGF-β signaling pathway [59] is the lead activator of other signaling pathways; the Shh, Notch and JAK-STAT signaling pathways play fine-tuning roles [60].

In our study, we found that 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 hair follicles and secondary acquired hair follicles. Recent studies have shown that HF morphogenesis is not only related to the proliferation and differentiation of HF cells but also affected by other cells around the HF, such as sebaceous glands and sweat glands [61–63]. PPARs are members of the nuclear hormone receptor family and have become recognized as important mediators of lipid metabolism in adipocytes and sebaceous glands [64, 65]. A group of different 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 [66, 67].

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, function verification experiments are needed to validate the candidate genes for the traits related to HF development identified in this study and evaluate their effectiveness in animal breeding. Previously, we have performed sequencing to analyze noncoding RNAs [1, 9] and methylation [68]. In future research, sequencing of sheep dermal papilla cells should be performed using single-cell RNA-seq data should be included to determine the molecular drivers of HF development. In general, the discovery of candidate genes related to sheep hair follicle development will help us to improve and utilize sheep breeds. Finally, to provide the necessary animal products for human beings is the ultimate wish of our animal husbandry workers.

In summary, we used RNA-seq to sequence the mRNAs of the six developmental stages of SBM HF and conducted bioinformatics analysis. We identified a total of 7879 DEGs and 12623 novel DEGs and identified 14 key candidate genes that are significantly related to HF morphogenesis. The DEGs selected in this study and the pathways involved in HF development (TGF-β, Wnt, Hippo, PI3K-Akt, MAPK, PPAR and Hedgehog signaling pathways). Our study provides a new reference of 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 of molecular breeding for the culture of fine wool sheet. These results are a valuable resource for biological investigation of fleece evolvement in animals and also supply some potential clues for understanding the molecular mechanisms of human hair development.

Animal selection and skin tissue collection

The SBM is a subtype the Merino breed of sheep in China and is famous for its excellent quality and yield of wool and high survival rate. The tested individuals were selected from the sheep herd located at 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 a SBM ram (3 years old; MFD, 19.0 ± 0.4 µm) and then managed at the same flock. The day of insemination was designated embryonic day 0 (E0). Embryos were collected from pregnant 12 ewes (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 cm²× 3 mm at D7 and D30. Three biological replicates were generated for each of six developmental stages representing groups 1 to 6 (G1 to G6). 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 the 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 Epicentre Ribo-Zero rRNA Removal Kit (Epicentre, 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

Cuffdiff (v2.1.1) was used to calculate the fragments per kb per million reads (FPKM) values of the mRNAs in each sample based on the TopHat BAM files and the reference GTF file [69]. These counts were used in analyses of differential gene expression using edgeR v3.16.0 software [70]. mRNAs that were differentially expressed between any two comparison groups of samples representing specific developmental stages were identified using edgeR. For biological replicates, mRNAs with a P-adjust ≤ 0.05 and an absolute value log2(fold change) value ≥ 2 were considered differentially expressed between the adjacent different comparison groups (comparisons: G2/G1, G3/G2 G4/G3, G5/G4, G6/G5, and G6/G1) [71]. The expression patterns of differentially expressed mRNAs were analyzed by systematic clustering to explore the similarities and relationships between the different libraries. Furthermore, the differentially expressed mRNAs 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 based on the 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 false discovery rates (FDR) [72]. Functional evidence was obtained on the basis of the relationships between the significant GO terms (FDR = 0.05) and the DEGs. An FDR threshold of 0.05 for genes involved in BPs or metabolic pathways was applied to filter significant candidate genes by function, pathway, and network effects [73]. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) 6.8 (https://david.ncifcrf.gov/tools.jsp) was used to perform functional annotation for the genes that were significantly (P < 0.05) enriched in the expression patterns [74].

Construction of a PPI network

Following 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 a confidence score > 0.5 were selected to build the PPI network and were using MCC algorithm to calculate hub genes and select top 120 visualized using Cytoscape software.

K-means analysis and WGCNA

DEGs was clustered with the R package K-means function where K = 10 within the cluster package according to the Euclidean distance. WGCNA [75, 76] was applied to the FPKM expression data. A co-expression network was constructed with a beta value of 4. We calculated the coefficients of gene dissimilarity, performed hierarchical clustering on the genes and then determined the gene modules by the dynamic tree cut method. Through clustering analysis, the modules close to each other were merged into new modules. Before WGCNA, we identified and filtered the selected gene set. We removed the low-quality genes that had 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 (https://david.ncifcrf.gov/) [77]. The highly connected genes in each module, 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 co-expression network. These genes can act as genetic buffers to reduce the impact 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 3.8.0 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 for RT–PCR are listed see Additional file 11: Table S7. Total RNA from samples that were 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 the RT–PCR were 95°C for 30 sec followed by 40 cycles of 95°C for 5 sec and 60°C for 30 sec. 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.

Hair follicle

DEGs

Differentially expressed genes

WGCNA

Weighted gene co-expression network analysis

PPI

Protein-protein interaction

SBM

Subo merino

BMP

Bone morphogenetic protein

TGF-β

Transforming growth factor beta

Gene ontology

KEGG

Kyoto Encyclopedia of Genes and Genomes

CAM

Cell adhesion molecule

HCM

Hypertrophic cardiomyopathy

ECM

Extracellular matrix

PPAR

Peroxisome proliferator-activated receptor

MAPK

Mitogen-activated protein kinase

mRNA

Messenger RNA

RT–PCR

Reverse transcription-polymerase chain reaction

Newborn

PBS

Phosphate-buffered saline

RIN

RNA integrity number

FPKM

Fragments per kilobase of exon model per million mapped fragments

Fold change

FDR

False discovery rates

DAVID

Database for Annotation, Visualization, and Integrated Discovery

STRING

Search Tool for the Retrieval of Interacting Genes/Proteins.

Acknowledgements

The authors gratefully acknowledge the constructive comments from three reviewers.

Authors’ contributions

JH and BZ 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.

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.

Availability of data and materials

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

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 D7 and D30, 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 cm² × 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.

Author details

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

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

³ Animal Science and Technology College, Beijing University of Agriculture, Beijing, China.

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

⁵ 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.

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.
Galbraith H. Fundamental Hair Follicle Biology and Fine Fibre Production in Animals. Animal. 2010; 4(9):1490–509. https://doi.org/10.1017/S175173111000025X.
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.
Foitzik K, Lindner G, Mueller RS, 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.
Plikus MV, Mayer JA, Cruz D, Baker RE, Maini PK, Maxson R, et al. Cyclic Dermal BMP Signalling Regulates Stem Cell Activation During Hair Regeneration. Nature. 2008; 451(7176):340–44. https://doi.org/10.1038/nature06457.
Sander G, Simon Bawden C, Hynd PI, Nesci A, Rogers G, Powell BC. Expression of the Homeobox Gene, Barx2, in Wool Follicle Development. J Invest Dermatol. 2000; 115(4):753–56. https://doi.org/10.1046/j.1523-1747.2000.00122.x.
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.
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.
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.
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.
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.
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.
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.
Danilenko DM, Ring BD, Yanagihara D, Benson W, Wiemann B, Starnes CO, et al. 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. Am J Pathol. 1995; 147(1):145–54. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1869891/pdf/amjpathol00043-0151.
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–19. https://doi.org/10.1016/j.jdermsci.2005.12.011.
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.
Li ZJ, Nieuwenhuis E, Nien W, Zhang X, Zhang J, Puviindran V, et al. Kif7 Regulates Gli2 through Sufu-dependent and -Independent Functions During Skin Development and Tumorigenesis. Development. 2012; 139(22):4152–61. https://doi.org/10.1242/dev.081190.
Pantazi E, Gemenetzidis E, Trigiante G, Warnes G, Shan L, Mao X, et al. GLI2 Induces Genomic Instability in Human Keratinocytes by Inhibiting Apoptosis. Cell Death Dis. 2014; 5(1): e1028. https://doi.org/10.1038/cddis.2013.535.
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.
Zhang C, Li Y, Qin J, Yu C, Ma G, Chen H, et al. 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.
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.
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.
Wegner J, Loser K, Apsite G, Nischt R, Eckes B, Krieg T, et al. Laminin A5 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.
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. 2013; 140(3):583–93. https://doi.org/10.1242/dev.085118.
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.
Arikan V, Cumaogullari O, Ozgul BM, Oz FT. Investigation of SOSTDC1 Gene in Non-Syndromic Patients with Supernumerary Teeth. Medicina Oral Patología Oral y Cirugia Bucal. 2018:0–0. https://doi.org/10.4317/medoral.22520.
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.
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.
Samuelov L, Sprecher E, Tsuruta D, Tama S, Kloepper JE, Paus R. P-Cadherin Regulates Human Hair Growth and Cycling via Canonical Wnt Signaling and Transforming Growth Factor-β2. J Invest Dermatol. 2012. 132(10):2332–41. https://doi.org/10.1038/jid.2012.171.
Hibino T, Nishiyama T. Role of TGF-β2 in the Human Hair Cycle. J Dermatol Sci. 2004; 35(1):9–18. https://doi.org/10.1016/j.jdermsci.2003.12.003.
Inoue K, Aoi N, Yamauchi Y, Sato T, Suga H, Eto H, et al. TGF-beta is Specifically Expressed in Human Dermal Papilla Cells and Modulates Hair Folliculogenesis. J Cell Mol Med. 2009; 13(11–12):4643–56. https://doi.org/10.1111/j.1582-4934.2009.00739.x.
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.
Lv X, Sun W, Zou S, Chen L, Mwacharo JM, Wang J. Characteristics of the BMP7 Promoter in Hu Sheep. Animals. 2019; 9(11):874. https://doi.org/10.3390/ani9110874.
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.
Yasmin N, Bauer T, Modak M, Wagner K, Schuster C, Köffel R, et al. Identification of Bone Morphogenetic Protein 7 (BMP7) as an Instructive Factor for Human Epidermal Langerhans Cell Differentiation. J Exp Med. 2013; 210(12):2597–610. https://doi.org/10.1084/jem.20130275.
Lv X, Sun W, Zou S, Chen L, Mwacharo JM, Wang J. Characteristics of the BMP7 Promoter in Hu Sheep. Animals. 2019; 9(11):874. https://doi.org/10.3390/ani9110874.
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.
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.
Yang X, Zheng E, Ma Y, Chatterjee V, Villalba N, Breslin JW, et al. DHHC21 Deficiency Attenuates Renal Dysfunction During Septic Injury. Sci Rep. 2021; 11(1). https://doi.org/10.1038/s41598-021-89983-x.
Lezcano C, Ho J, Seethala RR. Sox10 and DOG1 Expression in Primary Adnexal Tumors of the Skin. Am J Dermatopathol Actions. 2017; 39(12):896–902. https://doi.org/10.1097/DAD.0000000000000872.
Slater NA, Googe PB. Sox10 Positive Breast Carcinoma Metastatic to the Skin. J Cutan Pathol. 2018; 45(5):373–74. https://doi.org/10.1111/cup.13114.
Wee LW, 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.
Zhang H, Zhang S, Zhao H, Qiao J, Liu S, Deng Z, et al. Ovine Hair Follicle Stem Cells Derived from Single Vibrissae Reconstitute Haired Skin. Int J Mol Sci. 2015; 16(8):17779–97. https://doi.org/10.3390/ijms160817779.
Ibrahim MR, Medhat W, El Fakahany H, Abdel Raouf H, Snyder EY. Deriving Keratinocyte Progenitor Cells and Keratinocytes from Human-Induced Pluripotent Stem Cells. Current Protocols in Stem Cell Biology. 2020; 54(1). https://doi.org/10.1002/cpsc.119.
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.
Lee YJ, Park SH, Park HR, Lee Y, Kang H, Kim JE. Mesenchymal Stem Cells Antagonize IFN-Induced Proinflammatory Changes and Growth Inhibition Effects via Wnt/beta-Catenin and JAK/STAT Pathway in Human Outer Root Sheath Cells and Hair Follicles. Int J Mol Sci. 2021; 22(9):4581. https://doi.org/10.3390/ijms22094581.
Jin H, Zou Z, Chang H, Shen Q, Liu L, Xing D. Photobio modulation Therapy for Hair Regeneration: A Synergetic Activation of beta-CATENIN in Hair Follicle Stem Cells by ROS and Paracrine WNTs. Stem Cell Rep. 2021; 16(6):1568–83. https://doi.org/10.1016/j.stemcr.2021.04.015.
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–26. https://doi.org/10.1111/exd.12785.
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.
Massagué J. How Cells Read TGF-β Signals. Nat Rev Mol Cell Bio. 2000; 1(3):169–78. https://doi.org/10.1038/35043051.
Headon DJ, Overbeek PA. Involvement of a Novel Tnf Receptor Homologue in Hair Follicle Induction. Nat Genet. 1999; 22(4):370–74. https://doi.org/10.1038/11943.
Kata B, Hamel PA. Alx4 Binding to LEF-1 Regulates N-CAM Promoter Activity. J Biol Chem. 2002; 277(2):1120–27. https://doi.org/10.1074/jbc.M109912200.
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.
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.
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.
Serrano CH, Ospina JP, Salazar PL, Cardona CN. 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.
Efrat AK, Torres IL, Schejter ED, Daniel SJ, 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–69. https://doi.org/10.1242/dev.02800.
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.
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.
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.
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.
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.
Kamberov YG, Karlsson EK, Kamberova GL, Lieberman DE, Sabeti PC, et al. A Genetic Basis of Variation in Eccrine Sweat Gland and Hair Follicle Density. Proc Natl Acad Sci U S A. 2015; 112(32):9932–37. https://doi.org/10.1073/pnas.1511680112.
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.
Smith KR, Thiboutot DM. Thematic review series: skin lipids. Antimicrobial lipids at the skin surface. J Lipid Res. 2008; 49(2):271–81. https://doi.org/10.1194/jlr.R700016-JLR200.
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.
Hanley K, Kömüves LG, Bass NM, He SS, Yan J, Crumrine D, et al. Stimulation of PPARalpha promotes epidermal keratinocyte differentiation in vivo. J Invest Dermatol. 2000; 115(3):353–60. https://doi.org/10.1046/j.1523-1747.2000.00073.x.
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.
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.
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.
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.
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–29. https://doi.org/10.1038/75556.
Cardoso TF, Cánovas A, Canela-Xandri O, González-Prendes R, 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. 2017; 7(1). https://doi.org/10.1038/srep40005.
Cánovas A, Pena RN, Gallardo D, Ramírez O, Amills M, Quintanilla R, et al. 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.
Langfelder P, Horvath S. WGCNA: An R Package for Weighted Correlation Network Analysis. 2008; 9(1):559. https://doi.org/10.1186/1471-2105-9-559.
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.
Da WH, Sherman BT, Lempicki RA. Systematic and Integrative Analysis of Large Gene Lists Using DAVID Bioinformatics Resources. Nat Protoc. 2008; 4(1):44–57. https://doi.org/10.1038/nprot.2008.211.

No competing interests reported.

Download PDF

Editorial decision: Major revision
15 Mar, 2022
Reviews received at journal
07 Mar, 2022
Reviewers agreed at journal
24 Feb, 2022
Reviewers invited by journal
24 Feb, 2022
Editor assigned by journal
24 Feb, 2022
Editor invited by journal
17 Feb, 2022
Submission checks completed at journal
17 Feb, 2022
First submitted to journal
07 Feb, 2022

You are reading this latest preprint version

Gene Network Analysis Reveals Candidate Genes Related with the Hair Follicle Development in Sheep

Status:

Version 1

Abstract

Figures

Background

Results

Overview of the sheep skin transcriptomic data

Analysis of genes associated with sheep skin

Enrichment analysis of DEGs

K-means analysis of DEGs

Weighted Gene Co-expression Network Analysis (WGCNA)

protein–protein interaction (PPI) network

Validation of the RNA-seq data

Discussion

Conclusion

Materials

Animal selection and skin tissue collection

Total RNA isolation, library construction, and sequencing

Differential expression analysis and clustering of mRNAs

Functional enrichment analysis and gene annotation

Construction of a PPI network

K-means analysis and WGCNA

Validation of RNA-seq data

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1