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)
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).
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.
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.
Varied body structure and organ function is an unavoidable and complex process associated with aging, and results in increased morbidity and mortality . 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 . 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 . 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 .
DNA methylation is an important epigenetic mechanism that is involved in regulation of various biological processes . Methylation in the 5′- flanking region of a gene plays a pivotal role in regulating expression , and epigenetic markers silence exogenous transposons and imprint genes . This mechanism may produce diverse phenotypes that are capable of responding to environmental changes during maturation . 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 . 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 .
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 . 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 ; 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.
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 , 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) . 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  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. .
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) . 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. .
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 . The KEGG database was used to test differentially expressed genes for statistically significant enrichment in biochemical pathways . Interaction networks for selected DMGs were analyzed using the String database (http://string-db.org/) .
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.
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.
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.
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.
To visualize the distribution of methylated sites as a function of genome location, we used the graphical tool Circos  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.
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.
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).
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.
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).
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.
Epigenetic modifications, such as changes in global and locus-specific DNA methylation, are suspected to play an important role in aging . One possible mechanism is DNA methylation, which is involved in age-dependent animal phenotypes and age-related physiological variation . 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 . 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 . 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 . 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 . 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 . SNP analysis revealed that common variations in WNT10A have pleiotropic effects on hair morphology . Moreover, other researchers found that WNT10A may participate in catagen by inhibiting keratinocyte proliferation and functioning as a WNT inhibitor . 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 . 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 . TBX3 is another gene that is involved in the WNT pathway . 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.
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.
Differentially methylated region
False discovery rate
Differentially methylated gene
Kyoto encyclopedia of genes and genomes
Polymerase chain reaction
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.
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.
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.
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.
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.
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.
Zhou FC. DNA methylation program during development. Front Biol. 2012;7(006):485–94. https://doi.org/10.1007/s11515-012-9246-1.
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.
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.
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.
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.
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.
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.
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.
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.
Rosenfeld CS. Animal models to study environmental epigenetics. Biol Reprod. 2010;82(3):473–88. https://doi.org/10.1095/biolreprod.109.080952.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Ö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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
We thank Professor Changxin Wu in the Department of Animal Genetics and Breeding at China Agricultural University for his help with this study.
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).
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
The authors declare that they have no competing interests.
List of differentially methylated regions among three comparisons (mon1 vs. mon24; mon1 vs. mon48; mon24 vs. mon48).
GO analysis of DMR-associated genes.
KEGG pathway analysis for DMR-associated genes.
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.
About this article
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