Skip to main content
Download PDF

Dynamic changes of genomic methylation profiles at different growth stages in Chinese Tan sheep

Journal of Animal Science and Biotechnology volume 12, Article number: 118 (2021) Cite this article

Abstract

Background

Tan sheep, an important local sheep breed in China, is famous for their fur quality. One-month-old Tan sheep have white, curly hair with beautiful flower spikes, commonly known as “nine bends”, which has high economic value. However, the “nine bends” characteristic gradually disappears with age; consequently, the economic value of the Tan sheep decreases. Age-related changes in DNA methylation have been reported and may be responsible for age-induced changes in gene expression. Until now, no genome-wide surveys have been conducted to identify potential DNA methylation sites involved in different sheep growth stages. In this study we investigated the dynamic changes of genome-wide DNA methylation profiles in Tan sheep using DNA from skin and deep whole-genome bisulfite sequencing, and compared the DNA methylation levels at three different growth stages: 1, 24, and 48 months old (mon1, mon24, and mon48, respectively).

Results

In this study, 11 skin samples from three growth stages (four for mon1, four for mon24, and three for mon48) were used for DNA methylation analysis and gene expression profiling. There were 52, 288 and 236 differentially methylated genes (DMGs) identified between mon1 and mon24, mon1 and mon48, and mon24 and mon48, respectively. Of the differentially methylated regions, 1.11%, 7.61%, and 7.65% were in the promoter in mon1 vs. mon24, mon24 vs. mon48, and mon1 vs. mon48, respectively. DMGs were enriched in the MAPK and WNT signaling pathways, which are related to age growth and hair follicle morphogenesis processes. There were 51 DMGs associated with age growth and curly fleece formation. Four DMGs between mon1 and mon48 (KRT71, CD44, ROR2 and ZDHHC13) were further validated by bisulfite sequencing.

Conclusions

This study revealed dynamic changes in the genomic methylation profiles of mon1, mon24, and mon48 sheep, and the percentages of methylated cytosines were 3.38%, 2.85% and 4.17%, respectively. Of the DMGs, KRT71 and CD44 were highly methylated in mon1, and ROR2 and ZDHHC13 were highly methylated in mon48. These findings provide foundational information that may be used to develop strategies for potentially retaining the lamb fur and thus improving the economic value of Tan sheep.

Background

Varied body structure and organ function is an unavoidable and complex process associated with aging, and results in increased morbidity and mortality [1]. Changes in DNA methylation profiles have been considered a cause of senescence since the Vaniushin group first reported an age-related decrease in methylcytosine levels in the somatic tissues of salmon [2]. Age-related changes in DNA methylation have been demonstrated in mammals, but the generality of this phenomenon in vertebrates remains unclear [3, 4]. In a previous study, CpG-rich areas of the genome became hypermethylated in rhesus monkeys, and there was typically decreased gene expression near transcription start sites [5]. In the study of the relationship between fleece development and age, diverse whole-genome methylation profiles have been found to characterize the two periods of hair follicle (HF) growth (anagen and telogen), which indicates that growth stages affect changes in methylation expression patterns, which in turn affects phenotype [6].

DNA methylation is an important epigenetic mechanism that is involved in regulation of various biological processes [7]. Methylation in the 5′- flanking region of a gene plays a pivotal role in regulating expression [8], and epigenetic markers silence exogenous transposons and imprint genes [9]. This mechanism may produce diverse phenotypes that are capable of responding to environmental changes during maturation [10]. A standard approach to measuring DNA methylation is bisulfite sequencing (BS-seq). BS-seq couples bisulfite conversion of DNA with next-generation sequencing to profile genome-wide DNA methylation at single base resolution [11, 12]. BS-seq was first used in Arabidopsis thaliana DNA methylome studies; since then, it has provided accurate, unbiased, and high-coverage DNA methylation landscapes in model animals such as humans, mice, and rats [9, 13,14,15]. In recent studies, methylation analysis in sheep focused on the methylation levels of both global and the specific genes. For example, utero-specific DNA methylation analysis showed that high DNA methylation in ewes could affect their offspring’s resistance to diseases [16]. In a specific genes study, highly expressed DNA methyltransferase 1 prevented the epithelial progenitor cells in the hair fleece bulb from over-proliferating. Over-proliferating drives differentiation; thus, high expression of DNA methyltransferase 1 can maintain normal hair fleece structure [17].

The Chinese Tan sheep, which is indigenous to China, is a short-tailed breed that provides an important source of high-quality fur. The lustrous curly fleece appears when Tan lambs are 1 month old and is characterized by a natural white color [18]. With aging, their fur quality decreases, which in turn affects their economic value [18, 19]. Moreover, we observed that age also impacts on fur quality of Tan sheep. Previous studies on Chinese Tan sheep mainly focused on genetic evaluation and phenotypes characteristic of fur [20]; very few studies have focused on deciphering the relationship between genome and phenotype at the molecular level, including genomic methylation levels, in Tan sheep at different ages. In previous studies, we assessed the transcriptomic and miRNA differences in skin tissues derived from Tan sheep with different degrees of curly fleece that were in different growth stages (1 and 48 months old) [7, 8]. Genome sequence analysis showed that no mutations occurred in candidate gene sequences at different stages. However, little is known about the biological process of DNA methylation regulation involved in different sheep growth stages.

In this study, we hypothesized that the dynamic changes of DNA methylation modulate the epigenetic changes of different growth stages in Tan sheep. Therefore, the objectives of this study were to determine differentially methylated regions (DMRs) in skin tissues from 1-, 24-, and 48-month-old Tan sheep and to evaluate the correlation between transcriptomic profiles and DNA methylation signatures. These results will provide new insights into the biological mechanisms underlying variation of the growth process and by providing foundational information that could be used to develop strategies for potentially retaining the lamb coat and thus improve the economic value of Chinese Tan sheep.

Materials and methods

Animals and sample preparation

Eleven female Chinese Tan sheep were used in this study. We sampled skin tissues from the shoulder of differently aged Tan sheep (Fig. 1) for DNA. Animals were divided into three age-based groups that represented different growth stages: 1 month (n = 4; mon1 group), 24 months (n = 4; mon24 group), and 48 months (n = 3; mon48 group). Animals were maintained under identical feeding conditions to minimize external factors. When they were 1, 24, and 48 months old, we cleaned their scapular region, from which the skin tissues were collected using sterilized scalpel blades; samples were immediately frozen in liquid nitrogen and stored at − 80 °C. All resulting wounds were treated with Yunnan Baiyao powder (Yunnan Baiyao Group Co., Ltd., China) to stop the bleeding. By following the manufacturer’s protocol, genomic DNA was isolated using a TIANamp Genomic DNA Kit (TIANGEN, Beijing, China). DNA integrity and quality were evaluated by agarose gel electrophoresis and spectrophotometry with a Nano-Drop spectrophotometer.

Fig. 1
figure 1

Wool phenotypes in Chinese Tan sheep at three growth stages (mon1, mon24, and mon48 groups)

Full size image

DNA bisulfite treatment and sequencing

After DNA samples were extracted and assessed for quality, bisulfite conversion was performed using the EZ DNA Methylation Direct Kit (Zymo Research Corporation, Irvine, CA, United States). A Covaris S220 sonicator (Covaris, Woburn, MA, United States) was used to fragment genomic DNA to a mean size of approximately 150 bp. After fragmentation, DNA fragments with recessed 3′-ends were repaired using dA to make the ends blunt. Adapters were added according to the manufacturer’s instructions. Following the protocol described in [10], the adapter-modified DNA was then reacted with sodium bisulfate. Bisulfite-treated DNA was PCR amplified, then subjected to paired-end sequencing using the Illumina HiSeq 2500 high-throughput sequencing system (Illumina, San Diego, CA, United States). Reads were 150 nucleotides in length.

BS-seq read mapping and methyl-cytosine analysis

Reads generated by BS-seq were mapped to the sheep reference genome (Oar_v3.1) from the Livestock Genomics website (http://www.livestockgenomics.csiro.au/) using Bismark (v0.16.3) [14]. In the BS-seq procedure, unmethylated cytosines in genomic DNA were converted into thymines after bisulfite treatment and PCR amplification, but methylated cytosines remained unchanged. The bisulfite conversion rate was calculated using Bismark [14] as the percentage of methylated reads in the total number of clean reads. Only one mismatch was permitted in the “seed” (the high-quality end of the read, default seed length was 28 bases) during alignment. Otherwise, default parameters were used.

The “bismark_methylation_extractor” script packaged with Bismark was used to extract methylation calls. The methylation status of each cytosine with ≥ 5-fold coverage was analyzed using the binomial distribution model. At each position, we used to read depth to test whether the number of detected cytosines exceeded the number expected because of sequencing error. The P-value was adjusted by the Benjamini-Hochberg method, and methylated cytosine residues were defined as those with a false-discovery rate (FDR) < 0.05. Next, we estimated the level of differential methylation in pairwise comparisons between the three groups (mon1 vs. mon24, mon1 vs. mon48, and mon24 vs. mon48). This analysis generated coverage statistics for methylated cytosine reads and their contexts (CG, CHG and CHH, where H can be A, T, or C, respectively). Because the methylation of a single C site cannot be discriminated by Bismark, we used the binomial distribution to confirm methylation at each C by examining sites with ≥ 5-fold coverage and FDR < 0.05. The methylation level at a C site was determined using the procedure described in Cokus et al. [21].

To identify differentially methylated CpG sites that are likely to function in gene regulation, we compared mapped CG-containing reads from mon1, mon24, and mon48 sheep. Sites were only considered if they were represented by at least 10 reads in each sheep. A site was classified as differentially methylated if the mon1, mon24, and mon48 versions of the site exhibited different methylation levels. More stringent criteria were then applied to identify methylation level differences [larger than 10% (two-tailed Fisher’s exact test, P < 0.05) with an FDR < 0.00005].

DMR identification and analysis

After removing low-quality reads, all cytosine sites with ≥ 10-fold coverage were aligned to the Ensembl sheep reference genome (Oar v3.1/2.62) using methylKit (v.0.9.2) and eDMR (v0.5.1) [22]. This algorithm identifies DMRs using the bimodal normal distribution of distances between cytosines to optimize the defined region. At least one differentially methylated cytosine and three CpGs were required per DMR, and the difference in methylation levels was required to be greater than 0.2 (0.3 for the CG context) with a Fisher’s exact test P < 0.05. The methylation level in a region was calculated as described by Cokus et al. [21].

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analyses of differentially methylated genes (DMGs)

DMGs were annotated using the GO and KEGG functional databases to infer gene function. GO enrichment analysis was conducted using the Wallenius non-central hyper-geometric distribution in the GOseq R package [14]. The KEGG database was used to test differentially expressed genes for statistically significant enrichment in biochemical pathways [15]. Interaction networks for selected DMGs were analyzed using the String database (http://string-db.org/) [13].

DMR validation using BS-seq

Because of the extreme phenotypes of 1- and 48-month-old Tan sheep, this study mainly focused on differences between the mon1 and mon48 groups. DNA was extracted from the skin tissues of seven female Tan sheep (four 1-month-old lambs and three 48-month-old adults). DNA was extracted from frozen skin tissues using extraction kits (Tiangen, Beijing, China). The extracted DNA (concentration: 50 ng/μL) was treated with the EpiTect Bisulfite Kit (QIAGEN, Valencia, CA, USA) according to the manufacturer’s protocol. The DNA samples were treated with bisulfite prior to a PCR, then subjected to pyrosequencing; the primers and conditions are listed in Table 1. The 5′-end of one primer from each pair was conjugated to biotin to allow binding of one DNA strand to streptavidin beads. The 50-μL PCR mixture contained 10 μL of PCR mix buffer (QIAGEN, Valencia, CA, USA), 1 pmol biotinylated forward primer, 1 pmol reverse primer, 2 μL bisulfite-treated genomic DNA, and water.

Table 1 Bisulfite-converted DNA PCR primers
Full size table

The degree of methylation was expressed as the proportion of 5-methylated cytosines (mC). Non-CpG cytosine residues were used as built-in controls to verify bisulfite conversion. Methylation measurements were standardized by processing batch number using a mean value of 0 and standard deviation of 1. The methylated DNA assays were designed to cover the greatest possible number of CpG sites within the promoter region; necessary PCR amplification length, target sequence length, and primers that avoided CpGs were taken into account. DNA methylation levels at multiple CpG sites were measured and the mean methylation values for each gene were calculated. Significance was determined using a Student’s unpaired t-test using SAS (v9.1; SAS Institute Inc., Cary, NC, USA); a value of P < 0.05 was considered significant.

Results

Whole-genome BS-seq and analysis

To identify the dynamic changes in genomic methylation profiles at different growth stages, 324.73 Gb, 310.81 Gb, and 229.8 Gb of raw reads were generated for the mon1, mon24, and mon48 groups, respectively. After low quality reads (Ns > 10%, low-quality sites > 40%, adapter contamination, and duplication pairs) were removed from the data set, around 500 million clean reads were left for each group. In the three comparisons, 69.97%–84.48% of reads could be mapped (Table 2), and reads were mapped across the genome. The mapped reads were used for subsequent analyses.

Table 2 Summary statistics for whole-genome bisulfite sequencing data for lambs (1- and 24-month-old) and adult sheep (48-month-old) (mon1, mon24, and mon48, respectively)
Full size table

Methylation in sheep was found in three sequence contexts: CG, CHG, and CHH. These contexts were present in similar proportions in each group, with CG as the most abundant context by a wide margin; the methylation proportions of the mCG contexts were 97.69%, 97.22%, and 97.42% for mon1, mon24, and mon48, respectively. The proportions of the other mCHG and mCHH contexts ranged from 0.48%–0.52% and 1.82%–2.25%, respectively.

Although methylation in different sequence contexts can be conveniently described using the statistics above, the frequency of DNA methylation varies by site, even within a single context. The frequency at a particular site can be estimated from the fraction of mapped reads that contain a methylcytosine residue at the appropriate position. This frequency was defined as the methylcytosine methylation level. The column diagram in Fig. 2 shows the distribution of mC levels for all methylated sites detected in each group. Although some sites were frequently methylated, the vast majority were methylated at relatively low levels.

Fig. 2
figure 2

Distribution of DNA methylation levels in mon1, mon24, and mon48 groups (*P < 0.05)

Full size image

To visualize the distribution of methylated sites as a function of genome location, we used the graphical tool Circos [23] to generate large-scale methylation maps organized by chromosome (Fig. 3, Additional Fig. S1). Cytosine methylation levels in the CG (red), CHG (green), and CHH (blue) contexts are shown separately for mon1, mon24, and mon48 methylomes.

Fig. 3
figure 3

Genome-scale view of the Tan sheep shoulder skin methylome. The red, green, and blue bars represent the methylation levels in the CG, CHG and CHH context, respectively; highlighted gray bands refer to normal annotated genes

Full size image

Motif analysis of nearby mC at different contexts

To understand the base-pair distribution of mC, we analyzed the relationship between sequence context and methylation preference. The mC sites were in a hypermethylation state, the CAG was the most common sequence motif in the CHG mC sites, and similar frequencies were discovered in the CHH context for the three groups (Additional Fig. S1).

Distribution of DNA methylation sites in functional genomic regions

Initial screening yielded 175,433 differentially methylated CpG sites between mon1 and mon24, 447,869 between mon1 and mon48, and 442,383 between mon24 and mon48. As shown in Fig. 4a, methylation levels for the CG context were similar in the three groups. Exon, UTR3, and gene regions contained the majority of mC sites. CG methylation levels in the promoter region were lower than those of the other elements; however, the levels showed a decreasing trend in intron regions. Methylation levels for the CHG and CHH contexts were higher in mon48 than the other two groups (Fig. 4b and c, respectively). These methylation levels were similar in mon1 and mon24, with lower methylation levels of the exon and UTR5 regions than other elements in the three groups.

Fig. 4
figure 4

mC distribution in different functional elements. a, b, and c show the trends in mC distribution of different functional elements for CG, CHG, and CHH contexts, respectively. A, C, or T. The abscissa represents the different regions of gene functional elements (intergenic, gene, exon, intron, promoter, UTR5, and UTR3), and the ordinate represents the ratio of mC

Full size image

Identification of DMRs among tan sheep of three growth stages

The analysis yielded 72 DMRs between mon1 and mon24, 340 DMRs between mon1 and mon48, and 276 DMRs between mon24 and mon48 (Additional file 1: Table S1). Most DMRs were in distal intergenic regions. Only 6, 25, and 20 DMRs were found in promoter regions between mon1 and mon24, mon1 and mon48, and mon24 and mon48, respectively. For DMRs that were in one of the seven described locations (intergenic, gene, exon, intron, promoter, UTR5, and UTR3), 63.89%, 72.1%, and 74.41% were in introns, and 1.11%, 7.61%, and 7.65% were in the promoter in mon1 vs. mon24, mon24 vs. mon48, and mon1 vs. mon48, respectively (Fig. 5).

Fig. 5
figure 5

DMR distribution among gene functional elements

Full size image

The length distributions of hypermethylated (methylation differences greater than 10%, P < 0.05) and hypomethylated (methylation differences less than 10%, P < 0.05) DMRs of the three growth stages are shown in Fig. 6. The length distributions were very similar for all three comparisons. The results indicate that hypomethylated DMRs may be slightly shorter on average than hypermethylated DMRs. In total, the DMR length shortened with age, regardless of hypermethylation or hypomethylation status. In our data, the methylation levels of DMRs were influenced by age, and most showed high methylation in mon48. We inferred that the organismal metabolism slows and gene expression weakens with age. More detailed results are listed in Additional Table S1.

Fig. 6
figure 6

Length distribution of hypermethylated and hypomethylated DMRs (the methylation differences greater and less than 10%, respectively; P < 0.05)

Full size image

Analyses of DMR-containing genes

GO analysis revealed that the CG type DMGs were significantly enriched in hair follicle morphogenesis (GO:0031069), hair follicle development (GO:0001942), and hair cycle processes (GO:0022405). KEGG analysis showed that the CG type DMGs were significantly enriched in the WNT (map 04310), MAPK (map 04010), and TGF-beta (map 04350) signaling pathways. Importantly, some DMGs were involved in biological processes important for cell differentiation (e.g., cell morphogenesis involved in differentiation, GO:0035120), which indicated that these genes may influence aging and fleece growth. For more detailed results from the GO and KEGG analyses of DMGs, see Additional Table S2, Additional Table S3 and Additional Fig. S2.

Correlation analysis of DMGs and sheep growth

To further understand the potential relationships between DNA methylation and growth stages, we performed an association analysis using two limiting factors: DMGs were enriched in fleece growth-related functions in GO analysis, and prior reports showed that these DMGs are related to growth. In our data, 51 genes were detected that met these two criteria (Table 3 and Additional Table S1). Some DMGs were hyper-methylated in mon1 that were hypomethylated in mon24 and mon48, and vice versa. Overall, methylation level showed a significant increasing trend with age in Chinese Tan sheep (P < 0.05) (Fig.7).

Table 3 Detailed information of DMGs associated with curly fleece formation
Full size table
Fig. 7
figure 7

Methylation levels of DMGs related to curly fleece growth in mon1, mon24 and mon48. The ordinate represents the average methylation levels of DMGs, the abscissa represents the different growth stages (**P < 0.01; *P < 0.05)

Full size image

DMR validation using BS-seq

The DMRs in KRT71, CD44, ROR2, and ZDHHC13 identified by whole-genome DNA methylation analysis were selected for validation using the bisulfite Sanger sequencing method. The evaluated methylated CpG sites were all within the promoter region of these four genes. Methylation levels of the selected DMRs were consistent with the sequencing results (Fig. 8), BS-seq revealed that KRT71 and CD44 methylation levels were higher in mon48 than mon1 (P < 0.01), whereas ROR2 and ZDHHC13 exhibited higher methylation levels in mon1 than mon48 (P < 0.01). The concordance between these obtained data indicated that the DMGs identified by whole-genome BS-seq were credible.

Fig. 8
figure 8

Methylation levels of four DMR-associated genes validated by pyrosequencing analysis between mon1 and mon48 (**P < 0.01)

Full size image

Discussion

Epigenetic modifications, such as changes in global and locus-specific DNA methylation, are suspected to play an important role in aging [24]. One possible mechanism is DNA methylation, which is involved in age-dependent animal phenotypes and age-related physiological variation [25]. A recent study revealed that DNA methylation patterns established in youth, in combination with other epigenetic marks, were able to accurately predict changes in transcript trajectories with aging [26]. In this study, methylation level showed a significant difference was observed between mon1 and mon48. Different methylation trends were found among Tan sheep at three growth stages, which confirmed the conclusion that the methylation patterns changed with age. The dynamic changes in DNA methylation of Tan sheep at different growth stages may represent regulatory mechanisms of changes in sheep biological development. DMGs play a major role in transcription regulation and are associated with either gene silencing or transcription elongation, depending upon the specific location of the methylated CpG sites in the gene [27]. In mammals, some transcription factors exhibit affinity to methylated DNA oligonucleotides outside their usual consensus motifs [28, 29]. In our study, 576 DMGs, which represented 1.42% of annotated genes within the current genome assembly, were identified between mon1 and mon24, mon1 and mon48, and mon24 and mon48. The ratio of DMGs in the promoter region between mon1 and mon48 was higher than that in the other two comparisons, which indicated that more genes silenced their function with age in Tan sheep. Therefore, these DMGs may be important post-selected genes for subsequent studies to address the relationship between methylation modifications and age in sheep.

GO and KEGG enrichment analyses revealed that the genes in DMRs were significantly enriched in the WNT, AMPK and MAPK signaling pathways; interestingly, most of these pathways are associated with wool development [30,31,32,33]. The WNT/β-catenin signaling pathway in particular plays an essential role in many aspects of development, and it is required for the initiation of curly hair development [34]. A previous study showed that diverse whole-genome methylation profiles characterized the two periods of HF growth (anagen and telogen), and increased DNA methylation may suppress HF development in anagen [35, 36]. IRF2 and STAT5A are transcription factors that serve as mediators that regulate transcriptional processes in HF growth and skin disease, and are affected by the methylation levels of binding genes [37]. We found 10 DMR-containing genes between mon1 and mon48 that were enriched in the WNT pathway, including WNT10A, KRT71, and TXB3, and some were already known to be involved in hair morphology. The functions of these DMRs require further validation in the future.

WNT10A was identified as a dermal papilla signature gene in mouse pelage follicles [38]. SNP analysis revealed that common variations in WNT10A have pleiotropic effects on hair morphology [39]. Moreover, other researchers found that WNT10A may participate in catagen by inhibiting keratinocyte proliferation and functioning as a WNT inhibitor [40]. In this study, we found that WNT10A was hypomethylated in mon1 but not mon24 or mon48. This result is consistent with the idea that WNT10A negatively affects curly hair growth [41]. KRT71 is also of great importance in curly fleece growth. Several studies indicated that KRT71 has positive effects on curly fleece growth [21, 42, 43]. Most previous studies modified the effects of KRT71 by producing DNA sequence or transcription regulation variants, and our data showed that the methylation levels of KRT71 changed in the three growth stages. Our BS-seq data showed that the methylation level of the KRT71 promoter region increased with age. We speculate that methylation modification of KRT71 is one reason for wool phenotype changes in Tan sheep. A previous study showed that TBX3 was necessary in embryonic cells for specifying epithelial cell fates in organs that require epithelial-mesenchymal interactions for their development [38]. TBX3 is another gene that is involved in the WNT pathway [44]. Our data showed that the methylation levels of TBX3 differed in mon1, mon24, and mon48. However, the methylation levels showed a decreasing trend with age, which requires further investigation.

Conclusions

This study reports the dynamic changes in genomic methylation profiles at different growth stages of Chinese Tan sheep. Differential DNA methylation profiles were observed in the three growth stages, and were as well associated with differential gene expression. We identified 51 DMGs related to age and fleece growth. Importantly, among DMGs, KRT71, WNT10A, and TBX3 were identified in this study as potentially important genes associated with curly fleece formation, and previous studies demonstrated that these genes play significant roles in regulating keratinocyte proliferation and migration. The data presented here highlight the importance of epigenetic profiles in different growth stages of Tan sheep to help elucidate the variation of gene expression in different growth stages and improve the economic value in sheep.

Availability of data and materials

Additional data can be found in supplementary files.

Abbreviations

BS-Seq:

Bisulfite sequencing

DMR:

Differentially methylated region

FDR:

False discovery rate

DMG:

Differentially methylated gene

GO:

Gene ontology

KEGG:

Kyoto encyclopedia of genes and genomes

PCR:

Polymerase chain reaction

mC:

5-methylated cytosines

HF:

Hair follicle

References

  1. Steegenga WT, Boekschoten MV, Lute C, Hooiveld GJ, de Groot PJ, Morris TJ, et al. Genome-wide age-related changes in DNA methylation and gene expression in human PBMCs. Age (Dordr). 2014;36(3):9648. https://doi.org/10.1007/s11357-014-9648-x.

    Article  CAS  Google Scholar 

  2. Berdyshev GD, Korotaev GK, Boiarskikh GV, Vaniushin BF. Nucleotide composition of DNA and RNA from somatic tissues of humpback and its changes during spawning. Biokhimiia (Moscow, Russia). 1967;32(5):988–93.

    CAS  Google Scholar 

  3. Wilson VL, Jones PA. DNA methylation decreases in aging but not in immortal cells. Science. 1983;220(4601):1055–7. https://doi.org/10.1126/science.6844925.

    Article  CAS  PubMed  Google Scholar 

  4. Shimoda N, Izawa T, Yoshizawa A, Yokoi H, Kikuchi Y, Hashimoto N. Decrease in cytosine methylation at CpG island shores and increase in DNA fragmentation during zebrafish aging. Age (Dordr). 2014;36(1):103–15. https://doi.org/10.1007/s11357-013-9548-5.

    Article  CAS  Google Scholar 

  5. King GD, Rosene DL, Abraham CR. Promoter methylation and age-related downregulation of klotho in rhesus monkey. Age (Dordr). 2012;34(6):1405–19. https://doi.org/10.1007/s11357-011-9315-4.

    Article  CAS  Google Scholar 

  6. Xiao P, Zhong T, Liu ZF, Ding YY, Guan WJ, He XH, et al. Integrated analysis of methylome and transcriptome changes reveals the underlying regulatory signatures driving curly wool transformation in Chinese Zhongwei goats. Front Genet. 2020;10:1263. https://doi.org/10.3389/fgene.2019.01263.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Zhou FC. DNA methylation program during development. Front Biol. 2012;7(006):485–94. https://doi.org/10.1007/s11515-012-9246-1.

    Article  CAS  Google Scholar 

  8. Daugela L, Nüsgen N, Walier M, Oldenburg J, Schwaab R, Osman EM, et al. Measurements of DNA methylation at seven loci in various tissues of CD1 mice. PLoS One. 2012;7(9):e44585. https://doi.org/10.1371/journal.pone.0044585.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009;462(7271):315–22. https://doi.org/10.1038/nature08514.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. KaminenAhola N, Ahola A, Maga M, Mallitt KA, Chong S. Maternal ethanol consumption alters the epigenotype and the phenotype of offspring in a mouse model. PLoS Genet. 2010;6(1):e1000811. https://doi.org/10.1371/journal.pgen.1000811.

    Article  CAS  Google Scholar 

  11. Frommer M, McDonald LE, Millar DS, Collis CM, Watt F, Grigg GW, et al. A genomic sequencing protocol that yields a positive display of 5-methylcytosine residues in individual DNA strands. Proc Natl Acad Sci U S A. 1992;89(5):1827–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Laird PW. Principles and challenges of genome wide DNA methylation analysis. Nat Rev Genet. 2010;11(3):191–203. https://doi.org/10.1038/nrg2732.

    Article  CAS  PubMed  Google Scholar 

  13. Cobb JE, Wong NC, Yip LW, Martinick J, Bosnich R, Sinclair RD, et al. Evidence of increased DNA methylation of the androgen receptor gene in occipital hair follicles from men with androgenetic alopecia. Br J Dermatol. 2011;165(1):210–3. https://doi.org/10.1111/j.1365-2133.2011.10335.x.

    Article  CAS  PubMed  Google Scholar 

  14. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2. https://doi.org/10.1093/bioinformatics/btr167.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Wang YZ, Wu JH, Ma X, Liu B, Dong Y. Single base-resolution methylome of the dizygotic sheep. PLoS One. 2015;10(11):e0142034. https://doi.org/10.1371/journal.pone.0142034.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Rosenfeld CS. Animal models to study environmental epigenetics. Biol Reprod. 2010;82(3):473–88. https://doi.org/10.1095/biolreprod.109.080952.

    Article  CAS  PubMed  Google Scholar 

  17. Sen GL, Reuter JA, Webster DE, Zhu L, Khavari PA. DNMT1 maintains progenitor function in self-renewing somatic tissue. Nature. 2010;463(7280):563–7. https://doi.org/10.1038/nature08683.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Kang XL, Liu G, Liu YF, Xu QQ, Zhang M, Fang MY. Transcriptome profile at different physiological stages reveals potential mode for curly fleece in Chinese tan sheep. PLoS One. 2013;8(8):e71763. https://doi.org/10.1371/journal.pone.0071763.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Liu YF, Zhang JB, Xu Q, Kang XL, Wang KJ, Wu KL, et al. Integrated miRNA-mRNA analysis reveals regulatory pathways underlying the curly fleece trait in Chinese tan sheep. BMC Genomics. 2018;19(1):360. https://doi.org/10.1186/s12864-018-4736-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Cui ZJ, Zhang YL, Jiang Y, Chen GN, Xu Z. Breeding report of tan sheep. Part 1: the relationship of ecology and reproduction of Chinese tan sheep. China Anim Husb Veter. 1962;4:1–5.

    Google Scholar 

  21. Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, et al. Shotgun bisulfite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008;452(7184):215–9. https://doi.org/10.1038/nature06745.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Sheng L, GarrettBakelman FE, Akalin A, Zumbo P, Levine R, To BL, et al. An optimized algorithm for detecting and annotating regional differential methylation. BMC Bioinformatics. 2013;14(5):S10. https://doi.org/10.1186/1471-2105-14-S5-S10.

    Article  Google Scholar 

  23. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19(9):1639–45. https://doi.org/10.1101/gr.092759.109.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Issa JP. Age-related epigenetic changes and the immune system. Clin Immunol. 2003;109(1):103–8. https://doi.org/10.1016/S1521-6616(03)00203-1.

    Article  CAS  PubMed  Google Scholar 

  25. Wang SH, Li F, Liu JW, Zhang YL, Zheng YJ, Ge W, et al. Integrative analysis of methylome and transcriptome reveals the regulatory mechanisms of hair follicle morphogenesis in cashmere goat. Cells. 2020;9(4):969. https://doi.org/10.3390/cells9040969.

    Article  CAS  PubMed Central  Google Scholar 

  26. Hadad N, Masser DR, Blanco-Berdugo L, Stanford DR, Freeman WM. Early-life DNA methylation profiles are indicative of age-related transcriptome changes. Epigenetics Chromatin. 2019;12(1):58. https://doi.org/10.1186/s13072-019-0306-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Hu S, Wan J, Su Y, Song Q, Zeng Y, Nguyen HN, et al. DNA methylation presents distinct binding sites for human transcription factors. eLife. 2013;2(2):e00726. https://doi.org/10.7554/eLife.00726.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Zhao M, Liang G, Wu X, Wang S, Zhang P, Su Y, et al. Abnormal epigenetic modifications in peripheral blood mononuclear cells from patients with alopecia areata. Br J Dermatol. 2012;166(2):226–73. https://doi.org/10.1111/j.1365-2133.2011.10646.x.

    Article  CAS  PubMed  Google Scholar 

  29. Spruijt CG, Gnerlich F, Smits AH, Pfaffeneder T, Jansen PTC, Bauer C, et al. Dynamic readers for 5-(hydroxy) methylcytosine and its oxidized derivatives. Cell. 2013;152(5):1146–59. https://doi.org/10.1016/j.cell.2013.02.004.

    Article  CAS  PubMed  Google Scholar 

  30. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat. Rev. Genet. 2012;13(7):484–92. https://doi.org/10.1038/nrg3230.

    Article  CAS  PubMed  Google Scholar 

  31. Petridis C, Navarini AA, Dand N, Saklatvala J, Baudry D, Duckworth M, et al. Genome-wide meta-analysis implicates mediators of hair follicle development and morphogenesis in risk for severe acne. Nat Commun. 2018;9(1):5075. https://doi.org/10.1038/s41467-018-07459-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Öztürk ÖA, Pakula H, Chmielowiec J, Qi J, Stein S, Lan L, et al. Gab1 and Mapk signaling are essential in the hair cycle and hair follicle stem cell quiescence. Cell Rep. 2015;13(3):561–72. https://doi.org/10.1016/j.celrep.2015.09.015.

    Article  CAS  Google Scholar 

  33. Nam KS, Oh S, Shin I. Ablation of CD44 induces glycolysis-to-oxidative phosphorylation transition via modulation of the c-Src-Akt-LKB1-AMPKα pathway. Biochem J. 2016;473(19):3013–30. https://doi.org/10.1042/BCJ20160613.

    Article  CAS  PubMed  Google Scholar 

  34. Nissimov JN, Das Chaudhuri AB. Hair curvature: a natural dialectic and review. Biol Rev. 2014;89(3):723–66. https://doi.org/10.1111/brv.12081.

    Article  PubMed  Google Scholar 

  35. Stadler MB, Murr R, Burger L, Ivanek R, Lienert F, Schler A, et al. DNA-binding factors shape the mouse methylome at distal regulatory regions. Nature. 2011;480(7378):490–5. https://doi.org/10.1038/nature10716.

    Article  CAS  PubMed  Google Scholar 

  36. Bock C, Beerman I, Lien WH, Smith ZD, Gu H, Boyle P, et al. DNA methylation dynamics during in vivo differentiation of blood and skin stem cells. Mol Cell. 2012;47(4):633–47. https://doi.org/10.1016/j.molcel.2012.06.019.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Kuramoto T, Hirano R, Kuwamura M, Serikawa T. Identification of the rat rex mutation as a 7-bp deletion at splicing acceptor site of the krt71 gene. J Vet Med Sci. 2010;72(7):909–12. https://doi.org/10.1292/jvms.09-0554.

    Article  CAS  PubMed  Google Scholar 

  38. Lee JR, Hong CP, Moon JW, Jung YD, Kim DS, Kim TH, et al. Genome-wide analysis of DNA methylation patterns in horse. BMC Genomics. 2014;15(1):598. https://doi.org/10.1186/1471-2164-15-598.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Kimura R, Watanabe C, Kawaguchi A, Kim YI, Park SB, Koutaro M, et al. Common polymorphisms in WNT10A affect tooth morphology as well as hair shape. Hum Mol Genet. 2015;9(9):2673–80. https://doi.org/10.1093/hmg/ddv014.

    Article  CAS  Google Scholar 

  40. Ewelina P, Karłowska-Pik J, Marcińska M, Sarah A, Jeppe DA, van den Berge M, et al. Evaluation of the predictive capacity of DNA variants associated with straight hair in Europeans. Forensic Sci. Int.Genet. 2015;19:280–8.

    Google Scholar 

  41. Cadieu E, Neff MW, Quignon P, Walsh K, Chase K, Parker HG, et al. Coat variation in the domestic dog is governed by variants in three genes. Science. 2009;326(5949):150–3. https://doi.org/10.1126/science.1177808.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Gandolfi B, Outerbridge CA, Beresford LG, Myers JA, Pimentel M, Alhaddad H, et al. The naked truth: sphynx and Devon rex cat breed mutations in krt71. Mamm Genome. 2010;21(9-10):509–15. https://doi.org/10.1007/s00335-010-9290-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Carroll LS, Capecchi MR. Hoxc8 initiates an ectopic mammary program by regulating Fgf10 and Tbx3 expression and Wnt/β-catenin signaling. Development. 2015;142(23):4056–67. https://doi.org/10.1242/dev.128298.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Price FD, Yin H, Jones A, van Ijcken W, Grosveld F, Rudnicki MA. Canonical Wnt signaling induces a primitive endoderm metastable state in mouse embryonic stem cells. Stem Cells. 2013;31(4):752–64. https://doi.org/10.1002/stem.1321.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Professor Changxin Wu in the Department of Animal Genetics and Breeding at China Agricultural University for his help with this study.

Funding

This research was supported by the talent cultivation and developmental support program of China Agricultural University, an award to study the cultivation of high-quality mutton sheep varieties (or lines) from Ningxia province (NXNYYZ20150101), and the Natural Science Foundation of Hebei Province of China for Youths (C2019402261).

Author information

Author notes

  1. Yufang Liu and Qiao Xu contributed equally to this work.

Authors and Affiliations

  1. Department of Animal Genetics and Breeding, National Engineering Laboratory for Animal Breeding, MOA Laboratory of Animal Genetics and Breeding, College of Animal Science and Technology, China Agricultural University, No. 2 Yuanmingyuan West Rd, Beijing, 100193, People’s Republic of China

    Yufang Liu, Qiao Xu, Jve Wang & Meiying Fang

  2. College of Life Sciences and Food Engineering, Hebei University of Engineering, Handan, 056021, People’s Republic of China

    Yufang Liu & Ying Bai

  3. Biotechnology Institute, Nanchang Normal University, Nanchang, 330029, People’s Republic of China

    Qiao Xu & Dengzhen Feng

  4. College of Agriculture, Ningxia University, Yinchuan, 750021, People’s Republic of China

    Xiaolong Kang

  5. College of Animal Science and Veterinary Medicine, Henan Agricultural University, Zhengzhou, 450002, People’s Republic of China

    Kejun Wang

  6. Beijing Key Laboratory for Animal Genetic Improvement, Beijing, 100193, People’s Republic of China

    Meiying Fang

Authors
  1. Yufang Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Qiao Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xiaolong Kang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Kejun Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Jve Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Dengzhen Feng
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Ying Bai
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Meiying Fang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

MYF designed the study. YFL, QX, XLK, KJW, JW, DZF, and YB collected the samples. YFL generated the data. MYF and YFL analyzed the data. YFL, MYF, QX, and XLK contributed to the writing of the manuscript. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Meiying Fang.

Ethics declarations

Ethics approval and consent to participate

Experimental procedures were approved by the Animal Welfare Committee in the State Key Laboratory for Agrobiotechnology at China Agricultural University (approval number XK257) and this study was carried out in strict accordance with the guidelines and regulations established by this committee.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Supplementary Information

Additional file 1 Table S1.

List of differentially methylated regions among three comparisons (mon1 vs. mon24; mon1 vs. mon48; mon24 vs. mon48).

Additional file 2 Table S2.

GO analysis of DMR-associated genes.

Additional file 3 Table S3.

KEGG pathway analysis for DMR-associated genes.

Additional file 4 Fig. S1.

Motif of mCG sites in the genome. Methylation preferences in 9 bp spanning CG, CHG, and CHH methylcytosine sites. mon1 (mon1–1, mon1–2, mon1–3, mon1–4), curly fleece. mon24 (mon24–1, mon24–2, mon24–3, mon24–4), intermediate phenotype. mon48 (mon48–1, mon48–2, mon48–3), uncurled fleece. H = A, C, or T. The abscissa is the base number of the methylation site and the total height of each position is the sequence conservation of the base, which represents the relative frequency of the base at that position.

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Y., Xu, Q., Kang, X. et al. Dynamic changes of genomic methylation profiles at different growth stages in Chinese Tan sheep. J Animal Sci Biotechnol 12, 118 (2021). https://doi.org/10.1186/s40104-021-00632-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40104-021-00632-9

Keywords

Journal of Animal Science and Biotechnology

ISSN: 2049-1891