Skip to main content

Analysis of DNA methylation profiles during sheep skeletal muscle development using whole-genome bisulfite sequencing



DNA methylation is an epigenetic regulatory form that plays an important role in regulating the gene expression and the tissues development.. However, DNA methylation regulators involved in sheep muscle development remain unclear. To explore the functional importance of genome-scale DNA methylation during sheep muscle growth, this study systematically investigated the genome-wide DNA methylation profiles at key stages of Hu sheep developmental (fetus and adult) using deep whole-genome bisulfite sequencing (WGBS).


Our study found that the expression levels of DNA methyltransferase (DNMT)-related genes were lower in fetal muscle than in the muscle of adults. The methylation levels in the CG context were higher than those in the CHG and CHH contexts, and methylation levels were highest in introns, followed by exons and downstream regions. Subsequently, we identified 48,491, 17, and 135 differentially methylated regions (DMRs) in the CG, CHG, and CHH sequence contexts and 11,522 differentially methylated genes (DMGs). The results of bisulfite sequencing PCR (BSP) correlated well with the WGBS-Seq data. Moreover, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional annotation analysis revealed that some DMGs were involved in regulating skeletal muscle development and fatty acid metabolism. By combining the WGBS-Seq and previous RNA-Seq data, a total of 159 overlap genes were obtained between differentially expressed genes (DEGs) and DMGs (FPKM > 10 and fold change > 4). Finally, we found that 9 DMGs were likely to be involved in muscle growth and metabolism of Hu sheep.


We systemically studied the global DNA methylation patterns of fetal and adult muscle development in Hu sheep, which provided new insights into a better understanding of the epigenetic regulation of sheep muscle development.


Mutton is a popular meat globally, owing to its low cholesterol, low fat, and high protein content. However, the slow growth rate, low slaughter rate, and low meat yield of sheep in many countries, including China, constitute an important bottleneck that must be addressed to improve the efficiency of large-scale lamb meat production.

The skeletal muscle development is closely related to meat yield and quality in animals reared for meat. The development and growth of muscle involve the proliferation, fusion, and differentiation of myoblast cells into muscle fibers [1]. These processes are affected not only by genotype, but also a set of complicated epigenetic regulatory mechanisms, including DNA methylation. At present, although the mechanism involved in muscle development have been studied at the signaling pathway, transcriptional, and translational levels [2, 3], less is known of the associated epigenetic regulatory mechanisms.

DNA methylation is an epigenetic regulatory mechanism that mediates numerous biological processes such as growth, development, and genomic imprinting [4]. Whole Genome DNA methylation changes in the skeletal muscle have been analyzed based on differentpig breeds, with the results highlighting the differentially methylated regions in the promoter are highly correlated with known obesity-related genes and novel genes, eg. FTO, ATP1B1, COL8A2 and so on [5]. Genome-wide DNA methylation profiling in skeletal muscle tissues of aging pigs showed that DNA methylation play a key role in improving proteolysis that is related to muscle function [6]. A comparative analysis of whole genome DNA methylation regulation of gene expressionat the level of transcription in muscles of Japanese Black and Chinese Red Steppes cattle identified several genes associated with DMRs that is related to muscle development [7]. These studies indicate DNA methylation play an important roles in muscle development.

However, little is known about the expression patterns and potential value of DNA methylation in skeletal muscle development of Hu sheep, a Chinese endemic species bred for its meat and skin. The number of sheep muscle fibers increases rapidly at 75–120 d of gestation, following which myofibers grow to fuse and hypertrophy after birth [8]. It is necessary to understand the dynamics of DNA methylation profiles in sheep muscle during these processes. Whole-genome bisulfite sequencing (WGBS) is the most comprehensive DNA methylation sequencing methods available, achieving single-base resolution through bisulfite conversion. WGBS have excellent specificity and non-sensitivity, and can obtain almost complete information of methylcytosine [9]. In our study, we systematically analyzed the DNA methylation profiles in sheep muscle at two key developmental stages (110-day fetus and two-year-old adult) using WGBS technology, thereby expanding the sheep muscle methylome catalog.


DNMTs expression levels

The expression levels of DNMTs (DNMT1, DNMT3A, and DNMT3B) in LD muscle of fetal and adult sheep were first analyzed by Quantitative reverse transcription-PCR (qRT-PCR). The expression levels of DNMT1, DNMT3A, and DNMT3B in the LD muscle of adult sheep were significantly lower than those in fetal LD muscle (Fig. 1) (P < 0.05).

Fig. 1
figure 1

The mRNA expression levels of DNA methyltransferases (DNMTs) determined by qRT-PCR. The relative expression of DNMTs in ovaries was detected by qRT-PCR. The experiment was performed using three biological repeats and three technical repeats. The relative expression levels were normalized to that of GAPDH. The results are expressed as means±SEM relative to the fetal samples and the ordinate represents log10-transformed values. **P < 0.01

Genome-wide DNA methylation profiling

Global DNA methylation analysis of the LD muscle was performed by WGBS with 30× genome coverage and > 99% conversion efficiency. A total of 78.17 and 75.90 Giga raw bases were generated on average for fetal and adult muscle, respectively. After filtering out low-quality data, approximately 230 million clean reads were generated for each group, with the Q30 of clean, full-length reads ranging from 90.86 to 93.01%. The mapped reads were used for subsequent analysis as the rates ranged from 69.46 to 72.21%. Details of the quality of sequencing data are shown in Table 1.

Table 1 Sequencing data by whole genome bisulfite sequencing (WGBS) for sheep Fetus and Adult stages

All methylated genomic C sites were approximately 3.5% in each group (Table 1). The CG, CHH, and CHG (where H is A, C, or T) methylation levels were different. We found genome-wide methylated cytosine (mC) levels of 88.87 ± 0.67% for CG, 2.58 ± 0.16% for CHG, and 8.55 ± 0.52% for CHH in fetal samples, and 85.33 ± 0.95% for CG, 3.31 ± 0.21% for CHG, and 11.36 ± 0.74% for CHH in adult samples, and proportions of these contexts were similar in each group (Fig. 2).

Fig. 2
figure 2

The average ratio of DNA methylation types in fetal and adult genomes of Hu sheep. The blue, orange, and gray colors represent methylated (m) CG, mCHG, and mCHH, respectively

A violin graph was drawn with points representing different levels of methylation. The CG methylation levels were high with wide sections in the violin graph (Fig. 3a), but CHG and CHH methylation levels were low with narrow sections in the violin graph (Fig. 3b and c). And then chromosome methylation maps for fetus and adult samples were plotted. The results showed that most chromosomal cytosine hypermethylation was found in the CG context and that the chromosomal mC sites were different between the fetal and adult stages (Additional file 9).

Fig. 3
figure 3

Violin plot for the overall distribution of methylation levels for different methylation types. (a) CG, (b) CHG, and (c) CHH. Fetus (fetus 1, fetus 2, fetus 3), adult (adult 1, adult 2, adult 3). H = A, C or T. The abscissa represents the different samples, the ordinate represents the level of methylation of the samples; the width of each violin represents the density of the point at that methylation level, while the boxplot shows the methylation levels in each violin

To further compare the genome-wide distribution and the methylation levels of various functional genomic elements between the two developmental stages, we analyzed the methylation status of six different regions, including the promoters, 5′UTRs (untranslated regions), exons, introns, 3′UTRs and distal intergenic. No significant differences were observed among the different genetic elements for the three mC contexts. Overall, the methylation levels in the CG context were higher than those in the CHG and CHH contexts, where the CHH context was hypomethylated and stable in all the functional elements and the CHG context was almost entirely unmethylated. The DNA methylation levels in the CG context were highest in introns, followed by exons (except the first exon) and downstream regions, with sites near the transcription start site (TSS) showing the lowest level. The methylation levels gradually decreased from the promoters to the TSSs and increased from the TSSs to the introns. More detailed information is listed in Fig. 4 and Additional file 3.

Fig. 4
figure 4

DNA methylation levels across genomic elements in fetal and adult sheep. (a) Adult, (b) fetus. The abscissa represents different genomic elements, with a, b, c, d, e, f, and g denoting upstream, first exon, first intron, inner exon, inner intron, last exon, and downstream, respectively. The left ordinate represents the methylation levels of CG/CHG contexts, and the right ordinate represents the methylation levels of the CHH context. The dotted, green, vertical line represents the transcription start site (TSS), and the red, orange, and blue solid lines represent CG, CHH, and CHG, respectively

Characterization of DMRs

We had identified 48,491 differentially methylated CG regions, 17 differentially methylated CHG regions, and 135 differentially methylated CHH regions. Among the DMRs, 21,640 (CG:21528 + CHG:10 + CHH:112) were hypermethylated and 26,937 (CG:26943 + CHG:7 + CHH:23) hypomethylated. The DMRs were mostly located at distal intergenic regions, followed by introns, exons, and regulatory regions such as promoters, 5′UTRs, and 3′UTR. In the CG context, only 41,151, and 1250 DMRs were in 5′UTRs, 3′UTRs, and promoters, respectively. More detailed information is listed in Fig. 5.

Fig. 5
figure 5

Identification of differentially methylated regions (DMRs) among the fetal and adult sheep muscle samples. (a) The number of differentially methylated regions between different methylation types. Histograms show the distribution numbers of DMRs in different genomic elements in the CG (b), CHG (c), and CHH (d) contexts

To detect relatedness between samples individuals, pearson’s correlation coefficient (r2) was used as the evaluation index. The results showed the correlation between replicates in each group was high, thus indicating that further data analysis is reasonable (Additional file 10). In addition, as shown in the heat maps in Fig. 6, we analyzed genome-wide methylation in sheep at the fetal and adult stages using hierarchical clustering. The results showed a clear separation between the two developmental stages. More detailed DMR results are listed in Additional file 4.

Fig. 6
figure 6

Heat map cluster analysis of differentially methylated regions (DMRs) in different gene functional regions. In the heat map, highly methylated loci are displayed in red and sparsely methylated loci in blue. In addition, the red, yellow, green, turquoise, blue, purple, and pink colors indicate upstream, first exon, first intron, inner exon, inner intron, last exon, and downstream, respectively, and are shown above the heatmap

Validation of WGBS data by bisulfite sequencing

To verify the reliability of the WGBS-Seq data, four regions were randomly selected for bisulfite sequencing PCR (BSP). Although the differences in methylation levels among the DMRs (Table 2 and Fig. 7) validated by BSP were lower than those obtained by WGBS, trends were consistent, and the differences might be due to differences in methylation levels of different animals in each stage. On the whole, the BSP results agreed well with the WGBS data, indicating that the WGBS data were reliable and suitable for further study.

Table 2 DMR methylation levels of DLK1, KLHL31, FADS2 and RTL1 during muscle development of Hu sheep
Fig. 7
figure 7

Validation of the whole-genome bisulfite sequencing (WGBS) data by bisulfite sequencing PCR. *P < 0.05, **P < 0.01

GO and KEGG enrichment analysis of DMGs

To explore the change in the methylation status of genes under muscle development, the GO and KEGG databases were used to annotate 11,522 DMGs detected in the DMRs. Because most of the DMGs were of the CG context (more than 95%), we focused on CG methylation for the DMG functional enrichment analysis. Based on the GO database, the terms that play a key role in muscle growth and are significantly enriched (corrected P < 0.05), including embryonic skeletal system development, skeletal muscle cell differentiation, and skeletal muscle tissue development. The 20 most significantly differentially enriched muscle development-related GO terms for DMGs between fetal and adult LD samples are shown in Fig. 8a. According to the KEGG pathway analysis, DMGs were significantly enriched in the Hippo, cAMP, PI3K-Akt, calcium, and MAPK signaling pathways. The 21 muscle development-related KEGG terms with a corrected P-value < 0.05 for the DMGs are listed in Fig. 8b. The results suggest that these DMGs, which are influenced by DNA methylation, can affect muscle development. More detailed results of the COG (cluster of orthologous groups of proteins), GO, and KEGG analyses of CG, CHG, and CHH methylation are shown in Additional file 11-13 and Additional file 5.

Fig. 8
figure 8

Enrichment analysis of CG-type differentially methylated genes. (a) Counts of DMGs enriched for the top 20 muscle development-related Gene Ontology (GO) terms. The abscissa represents the number of DMGs, and the ordinate shows the GO pathway terms. (b) The scatterplot of 21 muscle development-related Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The ordinate represents the enriched pathways, and the abscissa represents the Rich factor of the corresponding pathways; the size of the spots represents the number of genes related to DMRs enriched in each pathway, while the color of the spot represents the corrected P-value for each pathway. The Rich factors indicate the ratio of the number of DMGs mapped to a certain pathway to the total number of genes mapped to this pathway. Greater Rich factor means greater enrichment. DMRs: differentially methylated regions; DMGs: differentially methylated genes

Screening of potential functional differential methylation genes involved in muscle development

To identify key genes involved in the regulation of skeletal muscle development, we set three limiting factors to perform an association analysis. First, we screened 1914 candidate genes known to be associated with muscle development through GO and KEGG functional enrichment analysis of significant DMGs (Additional file 6-7). Second, using our previous RNA-seq data between fetus and adult, we further screened 525 overlap genes between DMGs and differentially expressed genes (DEGs). Third, we identified 159 genes with fragments per kilobase of exon model per million reads mapped (FPKM) values > 10 and a fold change > 4, and their interaction network was generated using STRING software. In total, we identified 118 candidate genes that interacted with each other. As shown in Fig. 9, ADIPOQ, CCNA2, ITGA2, MYOG, MAPT, DIAPH1, NR4A1, DLK1, and COL1A2 were identified as hub genes in the interaction network related to the muscle development pathway. More detailed results on the abovementioned genes are listed in Additional file 8.

Fig. 9
figure 9

Construction of the network of differentially methylated genes (DMGs) related to muscle development. Analysis of the interaction between DMGs related to muscle development using STRING software according to the interplay index (confidence > 0.7). The interplay index between genes was represented by edge width and transparency. Dark and wide edges indicated high confidence

The regulatory effect of differential gene methylation on the development of sheep muscle

To investigate the effect of DNA methylation on gene expression levels, we compared the trend between gene expression and methylation levels using the FPKM value for the RNA-seq data and the difference in methylation levels between fetal and adult WGBS-seq data samples. The results showed that the DMRs of ADIPOQ and ITGA1 genes were in distal intergenic and intron respectively, trends of DNA methylation levels of those genes were consistent with those of expression. Furthermore, qRT-PCR results showed that the expression level of the ADIPOQ gene was upregulated with development, while that of ITGA1 was downregulated. In addition, DIAPH1 and NR4A1 genes have DMRs in intron and distal intergenic besides DMRs in promoter regions, the DNA methylation levels of the DIAPH1 and NR4A1 genes were downregulated at the adult stage, which was the opposite of that observed for their expression levels. Furthermore, the qPCR results showed that the expression levels of the DIAPH1 and NR4A1 genes were significantly upregulated at the adult stage. DLK1 genes have more DMRs in exon and distal intergenic besides 1 DMR in promoter regions, and CCNA2 have 1 DMR in intron regions. The DNA methylation levels of the DLK1 and CCNA2 genes were upregulated at the adult stage, which was the opposite of that observed for their expression levels. Furthermore, the qPCR results showed that the expression levels of DLK1 and CCNA2 were downregulated at the adult stage. As MAPT, MYOG and COL1A2 genes may have contained more DMRs in gene body (exon, intron and distal intergenic) and promoter, the DNA methylation levels were inconsistent. And the expression of the MAPT gene tended to be upregulated, while that of the MYOG and COL1A2 genes tended to be downregulated at the adult stage. The qPCR results were in good agreement with the RNA-seq data. However, the levels of DNA methylation in the promoter regions of MAPT, DIAPH1, NR4A1, and DLK1 were the opposite of that observed for their expression levels. The results indicated that DNA methylation in promoter regions of MAPT, DIAPH1, NR4A1, and DLK1 affected their gene expression levels during skeletal muscle development. While the effect of DNA methylation in gene body regions on ADIPOQ, DIAPH1, CCNA2, ITGA1 and COL1A2 genes expression was variable. More detailed information is listed in Table 3 and Fig. 10.

Table 3 DMGs most likely involved in muscle development
Fig. 10
figure 10

qRT -PCR was performed to detect the relative mRNA expressions of candidate genes in longissimus dorsi (LD) muscle. GAPDH was used as internal control. Values are expressed as means ± SEM of three replicates. *P < 0.05, **P < 0.01


DNA methylation is a epigenetic regulation form with important roles in gene expression and tissue development [10]. Although muscle DNA methylation has been analyzed in cattle [11], pigs [12], humans [13], mice [10], and sheep [14], genome-wide DNA methylation analysis of sheep muscle has only been performed for specific types of metabolism, in different breeds, at one developmental stage. To our knowledge, this is the first systematic comparison of genome-wide DNA methylation profiles of the LD muscle between fetal and adult Hu sheep.. In our previous study, we confirmed through transcriptome analysis that numerous genes and pathways related to growth and development were differentially expressed between the fetal and adult stages in Hu sheep [15]. Therefore, the genome-wide DNA methylation profile of sheep muscle were investigated by WGBS to elucidate the relationship between differential muscle development and DNA methylation.

DNA methylation is mainly caused by the DNA methyltransferase (DNMT) family, including DNMT1, DNMT3A, and DNMT3B. DNMT1 mainly plays a role in maintaining DNA methylation, while DNMT3A and DNMT3B catalytic de novo methylation. In our study, we found that the expression levels of DNMT1 and DNMT3A in the fetus were higher than in the adult, indicating that DNMTs may play an important role in regulating the transcription of genes related to sheep muscle development.

The muscle of genome-wide methylation patterns were similar between fetus and adult in functional genomic regions. However, there were differences among the three mC contexts which might be related to differences in the sequences of the different genetic elements [7]. After filtering the WGBS-seq raw data, we obtained 230 million clean reads per sample. The rate of uniquely mapped reads ranged from 69.46 to 72.21%, higher than the values found in skeletal muscle satellite cells of mice by MeDIP [10]. Approximately 3.5% of cytosine sites were methylated, with the highest proportion of CG methylation. These results were similar to those found in other species [16] and tissue [17]. Among the gene functional regions, TSSs presented the lowest methylation levels, which was consistent with the results found in sheep ovaries by WGBS [17].

In our study, 48,643 DMRs and 11,522 genes related to these DMRs were identified. Most of the DMRs were small fragments (50–1000 nucleotides; > 90% of the DMRs), suggesting that methylation changes in small regions might be involved in the regulating gene expression. The DMRs were only a small proportion in promoter regions, 3′ UTRs, and 5′UTR, and mainly concentrated in distal intergenic regions (> 60%) with s. To further validate the sequencing results, we randomly selected four regions by BSP, and found that the BSP results were consistent with the sequencing data.

We tried to reveal the roles of DMGs by functional annotations. GO analysis of the DMGs showed that some terms that play a key role in muscle growth, such as skeletal muscle tissue development, embryonic skeletal system development, skeletal muscle cell differentiation, and myoblast fusion. Furthermore, KEGG analysis confirmed that several genes differentially methylated between the two muscle developmental stages might be related to muscle growth. Including Hippo, cAMP, PI3K-Akt, calcium, and MAPK signaling pathways. The PI3K-Akt signaling pathway is critical for skeletal muscle protein synthesis and degradation [18]. The calcium signaling pathway is the key pathway exerting allosteric regulation on many proteins, including through ion channel activation or by acting as a secondary messenger, which could directly affect skeletal muscle metabolism [7]. As muscle is a major metabolic tissue, it was not unexpected that metabolic pathways, including carbon metabolism, insulin secretion, and insulin signaling, were enriched in our study. The bioinformatic results showed that the DMGs related to these regulatory processes show significant differences between fetus and adult, indicating that they may be important for myofiber growth and muscle metabolism. However, how DNA methylation affects gene expression and how these genes work together are still poorly understood.

We generated DMG interaction networks to determine whether DMGs play a determinative role in sheep muscle function. Network analysis showed that ADIPOQ, CCNA2, ITGA2, MYOG, MAPT, DIAPH1, NR4A1, DLK1, and COL1A2 were the key nodes.

DIAPH1 participates in the regulation of cell division, transcriptional activity of serum response factor, and cell motility [19]. DIAPH1 also play a role in signal transduction in smooth muscle cells [20], and DIAPH1 blockade could reduce cardiac muscle cell damage after myocardial infarction [21]. Dlk1 encodes a transmembrane protein [22], and has been shown to inhibit myoblast proliferation and enhance differentiation when overexpressed in cell culture [23]. The level of DLK1 mRNA was upregulated in LD in callipyge lambs [24]. Methylation of the DLK1 gene promoter region increased the invasive ability of non-small cell lung cancer cells [25]. Cyclin A2 (CCNA2) is ubiquitously expressed, and plays an important role in controlling progression through mitosis [26]. Cardiomyocyte mitosis and number were increased in CCNA2-treated pigs [27]. Recent studies have shown that the ITGA1 gene encodes the integrin receptor alpha 1 subunit and negatively regulates cell proliferation [28]. Moreover, ITGA1 is important for late muscle differentiation and proliferation in the pig [29]. However, little is known about how DNA methylation of DIAPH1, DLK1, ANGPTL4, CCNA2, and ITGA1 regulates mammalian muscle growth. Several studies have shown that dynamic changes in DNA methylation patterns persist during development and cell differentiation [30, 31]. This suggests that these genes might play important roles in muscle proliferation, and differentiation and regulation of these genes through DNA methylation might be one of the mechanisms determining differential muscle development between the fetus and adult in sheep.

Many of the genes differentially methylated between fetal and adult sheep identified in the present study are also involved in insulin secretion and lipid and glucose metabolism. ADIPOQ (encoding adiponectin) is an important regulator of lipid metabolism and glucose and is exclusively secreted from adipose tissue [32]. Short-term fasting can cause significant changes in DNA methylation in the ADIPOQ gene promoter in adipose tissue [33]. Moreover, ADIPOQ methylation levels were increased by saturated fatty acid overfeeding [34]. In our study, only one DMR (chr1: 198577572–198,577,648, strong hypomethylation in adult sheep) was related to ADIPOQ and located in the distal intergenic region, indicating that ADIPOQ expression may be influenced by this DMR in the same manner that intragenic DNA methylation can downregulate IGF2 gene expression in cattle [35]. Importantly, the mRNA expression of ADIPOQ may determine the responses of follicular to gonadotropins, thereby inducing ovum release. Mutations in the MAPT gene are associated with amyotrophic lateral sclerosis [36], while skeletal muscle autophagy can be regulated through MAPT [37]. NR4A1, a transcription factor, was shown to regulate the expression of genes involved in wasting the mitochondrial energetic budget and activating the AMPK catabolic pathway [38]. Myofiber size and muscle mass decrease in mice knockout Nr4a1 [39]. In contrast, NR4A1 can promote the expression of genes that control glucose uptake and glycolysis in skeletal muscle [40]. Previous studies have shown that DNA methylation in promoter region were negatively correlated with the NR4A1gene expression level [41]. In addition, the COL1A2 gene encodes the pro-alpha2 chain of type I collagen, a protein found in most connective tissues. Mutations in COL1A2 are associated with myopathy [42]. Combined, these results and those from our study indicate that DNA methylation can affect muscle development in sheep.

A complex relationship exists between gene DNA methylation and gene expression levels [7]. Although DNA methylation in promoter regions can inhibit gene expression [43], how DNA methylation within the gene body affects gene expression is poorly understood [44]. In our study, hypermethylation of promoter regions inhibited DLK1 gene expression, while hypomethylation of promoter regions induced the expression of MAPT, DIAPH1, and NR4A1, consistent with previously reported results [45]. Therefore, we concluded that MAPT, DIAPH1, NR4A1, and DLK1 were the key regulatory genes during skeletal muscle development, and their DNA methylation status may be the key functional regulators of muscle development. Methylation of these genes may partially contribute to the significant variation observed in muscle development and lipid metabolism between fetal and adult sheep. Nonetheless, the epigenetic mechanisms involved in the regulation of these genes and genetic regions involved in muscle development and lipid metabolism require further study.


In conclusion, we have provided the first systematic description of genome-wide DNA methylation patterns of sheep muscle at the fetal and adult stages. We investigated several novel and important DMRs/DMGs and pathways related to muscle development in sheep. The results provide valuable data for further understanding the genetic and epigenetic mechanisms that control economic traits in sheep, which could be used to mark assisted selection procedures to promote the growth of skeletal muscle of sheep.


Animals and sample collection

The six ram at fetus (110 days old, 1.36 ± 0.14 kg) and adult (2 years old, 77.98 ± 3.19 kg) stages (n = 3) were supplied from Taizhou Hailun Sheep Industry Co., Ltd. (Taizhou, China). The sheep were raised under the same conditions, with natural light and free access to food and water. All animals were fasted overnight and were then euthanized by captive bolt stunning and exsanguination. The LD muscle samples were collected from between the 12th and 13th thoracic vertebrae of the right side at the fetus and adult stages, immediately frozen in liquid nitrogen, and stored at − 80 °C until use.

Library preparation

DNA was isolated from LD muscle samples using a DNA extraction kit (Tiangen, Beijing, China). The DNA concentration and quality were determined by NanoDrop (NanoDrop Technologies, Wilmington, DE, USA) and agarose gel electrophoresis. Three DNA libraries were constructed for each group. Equal amounts of genomic DNA (2 μg per sample) were fragmented to 400–500 bp by ultrasonication, followed by adenylation and end-repair. The selected fragments were treated with bisulfite and then amplified by PCR to generate the sequencing libraries.

WGBS and identification of DMRs

The library was sequenced using an IlluminaHiSeqTM2500 platform (Biomarker Technologies, Beijing, China). The peak signal was transformed into sequence data by base calling, following which the raw reads were quality-filtered to obtain the clean reads. First, reads were trimmed of the 3′ adapter sequence. Then, reads with > 10% unknown bases (N) and those of low quality (more than 50% of bases with a PHRED score ≤ 5) were removed. We also calculated the Q30 and GC content.

The clean reads were aligned to the sheep reference genome (Oar_v3.1) and the bisulfite mapping of methylation sites was performed using Bismark software. The duplicates were reads that aligned with the same region of the genome, and can estimated the sequencing depth and coverage. The bisulfite conversion rate is the percentage of methylated clean reads to the total number of clean reads in the genome. The binomial distribution test for each C site was used to confirm C-site methylation by screening conditions for coverage ≥4× and false discovery rate (FDR) < 0.05.

To identify the differentially methylated regions (DMRs) between fetal and adult samples, we referenced the model of [16] to estimate the methylation level. All C sites with read coverage > 10× were used for DMR analysis with MOABS [46]. Subsequently, DMRs were defined by the presence of at least three methylation sites in the region, and in which the difference in methylation levels was > 0.2 (> 0.3 for the CG context) and the P-value from Fisher’s exact test was < 0.05.

Functional enrichment analysis

The DMR-related genes (DMGs) were compared against functional databases such as GO and KEGG by BLAST for annotation of gene function. GO enrichment analysis of the DMGs was implemented by the GOseq R packages based on the Wallenius non-central hypergeometric distribution [47]. KOBAS software was used to test the significance of DMR-related gene enrichment in the KEGG pathway analysis [48]. Pathways with a corrected P-value < 0.05 were considered to be significantly enriched. The STRING database was used to analyze interaction networks of selected DMGs ( [49].

Quantitative reverse transcription-PCR

The expression levels of DNA methyltransferase-related genes and validate the DMGs by qRT-PCR. Total RNA was isolated from LD muscles using Trizol reagent (Invitrogen, Carlsbad, USA). cDNA was reverse transcribed from total RNA using the PrimeScript RT kit (Takara, Dalian, China). qPCR was performed on a StepOnePlus Real-Time PCR System (Life Technologies, USA) using SYBR Green Master Mix (Roche Applied Science, Mannheim, Germany). The gene primers are listed in Additional file 1. The relative expression of each gene was normalized to that of GAPDH using the 2−ΔΔCt method [50].

Bisulfite sequencing PCR

The bisulfite sequencing PCR was used to validate DNA methylation levels of selected candidate genes. Genomic DNA was modified with sodium bisulfite using the EZ DNA Methylation-Gold™ Kit (ZymoResearch, Los Angeles, USA). Then, bisulfite-converted gDNA was subjected to PCR amplification using Zymo Taq™ DNA polymerase (ZymoResearch). The PCR products were purified using a Gel Extraction Kit (Shenggong, Shanghai, China), ligated, and cloned into the pUC18-T vector (Shenggong). Fifteen clones of each sample were randomly selected for DNA sequencing. The quantification tool for methylation analysis was used to analyze bisulfite sequencing data (QUMA; Gene sequence-specific primers are listed in Additional file 2.

Association analysis

We previously screened many genes related to muscle development at two stages of Hu sheep (fetus and adult) using the Illumina platform, [15]. By association analysis of the differentially methylated genes and the differentially expressed genes, a set of differentially methylation DEGs at the intersection of the two was obtained. Negative correlations between DMR methylation level and the corresponding gene expression level were identified by correlation analysist (r with a negative value).

Statistical analysis

Statistical analyses were performed by the independent samples t-test with the SPSS 25.0 software package (SPSS Inc., Chicago, IL, USA). Results of the qRT-PCR data were expressed as means ± standard error of the mean (SEM) of three samples with three biological replicates. Differences were regarded as significant at P < 0.05.

Availability of data and materials

The full WGBS-data sets have been submitted to NCBI BioProject under Accession: PRJNA622418. The original data of RNA sequencing has been uploaded to GEO database, the number is GSE101669. Additional data can be found in Additional files.



Bisulfite sequencing PCR


Cluster of orthologous groups of proteins


Differentially expressed genes


Differential methylated regions related genes


Differential methylated regions


DNA methyltransferase


Fragments per kilobase of exon model per million reads mapped


Longissimus dorsi


False discovery rate


Gene Ontology


Glyceraldehyde-3-phosphate dehydrogenase


Kyoto encyclopedia of genes and genomes


Methylated cytosine


Quantitative reverse transcriptase PCR


Transcription Start Site


Untranslated regions


Whole-genome bisulfite sequencing


  1. Parker MH, Seale P, Rudnicki MA. Looking back to the embryo: defining transcriptional networks in adult myogenesis. Nat Rev Genet. 2003;4(7):497–507.

    Article  CAS  PubMed  Google Scholar 

  2. Lee DE, Bareja A, Bartlett DB, White JP. Autophagy as a Therapeutic Target to Enhance Aged Muscle Regeneration. Cells. 2019;8(2).

  3. Song C, Yang J, Jiang R, Yang Z, Li H, Huang Y, Lan X, Lei C, Ma Y, Qi X, Chen H. miR-148a-3p regulates proliferation and apoptosis of bovine muscle cells by targeting KLF6. J Cell Physiol. 2019.

  4. Smith ZD, Meissner A. DNA methylation: roles in mammalian development. Nat Rev Genet. 2013;14(3):204–20.

    Article  CAS  PubMed  Google Scholar 

  5. Li M, Wu H, Luo Z, Xia Y, Guan J, Wang T, Gu Y, Chen L, Zhang K, Ma J, Liu Y, Zhong Z, Nie J, Zhou S, Mu Z, Wang X, Qu J, Jing L, Wang H, Huang S, Yi N, Wang Z, Xi D, Wang J, Yin G, Wang L, Li N, Jiang Z, Lang Q, Xiao H, Jiang A, Zhu L, Jiang Y, Tang G, Mai M, Shuai S, Li N, Li K, Wang J, Zhang X, Li Y, Chen H, Gao X, Plastow GS, Beck S, Yang H, Wang J, Wang J, Li X, Li R. An atlas of DNA methylomes in porcine adipose and muscle tissues. Nat Commun. 2012;3.

  6. Zykovich A, Hubbard A, Flynn JM, Tarnopolsky M, Fraga MF, Kerksick C, Ogborn D, MacNeil L, Mooney SD, Melov S. Genome-wide DNA methylation changes with age in disease-free human skeletal muscle. Aging Cell. 2014;13(2):360–6.

    Article  CAS  PubMed  Google Scholar 

  7. Fang X, Zhao Z, Yu H, Li G, Jiang P, Yang Y, Yang R, Yu X. Comparative genome-wide methylation analysis of longissimus dorsi muscles between Japanese black (Wagyu) and Chinese Red Steppes cattle. Plos One. 2017, 12(8).

  8. Wei C, Wu M, Wang C, Liu R, Zhao H, Yang L, Liu J, Wang Y, Zhang S, Yuan Z, Liu Z, Hu S, Chu M, Wang X, Du L. Long noncoding RNA Lnc-SEMT modulates IGF2 expression by sponging miR-125b to promote sheep muscle development and growth. Cell Physiol Biochem. 2018;49(2):447–62.

    Article  CAS  PubMed  Google Scholar 

  9. Kurdyukov S, Bullock M. DNA Methylation Analysis: Choosing the Right Method. Biology (Basel). 2016;5(1).

  10. Zhang W, Zhang S, Xu Y, Ma Y, Zhang D, Li X, Zhao S. The DNA methylation status of Wnt and Tgfbeta signals is a key factor on functional regulation of skeletal muscle satellite cell development. Front Genet. 2019;10:220.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Huang Y, Sun J, Zhang L, Li C, Womack JE, Li Z, Lan X, Lei C, Zhang C, Zhao X, Chen H. Genome-wide DNA methylation profiles and their relationships with mRNA and the microRNA Transcriptome in bovine muscle tissue (Bos taurine). Sci Rep-Uk. 2014;4:1–17.

    Google Scholar 

  12. Ponsuksili S, Trakooljul N, Basavaraj S, Hadlich F, Murani E, Wimmers K. Epigenome-wide skeletal muscle DNA methylation profiles at the background of distinct metabolic types and ryanodine receptor variation in pigs. BMC Genomics. 2019;20.

  13. Gensous N, Bacalini MG, Franceschi C, Meskers CGM, Maier AB, Garagnani P. Age-related DNA methylation changes: potential impact on skeletal muscle aging in humans. Front Physiol. 2019;10:996.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Namous H, Penagaricano F, Del Corvo M, Capra E, Thomas DL, Stella A, Williams JL, Marsan PA, Khatib H. Integrative analysis of methylomic and transcriptomic data in fetal sheep muscle tissues in response to maternal diet during pregnancy. BMC Genomics. 2018;19(1):123.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Ren C, Deng M, Fan Y, Yang H, Zhang G, Feng X, Li F, Wang D, Wang F, Zhang Y. Genome-Wide Analysis Reveals Extensive Changes in LncRNAs during Skeletal Muscle Development in Hu Sheep. Genes (Basel). 2017;8(8).

  16. Lister R, Pelizzola M, Kida YS, Hawkins RD, Nery JR, Hon G, Antosiewicz-Bourget J, O'Malley R, Castanon R, Klugman S, Downes M, Yu R, Stewart R, Ren B, Thomson JA, Evans RM, Ecker JR. Hotspots of aberrant epigenomic reprogramming in human induced pluripotent stem cells. Nature. 2011;471(7336):68–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Zhang Y, Li F, Feng X, Yang H, Zhu A, Pang J, Han L, Zhang T, Yao X, Wang F. Genome-wide analysis of DNA methylation profiles on sheep ovaries associated with prolificacy using whole-genome bisulfite sequencing. BMC Genomics. 2017;18.

  18. Zhang Z, Li J, Liu J, Guo B, Leung A, Zhang G, Zhang B. Icaritin requires phosphatidylinositol 3 kinase (PI3K)/Akt signaling to counteract skeletal muscle atrophy following mechanical unloading. Sci Rep-Uk. 2016;6.

  19. Narumiya S, Tanji M, Ishizaki T. Rho signaling, ROCK and mDia1, in transformation, metastasis and invasion. Cancer Metastasis Rev. 2009;28(1–2):65–76.

    Article  CAS  PubMed  Google Scholar 

  20. Toure F, Zahm JM, Garnotel R, Lambert E, Bonnet N, Schmidt AM, Vitry F, Chanard J, Gillery P, Rieu P. Receptor for advanced glycation end-products (RAGE) modulates neutrophil adhesion and migration on glycoxidated extracellular matrix. Biochem J. 2008;416(2):255–61.

    Article  CAS  PubMed  Google Scholar 

  21. O'Shea KM, Ananthakrishnan R, Li Q, Quadri N, Thiagarajan D, Sreejit G, Wang L, Zirpoli H, Aranda JF, Alberts AS, Schmidt AM, Ramasamy R. The Formin, DIAPH1, is a key modulator of myocardial ischemia/reperfusion injury. Ebiomedicine. 2017;26:165–74.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Yu H, Waddell JN, Kuang S, Tellam RL, Cockett NE, Bidwell CA. Identification of genes directly responding to DLK1 signaling in Callipyge sheep. BMC Genomics. 2018;19.

  23. Waddell JN, Zhang P, Wen Y, Gupta SK, Yevtodiyenko A, Schmidt JV, Bidwell CA, Kumar A, Kuang SH. Dlk1 Is Necessary for Proper Skeletal Muscle Development and Regeneration. Plos One. 2010; 5(11).

  24. Fleming-Waddell JN, Wilson LM, Olbricht GR, Vuocolo T, Byrne K, Craig BA, Tellam RL, Cockett NE, Bidwell CA. Analysis of gene expression during the onset of muscle hypertrophy in callipyge lambs. Anim Genet. 2007;38(1):28–36.

    Article  CAS  PubMed  Google Scholar 

  25. Zhong Z, Ye Y, Guo W, He Y, Hu W. Relationship between DLK1 gene promoter region DNA methylation and non-small cell lung cancer biological behavior. Oncol Lett. 2017;13(6):4123–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Shekhar R, Priyanka P, Kumar P, Ghosh T, Khan MM, Nagarajan P, Saxena S. The microRNAs miR-449a and miR-424 suppress osteosarcoma by targeting cyclin A2 expression. J Biol Chem. 2019;294(12):4381–400.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Shapiro SD, Ranjan AK, Kawase Y, Cheng RK, Kara RJ, Bhattacharya R, Guzman-Martinez G, Sanz J, Garcia MJ, Chaudhry HW. Cyclin A2 Induces Cardiac Regeneration After Myocardial Infarction Through Cytokinesis of Adult Cardiomyocytes. Sci Transl Med. 2014; 6(224).

  28. Huang X, Huang T, Deng W, Yan G, Qiu H, Huang Y, Ke S, Hou Y, Zhang Y, Zhang Z, Fang S, Zhou L, Yang B, Ren J, Ai H, Huang L. Genome-wide association studies identify susceptibility loci affecting respiratory disease in Chinese Erhualian pigs under natural conditions. Anim Genet. 2017;48(1):30–7.

    Article  CAS  PubMed  Google Scholar 

  29. Zhao X, Mo D, Li A, Gong W, Xiao S, Zhang Y, Qin L, Niu Y, Guo Y, Liu X, Cong P, He Z, Wang C, Li J, Chen Y. Comparative Analyses by Sequencing of Transcriptomes during Skeletal Muscle Development between Pig Breeds Differing in Muscle Growth Rate and Fatness. Plos One. 2011;6(5).

  30. Stelzer Y, Wu H, Song Y, Shivalila CS, Markoulaki S, Jaenisch R. Parent-of-origin DNA methylation dynamics during mouse development. Cell Rep. 2016;16(12):3167–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Li J, Wu X, Zhou Y, Lee M, Guo L, Han W, Mo W, Cao W, Sun D, Xie R, Huang Y. Decoding the dynamic DNA methylation and hydroxymethylation landscapes in endodermal lineage intermediates during pancreatic differentiation of hESC. Nucleic Acids Res. 2018;46(6):2883–900.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Matsuzawa Y. Adiponectin. A key player in obesity related disorders. Curr Pharm Design. 2010;16(17):1896–901.

    Article  CAS  Google Scholar 

  33. Hjort L, Jorgensen SW, Gillberg L, Hall E, Brons C, Frystyk J, Vaag AA, Ling C. 36 h fasting of young men influences adipose tissue DNA methylation of LEP and ADIPOQ in a birth weight-dependent manner. Clin Epigenetics. 2017;9:40.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Perfilyev A, Dahlman I, Gillberg L, Rosqvist F, Iggman D, Volkov P, Nilsson E, Riserus U, Ling C. Impact of polyunsaturated and saturated fat overfeeding on the DNA-methylation pattern in human adipose tissue: a randomized controlled trial. Am J Clin Nutr. 2017;105(4):991–1000.

    Article  CAS  PubMed  Google Scholar 

  35. Huang Y, Zhan Z, Sun Y, Cao X, Li M, Wang J, Lan X, Lei C, Zhang C, Chen H. Intragenic DNA methylation status down-regulates bovine IGF2 gene expression in different developmental stages. Gene. 2014;534(2):356–61.

    Article  CAS  PubMed  Google Scholar 

  36. Origone P, Geroldi A, Lamp M, Sanguineri F, Caponnetto C, Cabona C, Gotta F, Trevisan L, Bellone E, Manganelli F, Devigili G, Mandich P. Role of MAPT in pure motor neuron disease: report of a recurrent mutation in Italian patients. Neurodegener Dis. 2018;18(5–6):310–4.

    Article  CAS  PubMed  Google Scholar 

  37. Li X, Hou Y, Wang X, Zhang Y, Meng X, Hu Y. Zhang Y. Biol Pharm Bull: To elucidate the inhibition of excessive autophagy of Rhodiola crenulata on exhaustive exercise-induced skeletal muscle injury by combined network pharmacology and molecular docking; 2019.

    Google Scholar 

  38. Weisova P, Anilkumar U, Ryan C, Concannon CG, Prehn JH, Ward MW. 'Mild mitochondrial uncoupling' induced protection against neuronal excitotoxicity requires AMPK activity. Biochim Biophys Acta. 2012;1817(5):744–53.

    Article  CAS  PubMed  Google Scholar 

  39. Tontonoz P, Cortez-Toledo O, Wroblewski K, Hong C, Lim L, Carranza R, Conneely O, Metzger D, Chao LC. The orphan nuclear receptor Nur77 is a determinant of Myofiber size and muscle mass in mice. Mol Cell Biol. 2015;35(7):1125–38.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Chao LC, Zhang Z, Pei L, Saito T, Tontonoz P, Pilch PF. Nur77 coordinately regulates expression of genes linked to glucose metabolism in skeletal muscle. Mol Endocrinol. 2007;21(9):2152–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Kasch J, Kanzleiter I, Saussenthaler S, Schurmann A, Keijer J, van Schothorst E, Klaus S, Schumann S. Insulin sensitivity linked skeletal muscle Nr4a1 DNA methylation is programmed by the maternal diet and modulated by voluntary exercise in mice. J Nutr Biochem. 2018;57:86–92.

    Article  CAS  PubMed  Google Scholar 

  42. Dalgleish R. The human type I collagen mutation database. Nucleic Acids Res. 1997;25(1):181–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Wagner JR, Busche S, Ge B, Kwan T, Pastinen T, Blanchette M. The relationship between DNA methylation, genetic and expression inter-individual variation in untransformed human fibroblasts. Genome Biology. 2014;15(2).

  44. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13(7):484–92.

    Article  CAS  PubMed  Google Scholar 

  45. Wu H, Zhang Y. Reversing DNA methylation: mechanisms, genomics, and biological functions. Cell. 2014;156(1–2):45–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Sun D, Xi Y, Rodriguez B, Park HJ, Tong P, Meong M, Goodell MA, Li W. MOABS: model based analysis of bisulfite sequencing data. Genome Biology. 2014;15(2).

  47. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology. 2010;11(2).

  48. Kanehisa M, Sato Y, Morishima K. BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and Metagenome sequences. J Mol Biol. 2016;428(4):726–31.

    Article  CAS  PubMed  Google Scholar 

  49. Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, Lin J, Minguez P, Bork P, von Mering C, Jensen LJ. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013;41(Database issue):D808–15.

    CAS  PubMed  Google Scholar 

  50. Wang L, You J, Zhong B, Ren C, Zhang Y, Meng L, Zhang G, Jia R, Ying S, Wang F. Scd1 mammary-specific vector constructed and overexpressed in goat fibroblast cells resulting in an increase of palmitoleic acid and oleic acid. Biochem Bioph Res Co. 2014;443(2):389–94.

    Article  CAS  Google Scholar 

Download references


We sincerely thank all the members in Feng Wang’s laboratory who contributed to sample collection and provided technical assistance.


This work was financially supported by the earmarked fund for Jiangsu Agricultural Industry Technology System (JATS [2019]419) and China Agricultural Research System (CARS-38). The funding bodies had no role in the study design, collection, analysis and interpretation of data, and in the writing of the manuscript.

Author information

Authors and Affiliations



YXF conducted the experiments and prepared the materials involved in this study. YXL and KPD performed the bioinformatic analysis. FW and ZZ designed and coordinated the research. YXF drafted the manuscript. FW, GMZ, and YLZ helped revise the manuscript. All authors read and approved the final version of the manuscript.

Corresponding author

Correspondence to Feng Wang.

Ethics declarations

Ethics approval

All experimental procedures including animal care and tissue sample collection were approved and carried out in accordance with the relevant guidelines set by the Ethics Committee of Nanjing Agricultural University, China (Approval ID: SYXK2011–0036).

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.

Primers for qRT-PCR.

Additional file 2.

Primers sequences of DNA methylation-related genes.

Additional file 3.

DNA methylation levels in gene functional elements in the Adult group and Fetus group.

Additional file 4.

Information of DMRs in the Fetus group and the Adult group.

Additional file 5.

Information about the COG, GO and KEGG analyses of DMGs. DMGs: Differential methylation genes.

Additional file 6.

The list of DMGs enriched for muscle development-related top 20 GO terms.

Additional file 7.

The list of DMGs enriched for muscle development-related KEGG terms.

Additional file 8.

Muscle development -related DMGs target genes in muscle between Fetus and Adult Hu sheep.

Additional file 9.

Plot of genome chromosome 5-methylcytosine map. A, CG type. B, CHG type. C, CHH type. H = A, C or T.

Additional file 10.

The correlation coefficients analysis of samples. The closer the number is to 1, the stronger the correlation.

Additional file 11.

GO and KEGG pathway analysis in CG type DMGs. A, GO analysis. B, KEGG analysis.

Additional file 12.

COG, GO and KEGG pathway analysis in CHG-type DMGs. A, COG analysis. B, GO analysis. C, top GO. D, KEGG analysis.

Additional file 13.

COG, GO and KEGG pathway analysis in CHH-type DMGs. A, COG analysis. B, GO analysis. C, top GO. D, KEGG analysis. E, top KEGG.

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 The Creative Commons Public Domain Dedication waiver ( 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

Verify currency and authenticity via CrossMark

Cite this article

Fan, Y., Liang, Y., Deng, K. et al. Analysis of DNA methylation profiles during sheep skeletal muscle development using whole-genome bisulfite sequencing. BMC Genomics 21, 327 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • DNA methylation
  • WGBS
  • Skeletal muscle
  • Development
  • Sheep