Comprehensive transcriptome and methylome analysis delineates the biological basis of hair follicle development and wool-related traits in Merino sheep

Abstract

Background: Characterization of the molecular mechanisms underlying hair follicle development is of paramount importance in the genetic improvement of wool-related traits in sheep and skin-related traits in humans. The Merino is the most important breed of fine-wooled sheep in the world. In this study, we systematically investigated the complexity of sheep hair follicle development by integrating transcriptome and methylome datasets from Merino sheep skin.

Results: We analysed 72 sequence datasets, including DNA methylome and the whole transcriptome of four gene types, i.e. protein-coding genes (PCGs), lncRNAs, circRNAs, and miRNAs, across four embryonic days (E65, E85, E105, and E135) and two postnatal days (P7 and P30) from the skin tissue of 18 Merino sheep. We revealed distinct expression profiles of these four gene types across six hair follicle developmental stages, and demonstrated their complex interactions with DNA methylation. PCGs with stage-specific expression or regulated by stage-specific lncRNAs, circRNAs, and miRNAs were significantly enriched in epithelial differentiation and hair follicle morphogenesis. Regulatory network and gene co-expression analyses identified key transcripts controlling hair follicle development. We further predicted transcriptional factors (e.g. KLF4, LEF1, HOXC13, RBPJ, VDR, RARA, and STAT3) with stage-specific involvement in hair follicle morphogenesis. Through integrating these stage-specific genomic features with results from genome-wide association studies (GWAS) of five wool-related traits in 7135 Merino sheep, we detected developmental stages and genes that were relevant with wool-related traits in sheep. For instance, genes that were specifically upregulated at E105 were significantly associated with most of wool-related traits. A phenome-wide association study (PheWAS) demonstrated that candidate genes of wool-related traits (e.g. SPHK1, GHR, PPP1R27, CSRP2, EEF1A2, and PTPN1) in sheep were also significantly associated with dermatological, metabolic, and immune traits in humans.

Conclusions: Our study provides novel insights into the molecular basis of hair follicle morphogenesis and will serve as a foundation to improve breeding for wool traits in sheep. It also indicates the importance of studying gene expression in the normal development of organs in understanding the genetic architecture of economically important traits in livestock. The datasets generated here are useful resources for functionally annotating the sheep genome, and for elucidating early skin development in mammals, including humans.

Keywords: DNA methylation; Developmental stage; Genome-wide association study; Hair follicle morphogenesis; Sheep; Transcriptome.

Fig. 1

Morphological observation of hair follicles across six developmental stages in sheep. a, b Haematoxylin-eosin (H&E) staining of skin for horizontal (magnification: × 40) and longitudinal (× 10) sections at embryonic day 65 (E65), E85, E105, E135, and postnatal day 7 (P7) and P30, respectively. The induction of hair follicles is initiated around E65 and characterized by crowded epidermal keratinocytes, hair placode (Pc) (red dashed lines) and dermal condensate (DC) (blue dashed lines). At E85, the number of primary follicles (PFs) increases significantly, and secondary follicles (SFs) start to form. At E105, PFs are regularly arranged, and the surrounding SFs, sebaceous gland (SG), outer root sheath (ORS), and inner root sheath (IRS) layers are clearly visible. At E135, the secondary-derived follicles branching from SFs appear. At P7, hair follicles mature with complete structure. The dermal papilla (DP), ORS, and IRS structures are clearly visible, and hair shafts emerge through the epidermis. At P30, the hair follicles enter into the mature anagen phase. DF, dermal fibroblast; EHG, elongated hair germ; TAC, transit amplifying cells; BSCP, bulge stem cell precursors. c Numbers of PF and SF, SF/PF ratio, body weight and length (average of three replicates; mean ± SD) of Merino sheep during hair follicles development

Fig. 2

General characteristics of molecular features in sheep skin tissue across developmental stages. a Principal component analysis (PCA) of all 18 samples based on the expression levels of four gene types and DNA methylation. PCG, protein-coding gene. b The distribution of gene expression and DNA methylation across number of stages. FPKM, fragments per kilobase of exon model per million mapped fragments; SRPBM, spliced reads per billion mapping; CPM, counts per million; FE, fold enrichment

Fig. 3.

Dynamic expression patterns of protein-coding genes (PCGs) during hair follicle development. a Heatmaps (left) shows the expression level (log₂(FPKM+ 1)) of top 10 upregulated stage-specific PCGs across developmental stages; heatmaps (right) shows the normalized consensus scores of significantly (FDR < 0.05) enriched Gene Ontology (GO) terms for all upregulated genes at each stage by the gene set enrichment analysis (GSEA). FPKM, fragments per kilobase of exon model per million mapped fragments. b Similar to a, but for the downregulated stage-specific PCGs. c, d Motifs of transcriptional factors (TFs) are significantly enriched in promoters of upregulated and downregulated stage-specific PCGs, respectively. e K-means clusters of all stage-specific PCGs. Genes in bold are those involved in the skin and epithelium development. Corresponding biological themes and the top enriched TF motifs are shown next to each cluster

Fig. 4

Stage-specific regulatory mechanisms of circRNAs and DNA methylation. a Heatmap shows the expression (log₂(SRPBM+ 1)) of upregulated stage-specific circRNAs across developmental stages, and the bubble plot shows the top three enriched gene ontology (GO) terms for parental genes of the upregulated circRNAs across stages. b K-means clustering of all stage-specific circRNAs. Corresponding biological processes are shown next to each cluster. c Heatmap shows the DNA methylation signals (log₂(FE + 1)) of stage-specific methylated regions (MRs) in promoters (1500 bp upstream and 500 bp downstream of transcriptional start sites, TSSs). FE, fold enrichment. The bubble plot shows the top five enriched gene ontology (GO) terms for genes with stage-specific MRs in promoter. d Distribution of DNA methylation levels around TSSs across all six skin developmental stages. e The methylation level of BRMS1 promoter decreases, whereas its expression level increases across developmental stages

Fig. 5

The regulatory networks of target genes of non-coding RNAs. a The regulatory networks of targets of three lncRNAs, i.e. TCONS_00394738, TCONS_00439958, and TCONS_00097544. b The regulatory networks of targets of three miRNAs, i.e. oar-miR-150, oar-miR-200c, and oar-miR-152. The colours of the edges are the directions (pink for upregulated; green for downregulated) of lncRNAs and miRNAs regulating target genes. The sizes of the circles correspond to the interaction degrees

Fig. 6

Competing endogenous RNA (ceRNA) network. a The top 200 circRNA-miRNA-mRNA and top 200 lncRNA-miRNA-mRNA interactions. The edges are the Pearson correlation coefficient (PCC) between genes. The sizes of the circles correspond to the connection degrees. b Top10 enriched Gene Ontology (GO) terms (biological processes, BP) in ceRNAs comparing to the genome background. The P values are computed by the hypergeometric test. c Sankey diagram of important candidate ceRNA pairs. Each rectangle represents a gene; degree of connection of each gene is directly proportional to rectangle size

Fig. 7

The weighted gene co-expression network analysis (WGCNA). a Bar chart shows numbers of each gene type detected in individual modules. Grey module includes genes that are not assigned to any co-expression modules. b Correlations between gene modules and developmental stages. The statistical significance of module-developmental stage relationship is corrected for multiple testing using the FDR method. The yellow stars denote FDR < 0.05. Each cell contains the correlation and the corresponding FDR value in bracket. c Heatmap shows the normalized gene expression for genes in the top five significant modules. Top enriched GO terms (biological process; BP), the normalized gene expression for module-enriched TFs, and the top representative sequence motif are shown next to each module

Fig. 8.

Integrative analysis of stage-specific molecular signatures and co-expression modules with GWAS signals of wool and growth traits. a The GWAS signal enrichment results across the stage-specific upregulated molecular features including PCGs, miRNA targets, circRNAs, circRNAs parental genes, lncRNAs, lncRNA targets, and DNA methylation. b Similar to a, but for stage-specific downregulated molecular features. c The GWAS signal enrichment results for 15 gene co-expression modules. Traits include mean fibre diameter (MFD), coefficient of variation of fibre diameter (CVFD), crimp number (CN), mean staple length (MSL), greasy fleece weight (GFW), and live weight (LW). Colour corresponds to enrichment degrees (− log₁₀FDR) calculated using the sum-based GWAS signal enrichment analysis, where (∗) means FDR < 0.1. d, e GHR; f PPP1R27. For left to right, the first plot is for the genetic variance explained by SNPs of SPHK1 (each dot is one SNP, and the highlighted region corresponds to gene region), the second for the expression patterns of SPHK1 across six developmental stages, the third for the expression patterns of SPHK1 across multiple tissues in sheep gene expression atlas, and the last for the results of phenome-wide association study (PheWAS) for SPHK1. e, f Similar to d, but for GHR and PPP1R27, respectively