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.