Comparative DNA methylome analysis of estrus ewes reveals the complex regulatory pathways of sheep fecundity

Background/aims Sheep are important livestock with variant ovulation rate and fertility. Dorset sheep is a typical breed with low prolificacy, whereas Small Tail Han sheep with FecB mutation (HanBB) have hyperprolificacy. Our previous studies have revealed the gene expression difference between the ovaries from Dorset and HanBB sheep contributes to the difference of fecundity, however, what leads to these gene expression difference remains unclear. DNA methylation, an important epigenetic process, plays a crucial role in gene expression regulation. Methods In the present study, we constructed a methylated DNA immunoprecipitation combined with high throughput sequencing (MeDIP-seq) strategy to investigate the differentially methylated genes between the Dorset and HanBB ovaries. Results Our findings suggest the genes involved in immune response, branched-chain amino acid metabolism, cell growth and cell junction were differentially methylated in or around the gene body regions. Conclusions These findings provide prospective insights on the epigenetic basis of sheep fecundity.


Introduction
Sheep (Ovis aries) is an important livestock in China. One of the major challenges in sheep raising industry is how to increase the sheep fertility. Fertility depends on many factors: estrus cycle, ovulation rate and litter size, etc.; among them, the fecundity of ewe plays a central role. Ovarian folliculogenesis, is precisely regulated by the interaction of the central nervous system, the pituitary and the ovary [1,2]. The environmental factors such as nutritional level and metabolic activity affect sheep reproduction well [3][4][5]. More importantly, the genetic background of ewe is in charge of its fertility. In the latest decades, a serial of studies have demonstrated the critical roles of particular genes in ovarian follicular growth. The ovine Booroola fecundity gene (FecB) was the first identified gene involved in ovulation rate and litter size, and the FecB mutant is believed to be coded by sheep bone morphogenetic receptor type 1 B (BMPR1B) gene [6,7]. BMPR1B belongs to the TGFbeta superfamily, and several other TGF-beta family members are also thought to be associated with the sheep fertility, including growth differentiation factor 9 (GDF9) and bone morphogenetic protein 15 (BMP15) [8,9]. These findings indicated the crucial role of genetic variations in sheep fertility.
Recent genomic studies provided a prospective view on the hereditary basis of sheep fertility. The first complete genome sequence of sheep was generated from two Texel individuals in 2014, and their findings indicated the key genes such as MOGAT2 and MOGAT3 involved in long chain fatty acid metabolism [10]. The sheep genome sequence difference could help to understand the genetic foundations of sheep fertility [11]. However, the DNA sequence is not the only thing to connect the gene and gene products. Epigenetics, how physiological phenotypic variations or environmental factors affect gene expression without gene sequence change, will also affect the sheep fertility. Typical epigenetic regulation includes non-coding RNA, chromatin variation and DNA methylation [12]. Several studies have focused on the differentially expressed microRNAs between sheep breeds with fertility variations [13][14][15][16]. In current study, we aim to explore the roles of DNA methylation in the epigenetic regulation of ewe fecundity.
Small Tail Han sheep is a widely bred sheep breed in China with year-round estrus. The Han sheep is wellknown for its hyperprolificacy, the mean litter size of Han sheep is 2.61 [17]. FecB mutant was identified in most of high prolific Han sheep [18]. On the other hand, Dorset sheep is a typical sheep breed with low fecundity, the mean litter size of Dorset sheep is only 1.45 [19]. A comparison of the hereditary background of Dorset sheep and Han sheep will provide important information to get a better grasp of sheep fertility. Ovary plays the central role in ewe fecundity. Our previous study has compared the ovarian mRNA and miRNA expression difference between Han sheep and Dorset [14]. However, the differentially expressed miRNAs cannot explain the mRNA expression difference completely. In the present study, in order to get a deeper understanding of the role of ovarian epigenetic regulation in sheep fertility, we construct a methylated DNA immunoprecipitation combined with high throughput sequencing (MeDIP-seq) strategy to identify the ovarian genes with different DNA methylation level between Dorset sheep and Han sheep with FecB mutant. As we know, this is the first comparative study to investigate the DNA methylome distribution in sheep ovary. Our findings will contribute to the understanding of the epigenetic mechanisms of sheep fertility.

Global mapping of DNA methylation in sheep ovarian genome
In the present study, we mapped the global DNA methylation status of ovarian genome collected from Dorset and Han sheep using MeDIP-seq. A total of 5.53 Gb and 5.63 Gb data were obtained from the ovarian samples of Dorset and HanBB sheep respectively, including 112,892, 144 (Dorset) and 114,963,972 (Han) raw reads. Of the total reads, 93.71 and 93.88% were mapped to the sheep reference genome (ftp://ftp.ncbi.nih.gov/genomes/Ovis_ aries) for Dorset and HanBB group, respectively. Among them, 55.18 and 55.82% were uniquely mapped to sheep genome, respectively.
Uniquely mapped MeDIP-seq reads were detected in all chromosomes of ewe. (GGA1-26 and chromosome X) (Fig. 1). Genome coverage analysis illustrated that most mapped regions had low sequencing depth (1)(2)(3)(4), and no more than 10% of mapped genome had higher sequencing depth (5)(6)(7)(8)(9)(10) in both groups. Additionally, the genome coverage of the CpG sites, and non-CpG methylation such as CHG, and CHH (with H being A, C, or T)) sites were similar between two groups. Around 35% of CpG sites, and 20% of CHG or CHH sites had only one sequencing depth in both groups, and only a small number of regions had high (5-10) sequencing depth in CpG, CHG and CHH sites (Fig. 2).
The mapped reads distributed on different genomic regions exhibits the methylation pattern. Over 16% of uniquely mapped reads were present at the gene body regions in both groups, followed by CpG islands (CGIs, around 11%) and introns (~9%). Further analysis illustrated that within the intragenic elements, the DNA methylation level gradually increased at the upstream region near the transcription start sites (TSSs), peaked at the middle of intragenic region and had another smaller peak just before the transcription termination site (TTS), and dropped dramatically and remained a level of DNA methylation similar as the TSS at the downstream 2 k area (Fig. 3). Many CpG islands located in repetitive elements, and we observed different methylation levels in different repetitive elements, with high methylation in LINE/RTE-BovB, LINE/L1 and SINE/BovA.

DNA methylation profile in sheep genome
To investigate the genome-wide DNA methylation profiles of sheep genome, the uniquely mapped reads were used to detect the methylated peaks, which means the areas in a genome that have been enriched with aligned reads with methylated cytosine. Using MACS1.4.0, we obtained 217, 427 and 208,732 methylated peaks in Dorset and HanBB group, covering 8.60 and 9.02% of sheep genome respectively. Most of peaks had 5 to 20 CpG sites. The peak distribution in different components of the genome were then further analyzed. The majority of peaks covered the coding sequence (CDS) regions followed by 2 kb downstream of TTS region, 2 kb upstream of TSS region and 3′-UTR. Intron and 5′-UTR had much fewer peaks (Fig. 4). Here, upstream2k were regarded as the proximal promoter.

Differential DNA methylation genes in sheep ovary
In the present study, we compared the differentially methylated regions (DMRs) between Dorset and HanBB sheep. Genes with methylation peaks identified in either the promoter or gene body regions were considered as methylated genes. A total of 10,821 differentially methylated genes were successfully identified.

MeDIP-seq data validation
We constructed a bisulfite sequencing strategy to validate the DMRs identified from MeDIP-seq. The result of gene BMPR1B-F from the bisulfite sequencing confirmed that the DNA methylation levels in gene BMPR1B were different between two sheep. (Fig. 6A, B).

GO and KEGG analysis of differentially methylated genes
To identify the biological functions of identified ovarian differentially methylated genes between Dorset and HanBB sheep, we first applied GO annotation analysis. GO annotation indicated that the methylated genes involved in three terms: biological process, cellular component, and molecular function. We mainly focused on the biological process. Compared with the Dorset sheep, among the genes with DMR in upstream2k, downmethylated genes in HanBB sheep were enriched in innate immune response, and up-methylated genes were involved in monocyte chemotaxis, positive regulation of axon extension and regulation of developmental growth. In 5′-UTR elements, no genes with DMR were enriched in any processes with significance. In CDS elements, the downmethylated genes were involved in developmental cell growth, synapse assembly, collateral sprouting and activation of MAPKK activity, while the up-methylated genes were enriched in regulation of histone methylation and DNA methylation on cytosine, including the Ovis aries DNA (cytosine-5-)-methyltransferase 3 beta (DNMT3B), transcript variant 1-6. In intron elements, which is the major region of DMR distribution, the down-methylated genes were related to cellular protein modification process, locomotion, ion transport and cell projection organization, as well as the up-methylated genes were also connected to cellular protein modification process and ion transport, besides, other processes such as G-protein coupled receptor signaling pathway and microspike assembly were involved. In 3′-UTR, the granulocyte differentiation were the only significant enriched process of genes with DMR, which were up-methylated in HanBB sheep. And in 2 kb downstream of TTS region, downmethylated genes were related to carbohydrate catabolic process, and up-methylated genes had no significant enriched processes. Taken together, the genes with DMR were primarily involved in biological processes such as cell signaling, protein modification, ion transport and immune regulation.
To further understand the functions of genes with DMR, we performed KEGG pathway analysis to predict the potential significant pathways (q < 0.05) related to differentially methylated genes. Compared with Dorset group, in HanBB sheep, the genes with up-methylated upstream2k region were associated with valine, leucine and isoleucine biosynthesis; the down-methylated 5′-UTR genes were related to T cell receptor signaling pathway; the genes with down-methylated CDS were involved in ECM-receptor interaction, tight junction, renin-angiotensin system and toxoplasmosis, as well as the up-methylated CDS were also related to ECMreceptor interaction, besides, the pathways of amoebiasis, cysteine and methionine metabolism and protein digestion and absorption were involved. In the intron elements, both the down-methylated and up-methylated genes were connected to tight junction, vascular smooth muscle contraction and long-term potentiation significantly, besides, the down-methylated genes were involved in ErbB signaling pathway. Taken together, the genes with DMR were enriched in pathways associated with cell-cell junction.

DNA methylation and gene expression level
In our previous study, we have explored the different transcriptome of ewe ovaries between Dorset and HanBB sheep (Additional file 2) [20]. To investigate the ovarian gene expression level in current samples, we constructed a RNA-seq strategy to detect the differentially expressed genes between Dorset and HanBB ovaries [20]. To validate the RNA-seq data, the mRNA expression level of krit1 were measured by real-time PCR, the PCR results were consistent with the RNA-seq data, suggests the RNA-seq results were reliable (Fig. 6C).
DNA methylation plays a critical role in gene expression modification. To explore whether the differentially methylated region in or around gene body affects the gene expression, the DNA methylation patterns were co-analyzed with the gene expression profiles. A total of 63 differentially expressed genes between Dorset and HanBB sheep ovaries were identified with DMR in or around the gene body, and the major DMR of these genes were located in the intron regions (Additional file 3). GO analysis indicated that "the hydrolase activity, acting on ester bonds" was the only significant molecular function term, while the biological process and cellular component had no

Discussion
The variation of sheep fecundity depends on the genetic background and gene expression level of ewe ovary. Our previous studies have revealed the differential ovarian gene expression between low-prolificacy Dorset sheep and high-prolificacy Small Tail Han sheep with FecB genotype [14]. FecB gene mutation leads to higher ovulation rate and litter size in sheep [18]. So the comparison between Dorset sheep and Han sheep with FecB mutation provided an excellent model to investigate the influence of genetic backgrounds of ovary on fertility. Our previous RNA-seq study have pointed out that those genes involved in cell proliferation, steroid metabolic processes, ribosome assemble and stress response may contribute to the hyperprolificacy of HanBB sheep [14]. However, what leads to the difference of these gene expression remain unclear. Epigenetics has been believed to play an important role in gene expression regulation. Epigenetics developed quickly in the latest decade, based on the development of next-generation sequencing technique [21]. Histone modification, DNA methylation and non-coding microRNAs are all possible to affect the gene expression and thereby induce physiological modification. Our previous studies have demonstrated the roles of microRNAs on sheep physiology such as adipogenesis and lipid metabolism, muscle development and reproduction [14,[22][23][24][25]. However, how the DNA methylation affects gene expression in sheep, particularly in ovary mature, ovulation and fecundity is still unknown. In the present study, we conducted a MeDIP-seq experiment to investigate the ovarian genome methylation difference between Dorset and HanBB sheep. Limited to our knowledge, this is the first study to explore the role of DNA methylome in sheep fecundity.
DNA methylation, refers to the methylated modification of 5′ position of cytosine residues in nucleotide, primarily occurred at the C base within CpG dinucleotides, the CpG rich region in genome sequence are thereby termed as CpG islands 30 years ago [26]. Recent studies found that the cytosine residues with in CHG or CHH (with H could be A, C, or T) sites could be methylated too [27]. In the current study, the genomic methylated cytosine in both CpG sites and CHG or CHH sites were measured by the technique known as MeDIP-seq. In MeDIP-seq, anti-5-methylcytosine antibody was used to recognize and pull-down the DNA fragments with methylated cytosines, and then these enriched methylated DNA fragments were sequenced by next-generation sequencing. In our present study, the sequencing reads covered all the chromosomes, suggests the MeDIP-seq experiment works well.
DNA methylation affects various crucial biological processes, such as embryonic development, cellular differentiation, tissue-specific gene expression, genomic imprinting, suppression of repetitive elements, carcinogenesis, aging, X chromosome inactivation, and chromosome stability, via regulating the gene expression and silencing [28][29][30]. The DNA methylation in different genomic regions have different effects on gene expression. The DNA methylation occurred in the promoter region normally leads to gene silencing via blocking the accessing and target site binding of DNA binding proteins or recruiting the methyl-CpG binding proteins to repress the gene transcription, whereas the demethylation of CpG islands in promoter region may result to gene expression [28,31]. Additionally, the DNA methylation in the gene body often promotes gene expression, which may partly attribute to transcript splicing regulation [32]. In our current study, the distribution of all methylated reads illustrated that the gene body are the most methylated regions in sheep ovarian genome, suggests only the intragenic regions were hyper-methylated. Those genes with methylation enriched peaks in or around gene body were thought to be DNA methylation affected genes. In these genes, the highly methylated gene elements were CDS, followed by 2 kb upstream of TSS, indicates both the coding exon regions and promoters were methylated in general. Then we analyzed the genes with DMRs in or around gene body between Dorset sheep and HanBB sheep. Interestingly, in these differentially methylated genes, most of DMRs were in enriched in intron, followed by CDS. Other DMRs were mainly in 5′-UTR and upstream2k regions. The GO and KEGG annotation analysis indicates that in potential promoter regions such as upstream2k, the differentially methylated genes were mainly related to immune regulation processes such as innate immune response and monocyte chemotaxis, and the branched-chain amino acids (BCAAs) metabolism were also involved. Healthy inflammatory response improves folliculogenesis, oocyte maturation and ovulation, whereas unhealthy chronic, low-grade inflammation (metaflammation) harms ovulation [33], our findings implied that ewe's ovulation may be regulated via DNA methylation modification. Nutrition intake such as BCAAs have important influence on sheep fertility [34], and our study implies that utilization of BCAAs may be a reason of different prolificacy between Dorset and HanBB sheep. In the gene body regions such as CDS and intron, the differentially methylated genes were related to cell growth, protein modification, ion transport, and particularly cell junction. The cell junction take a role in ovarian physiology, such as blood-follicle barrier and endothelial permeability, and affects the initiation of follicle growth [35,36]. The different levels of gene body methylation in cell junction related proteins provides some new information about the role of cell-cell junction in sheep fecundity.
The co-analysis of RNA-seq and MeDIP-seq revealed that many genes with higher expression level in HanBB sheep had lower methylated intron level, compared with Dorset sheep, such as krit1. DNA demethylation in intron may lead to up-regulated gene expression via increasing the intronic enhancer activity [37]. Krit1 takes part in SOD2-dependent catabolism of intracellular reactive oxygen species [38]. Many evidences indicate that oxidative stress has harmful influence on female reproduction (reviewed in [39]), the demethylated intron-associated enhanced expression of krit1 may contribute to the high fecundity in HanBB sheep.

Conclusion
To explore epigenetic mechanisms of the hyperprolificacy of Small Tail Han sheep with FecB genotype, we here constructed a MeDIP-seq based comparative DNA methylomics strategy to investigate the differentially methylated genes between the ovaries from HanBB and Dorset sheep, and our findings provide prospective insights on the differences of genomic DNA methylation level between hyperprolific HanBB sheep and relative low prolific Dorset sheep. The roles of DNA methylation in sheep fecundity are worthy further study.

Sheep sample preparation
The tissue samples used in the present study were the same as our previous study [14]. The adult Han ewes with the FecB mutation in the BMPR1B genotype BB were regarded as Han group, which has the high fecundity; meanwhile, The adult Dorset ewes were used as a control low-reproducing Dorset group. All animals were kept under similar conditions and provided water and food ad libitum. All ewes in the experiment were treated to synchronize estrus as previous studies [14]. Twentyfour hours after spontaneous estrus were detected, all ewes were euthanized and whole ovaries were excised to obtain better ovulation points on the surfaces of the ovaries. All samples were immediately snap-frozen in liquid nitrogen and stored at − 70°C for genome DNA and ovarian RNA extraction.

DNA extraction
Genomic DNA was extracted from the ovaries with the use of a phenol-chloroform extraction method, and DNA quality was evaluated by agarose gel electrophoresis and spectrophotometer. Genomic DNA isolated from the three Han sheep with FecB mutant were mixed in equal amounts to generate a pooled sample as high prolific Han sample, and genomic DNA isolated from the three Dorset sheep were mixed to generate the low prolific Dorset sample.

MeDIP-seq analysis
Genomic DNA mixtures were sonicated to gain fragments in the range 100-500 bp, and dAs were added to the 3′ ends of the DNA fragments. The adaptors were then ligated to the two ends of each fragment with paired-end DNA sample prep kit (Illumina) according to the manufacturer's protocol. The double-stranded DNA was denatured to single strands, and then the anti-5methylcytosine mouse monoclonal antibody (Calbiochem) was used to immunoprecipitate the DNA fragments containing regions of methylated cytosines with magnetic methylated DNA immunoprecipitation kit (Diagenode). The precipitated DNA fragments were then amplified by PCR and separated by agarose gel electrophoresis. The 200-300 bp fragments were excised and extracted using a gel extraction kit (28,706;Qiagen). The selected DNA fragments were quantified using an Agilent 2100 Analyzer. The quantified library was sequenced in an Illumina HiSeq 2000 at the BGI (Shenzhen, China).

MeDIP-seq sequence alignments and data analyses
Quality control of raw data from sequencing was performed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The reads containing adapter, or containing over 10% of N nucleotides, or over 50% of contained nucleotides had quality value Q < 21 were removed from the raw data, and all left reads were saved as clean reads in FASTQ files. The clean data were then aligned to the sheep reference genome (Ovis_ aries, ftp://ftp.ncbi.nih.gov/genomes/Ovis_aries) with SOAPaligner v 2.21 (http://soap.genomics.org.cn/) with no more than 2 bp mismatches. The reads distribution in sheep chromosomes and in the different genome components including upstream2k, 5′-untranslated region (UTR), 3′-UTR, CDS, introns, CpG islands (CGIs), downstream2k, and repeats were analyzed; as well as the genome coverage of the CG, CHG, and CHH sites under different sequencing depths and distributions of MeDIP-Seq reads in different CG density regions.
Methylation peaks (known as genome-wide highly methylated regions: regions with sequencing tags more than 20, and a p value < 1 × 10 − 5 ) scanning was conducted with MACS 1.4.0 (http://liulab.dfci.harvard.edu/ MACS/) [40]. The distribution of peaks in different components of sheep genome (upstream2k, 5'UTR, CDS, Intron, 3'UTR, and downstream2k) were analyzed. To identify the differentially methylated peaks between Dorset and Han group, peak regions were assembled and the number of reads of each sample was calculated. The numbers of reads in each peak were assessed using chi-square, and false discovery rate (FDR) statistical method was used to avoid false positive data. p ≤ 0.01 was considered significant. Only those peaks with reads covered in the same genomic components from both groups, and the difference in reads numbers between two groups greater than twofold with significance were identified as differentially methylated peaks. All genes within differentially methylated peaks were annotated using gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The terms and pathways with FDR < 0.05 were thought to be significant.

RNA-seq and real-time polymerase chain reaction (PCR) validation of differentially expressed genes (DEGs)
RNA-seq data used in the present analysis was published in our previous study. To verify the DEGs identified in RNA-seq, real-time PCR was applied, the procedure was similar as previous study [14]. Total RNA was extracted from each ovaries by using TRIZOL reagent (Invitrogen, Carlsbad, CA, USA) and reverse transcription was conducted using Thermo First cDNA Synthesis Kit (Sino-Gene) according to manufacturer instructions. The total RNA (1 μg) from each ovary sample was retrotranscribed within the presence of oligodeoxythymidylic acid primers. Real-time PCR was performed using the StepO-nePLUS (Applied Biosystems). Each qPCR sample was run in a 15-μl total volume with SG Green qPCR Mix (SinoGene). The following thermocycling conditions were used: 10 min at 95°C and 45 cycles of 95°C for 15 s and 60°C for 15 s, followed by a melting ramp from 60 to 95°C, holding for 45 s on the first step (60°C), followed by 15 s holds at 95°C. All reactions were performed in triplicates. GAPDH was used as the internal control. The data from the real-time PCR was analyzed by the 2 -ΔΔCT method [41]. The primer for gene krit1 was as followed, Forward: 5′-CAGAATGGATGGAT CATACCG-3′, Reverse: 5′-CTGCGTTTCTTGAGAG AGAC-3′.

Bisulfite sequencing
To validate the MeDIP-seq data, the gene BMPR1B-F was selected for bisulfite sequencing. 500μg of individual DNA sample from each group was treated with the EZ DNA Methylation-Gold Kit (ZYMO) and then incubated at 95°C for 10 min and 2.5 h at 64°C for bisulfite conversion. After the reaction, the modified DNA samples were extracted for bisulfite sequencing PCR, and the forward primer: 5′-ATAGGGTTGTTGGGTAGTGAGT T-3′, and reverse primer: 5′-AATCTTTCTTTCCCTA CCCTCTAT-3′. The PCR were performed in a 25 μL reaction mixtures that contain 2 μL modified DNA template, 1 μL primers for forward and reverse strand each, 0.15 μL Taq enzyme, 2.5 μL 10X PCR buffer, 1 μL Mg 2+ (50 mM), 0.5 μL dNTP (10 mM) and 17.85 μL ddH2O with the following program: 94°C for 2 min, followed by 10 cycles of denaturation at 94°C for 30 s, annealing at 60°C to 50°C Δ1°C temperature ladder for 30 s, extension at 72°C for 30 s, and then followed by 30 cycles of 94°C for 30 s, 50°C for 30 s, 72°C for 30 s, and a final extension at 72°C for 5 min. The PCR products were purified and subcloned into pEASY -T1 vector (transgen) for sequencing.

Statistical analysis
All data are presented as the mean ± SD. The significance of differences in data between the groups was determined by one-way ANOVA analysis of variance followed by the student's t-test for equality of variances using SPSS 17.0 (IBM, USA). Differences were considered statistically significant at p < 0.05.