Genome-wide analysis of DNA Methylation profiles on sheep ovaries associated with prolificacy using whole-genome Bisulfite sequencing

Abstract

Background

Ovulation rate and litter size are important reproductive traits in sheep with high economic value. Recent work has revealed a potential link between DNA methylation and prolificacy. However, a genome-wide study that sought to identify potential DNA methylation sites involved in sheep prolificacy indicated that it is still unknown. Here, we aimed to investigate the genome-wide DNA methylation profiles of Hu sheep ovaries by comparing a high-prolificacy group (HP, litter size of three for at least 2 consecutive lambings) and low prolificacy group (LP, litter size of one for at least 2 consecutive lambings) using deep whole-genome bisulfite sequencing (WGBS).

Results

First, our results demonstrated lower expression levels of DNA methyltransferase (DNMT) genes in the ovaries of the HP group than that in the ovaries of the LP group. Both groups showed similar proportions of methylation at CpG sites but different proportions at non-CpG sites. Subsequently, we identified 70,899 differential methylated regions (DMRs) of CG, 16 DMRs of CHG, 356 DMRs of CHH and 12,832 DMR-related genes(DMGs). Gene Ontology (GO) analyses revealed that some DMGs were involved in regulating female gonad development and ovarian follicle development. Finally, we found that 10 DMGs, including BMP7, BMPR1B, CTNNB1, FST, FSHR, LHCGR, TGFB2 and TGFB3, are more likely to be involved in prolificacy of Hu sheep, as assessed by correlation analysis and listed in detail.

Conclusions

This study revealed the global DNA methylation pattern of sheep ovaries associated with high and low prolificacy groups, which may contribute to a better understanding of the epigenetic regulation of sheep reproductive capacity.

Background

DNA methylation is an epigenetic regulatory mechanism that plays a significant role in mediating biological processes such as gene expression, genomic imprinting, cell differentiation and embryogenesis, as well as in determining phenotypic plasticity in organisms [1]. DNA methylation occurs at the cytosine residue of CpG dinucleotides, which are unevenly distributed throughout the genomes. Recently, the whole genome methylation of genes involved in vital biological functions has been extensively examined in mammalian species [2, 3] by using advanced high-throughput sequencing technologies. The reproductive efficiency of sheep in terms of litter size has an important impact on the economic returns of farmers [4]. Reproductive traits typically have low to medium heritability and do not exhibit a noticeable response to phenotypic selection [4]. Therefore, investigation of the genetic information associated with reproductive ability could efficiently enhance selection. Prolificacy is controlled by ovarian folliculogenesis, a process that is highly regulated by precise proliferation and differentiation events. Recent research has shifted focus to how DNA methylation regulates the initiation of ovarian and sexual maturation. Evidence has revealed that ovarian maturation is regulated by DNA methylation [5, 6]. By profiling the methylome of porcine ovaries, researchers examined the methylation changes during the process of sexual and ovarian maturation in pigs [7]. Similar studies also found that DNA methylation alterations influenced gene expression profiles in the goat hypothalamus during the onset of puberty [8]. Despite these findings, our understanding of DNA methylation patterns associated with prolificacy remains limited. Hu sheep are widely recognized as having early sexual maturity and high prolificacy; however, in recent years, much attention has been focused on meat rather than reproductive traits during the selection process. Reproduction is a complex process, and traits such as litter size are affected by many minor genes and some major genes; thus, understanding the role of DNA methylation in gene function is necessary.

Of the four DNA methylation sequencing technologies: methylated DNA binding domain sequencing, methylated DNA immunoprecipitation sequencing, representation bisulfite sequencing (RRBS) and whole-genome bisulfite sequencing (WGBS), WGBS is similar to whole genome sequencing, except for one detail: bisulfite conversion. It achieves single-base resolution through bisulfite conversion and does not have a base preference as it can yield almost complete information about methylcytosine with excellent specificity and non-sensitivity [9]. Therefore, WGBS is the most comprehensive of the existing methods. In this study, we investigated DNA methylation profiles in the ovaries of high and low prolificacy sheep at 3 years of age during the estrus stage using WGBS technology. Our research systematically analyzed the DNA methylation patterns potentially involved in litter size. In addition, our findings would advance knowledge and understanding of the sheep methylome.

Methods

Animals and tissue collection

A total of 6 non-pregnant ewes with identical lambing records (3 records) were selected and divided into a high prolificacy group (HP: n = 3, litter size = 3) and a low prolificacy group (LP: n = 3, litter size = 1). Progesterone sponges were intravaginally implanted in ewes for estrus synchronization and were removed after 11 days. Then the estrus status of the ewes was checked every day. Ewes were slaughtered within 12 h during the estrus stage, and ovaries were immediately collected and snap-frozen in liquid nitrogen immediately and then stored at −80 °C until use.

DNA, bisulfite treatment and RNA/cDNA preparation

Each whole ipsilateral ovary per sheep was collected, and the genomic DNA was extracted using a Genomic DNA kit (TIANamp, Cat.#DP304–02), and then bisulfite conversion was performed using the EZ DNA Methylation Direct Kit (Zymo Research Corporation, Cat#D5020). Total RNA was isolated using the TRIzol reagent (Invitrogen Corp, Cat.#15,596–026) and dissolved in RNase-Free water (QIAGEN, Cat.#129,112). The quality and quantity of DNA and RNA were determined using a NanoDrop instrument (NanoDrop Technologies, Wilmington, DE, USA). cDNA samples were synthesized from total RNA using a reverse transcription (RT) reagent kit with gDNA Eraser (Takara, Cat.#RR047A). All operations were conducted following the manufacturer’s recommended instructions.

WGBS library preparation and data analysis

Three samples from each group were selected for WGBS sequencing. Genomic DNA was fragmented by ultrasonication. The fragments were then end-repaired, 3′-end-adenylated and ligated with adapters. Agarose gel electrophoresis was used to select fragments of 400–500 bp in length. The selected fragments were treated with bisulfite and subjected to PCR amplification to form the sequencing library. Then, the qualified library was sequenced using an IlluminaHiSeq™2500 system (Biomarker Technologies, Beijing, China). The peak signal produced by the Illumina HiSeq was transformed into base sequence by base calling as Raw Data or Raw Reads. The Raw Reads were then filtered for subsequent information analysis to ensure the quality of information analysis, including the removal of reads that have adapters and filtration of reads with more than 10% N content or more than 50% low quality bases. The final filtered data are called clean reads.

Mapping reads to known genome

The sequencing reads need to be aligned with the reference genome (Oar_v3.1) before conducting the methylation analysis. Bismark software was used to perform a comparison of the alignments of bisulfite-treated reads to a reference genome using the default parameters. Reads that aligned with the same region of the genome were taken as the duplicate number. And the duplicate number was used to summarize the sequencing depth and coverage. The conversion rate of bisulfite was calculated as the percentage of the methylated clean reads as a percentage of the total number of clean reads in the lambda genome by using Bismark software. As unmethylated cytosine from the genome was converted into T after bisulfite treatment and PCR amplification, but methylated cytosine remained unchanged. Bismark was able to extract information about genome cytosine sites from the results of the comparison of the clean reads with the reference genome and thereby acquire cytosine site coverage statistics and the number of different types (CG as CpG, CHG and CHH) of methylated cytosine reads. As the methylation single C site cannot be discriminated by Bismark, we used the binomial distribution test for each C site to confirm the methylated C site by screening conditions for coverage ≥4× and false discovery rate (FDR) < 0.05.

Estimating methylation levels and the identification of DMRs

All cytosine sites with reads coverage >10X were used for DMRs analysis with MOABS [10]. First, to detect the different methylated C sites in a region, we defined C _i as the number of supporting methylation reads at a single C site, T _i as the number of supporting unmethylation reads at a single C site, i as the position of C, and n as the total number of C positions. The methylation level of a C site was counted as follows [11]:

$$ \mathrm{Methylation}\ \mathrm{level}\ \mathrm{of}\ \mathrm{C}\ \mathrm{site}={C}_i/\left({C}_i+{T}_i\right) $$

The binomial distribution test was used to determine whether the C site was methylated. Subsequently, DMRs were defined that three different methylation sites at least in the region, and in which the difference in methylation levels was greater than 0.2 (0.3 for CG type) with p value from Fisher’s exact test of less than 0.05. The methylation level of regions was counted as follows [11]:

$$ \mathrm{Methylation}\ \mathrm{level}\ \mathrm{of}\ \mathrm{region}=\frac{1}{n}\times \sum_{i=1}^n\frac{C_i}{C_i+{T}_i} $$

Bioinformatics analysis of DMGs

The DMGs were compared with functional databases such as GO, COG (Cluster of Orthologous Groups of proteins) and KEGG (Kyoto Encyclopedia of Genes and Genomes) by BLAST to obtain the annotation of these genes for analyzing gene function. The GO enrichment analysis was implemented by the Wallenius non-central hyper-geometric distribution in the GOseq R packages [12]. KOBAS software was used to test for statistically significant enrichment of differentially expressed genes in KEGG pathways [13]. The interaction networks of selected DMGs were analyzed using the String database (http://string-db.org/) [14].

Quantitative reverse transcriptase PCR

Quantitative reverse transcriptase PCR (qRT-PCR) was used to examine DNA methyltransferase gene expression levels and validate the DMGs from the sequencing results by detecting the mRNAs expression level. Ten DMGs were randomly selected, and the specific primers used in the qRT-PCR are listed in Additional file 1: Table S1.

qRT-PCR was performed on a StepOnePlus Real-Time PCR System (Life Technologies, USA). Each PCR (a reaction volume of 20 μL) system included 10 μL of SYBR Green Master mix (Roche Applied Science, Mannheim, Germany), 0.6 μL each of 10 μM forward and reverse primer, 1 μL of reverse transcription product and 7.8 μL of RNase-Free water. The comparative quantification of each of the results was standardized to GAPDH by the 2^-ΔΔCt method.

Statistical analysis

All data were analyzed using the SPSS 24.0 software package (SPSS, Chicago, IL, USA). The qRT-PCR results are expressed as the mean ± standard error, and each group contained three samples and the experiments were repeated at least 3 times. The different mRNA expression levels of genes in the HP and LP groups were compared using the independent-samples t-test. Differences were considered significant at P < 0.05.

Results

Expression levels of DNMTs

The expression levels of DNMTs (DNMT1, DNMT3A and DNMT3B) in ovaries were first analyzed by qRT-PCR in the HP and LP sheep. As Fig. 1 shows, DNMT1 and DNMT3A were expressed at significantly lower levels in the HP group than in the LP group (P < 0.05).

Fig. 1

The mRNA expression level of DNMTs as determinded by qRT-PCR. The relative expressions of DNMTs in ovaries was detected by qRT-PCR. The experiment was performed using 3 biological repeats and 3 technical repeats. The relative expression levels were normalized to the expression amount of GAPDH. The results are expressed relative to the LP group as the mean ± SEM and the ordinate represents log10-transformed values. **, P < 0.01

DNA methylation mapping and patterns

A total of 63.79 G and 66.72 G raw bases were generated on average for the HP group and the LP group respectively. After data filtering, approximately 200 million clean reads were generated for each group. These reads were detected in all chromosomal regions for each group. The mapped reads were used for subsequent analysis as the rates were from 71.36% to 74.68% (Table 1).

Table 1 Whole genome DNA bisulfite sequencing data

In each group, approximately 3.5% of all genomic C sites were methylated (Table 1). Methylation in sheep was found to exist in three sequence contexts: CG, CHG (where H is A, C or T), and CHH. These contexts were present in similar proportions in each group, and we found overall genome-wide methylated cytosine levels of 89.78% CG, 2.46% CHG, 7.76% CHH in the HP group and 88.60% CG, 2.66% CHG, 8.74% CHH in the LP group (Fig. 2).

Fig. 2

The average ratio of DNA methylation types in the genomes of HP and LP. HP, High Prolificacy. LP, Low Prolificacy, H = A, C or T. The green, blue and yellow colors represent mCG, mCHG and mCHH, respectively

Sequence preferences analysis for methylation

A violin graph was plotted with dots representing different methylation levels, and we found that the methylation levels were high with wide sections in the violin plot for CG methylation types (Fig. 3a), but the methylation levels were low with narrow sections in the violin plot for CHG and CHH methylation types (Fig. 3b and c). Then, we plotted chromosome methylation maps for each sample. The results showed that most hypermethylation cytosine were of the CG type in chromosomes and that the mC (methylated cytosine) sites were different on the chromosome like chromosome 18 between the two groups (Additional file 2: Figure S1).

Fig. 3

Violin plot for the overall distribution of methylation levels for different methylation types. a, CG. b, CHG. c, CHH. HP (J07, J08, J09), High Prolificacy. LP (J10, J11, J12), Low Prolificacy. H = A, C or T. The abscissa represents the different samples, the ordinate represents the level of methylation of the samples; and the width of each violin represents the density of the point at that methylation level; while the boxplot shows the methylation levels in each violin

Furthermore, we analyzed the relationship between sequence context and methylation preference. We calculated the percentage methylation for all possible 9-mer sequences in which either the methylated cytosine was in the fourth position (allowing analysis of the three nucleotides upstream of CG, CHG, and CHH methylation), the mC sites were in a hypermethylation state, the CAG was the most common sequence motif in the CHG mC sites, and different frequencies of the CHH contexts was discovered for the two groups (Fig. 4).

Fig. 4

Methylation preferences in 9 bp spanning CG, CHG, and CHH methylcytosine sites. HP (J07, J08, J09), High Prolificacy. LP (J10, J11, J12), Low Prolificacy. H = A, C or T. The abscissa is the base number of the methylation site, the total height of each position is the sequence conservation of the base, which represents the relative frequency of the base at that position

DNA methylation levels of different functional regions

We divided all mC into specific gene features: promoter, 5’UTR (Untranslated Regions), exons, introns, and 3’UTR. The methylation levels were evaluated in these functional regions. As shown in Fig. 5, the trend of methylation levels in the specified regions of the two groups were similar in the two groups, and the methylation levels for the CG type were higher than those for the CHG and CHH types. Intron regions, exons (except for the first exon) and downstream regions are the major components of mC containing sites (Additional file 3: Table S2). Moreover, the methylation levels of CG in the first exon were lower than those of the other elements except for the upstream region, as the levels showed a downward trend in upstream region, and the methylation levels of CG sites near the TSS were lower than those in the first exon. In addition, high levels of DNA methylation were detected in inner exons and introns, and the levels of methylation decreased gradually from the promoters to the TSS (Transcription Start Site) and increased from the TSS to the introns; the CHH type was hypomethylation and stable in each functional element; the CHG type was almost entirely unmethylated. More detailed information is listed in Additional file 4: Table S3.

Fig. 5

The trend in CGI distribution in different functional elements. a, J07. b, J08. c, J09. d, J10. e, J11. f, J12. HP (J07, J08, J09), High Prolificacy. LP (J10, J11, J12), Low Prolificacy. H = A, C or T. The abscissa represents the different regions of gene functional elements that a, b, c, d, e, f, g denote upstream, first exon, first intron, inner exon, inner intron, last exon and downstream, respectively. The left ordinate represents the methylation levels of CG/CHG, and the right ordinate represents the methylation levels of CHH. The dotted, green, vertical line represents the TSS, and the red, orange and blue solid lines represent CG, CHH and CHG, respectively, which show the methylation levels fluctuating in the different regions

Annotation of methylation CGI regions

We counted the quantity of hypermethylation CGI (CG islands) regions (hypermethylation CGI definition: methylation level over 0.7 except when the proportion of C sites with high confidence was less than 0.1) and annotated these with gene functional elements (with 3000 bp 5′ to the TSS and 3′ to transcription termination site as the gene upstream and downstream functional regions respectively). As Fig. 6 shows, approximately 68% hypermethylation CGI was distributed in distal intergenic regions, within 1.5% hypermethylation CGI was distributed in UTR and there was no significant difference between the two groups (P > 0.05).

Fig. 6

The ratio of hypermethylation CGI distribution in different functional elements of gene. HP, High Prolificacy. LP, Low Prolificacy. The charcoal, blue, orange, gray, yellow, hyacinthine, green, dark blue and brown colors represent distal intergenic, promoter (<=3 kb), 5’UTR, 3’UTR, 1st exon, other exon, 1st intron, other intron and downstream regions (<=3 kb), respectively

DMRs analysis for the HP and LP groups

DMRs were detected in the two groups and were annotated into gene functional elements according to different methylation types. In total, 70,899 CG DMRs, 16 CHG DMRs and 356 CHH DMRs were identified, most of which were in distal intergenic regions, with only 33 and 162 DMRs were in the 5’UTR and 3’UTR respectively. For all methylation types, the ratio of DMRs located in introns was the highest except for those in distal intergenic regions (Fig. 7). The difference in methylation levels between the two groups in chr12(74369711–74,369,728), chr14(9182359–9,182,369), chr12(69412005–69,412,018), chr3(176960930–176,960,977) and chr1(230630159–230,630,190) had exceed 90% (Additional file 5: Table S4). A heat map was generated using a cluster analysis of DMRs for the HP and LP groups (Fig. 8). More detailed DMRs results are listed in Additional file 5: Table S4.

Fig. 7

The ratio of DMRs with different methylation types in different gene functional regions. a, CG. b, CHG. c, CHH. H = A, C or T. The blue, dark blue, green, khaki, red, yellow, orange, purple, light yellow and brown colors represent promoter (<=1 kb), promoter (1-2 kb), promoter (2-3 kb),5’UTR, 3’UTR, 1st exon, other exon, 1st intron, other intron and distal intergenic regions, respectively

Fig. 8

Heat map cluster analysis of DMRs in different gene functional regions. Each column represents an individual DMR and each row represents one group. The colors in each block from blue to white to red sequentially represents the methylation ratio from 0 to 0.5 to 1, respectively. In addition, the red, yellow, green, turquoise, blue, purple and pink colors represent the 3’UTR, first intron, inner exon, inner intron, last intron, promoter and distal intergenic regions respectively which are shown above the heatmap

Verification of sequencing results

To further validate the sequencing results, 10 DMGs from the sequencing results were randomly selected for detection by qRT-PCR. As shows in Fig. 9, the GPNMB, ELK4, BACH1, CDIPT levels were significantly lower and the SCYL1 levels were significantly higher in the HP than in the LP group (P < 0.05), and the ABCG2, mTOR, STK3, ACVR1 and PSMD7 levels were not significantly different between the two group (P > 0.05). More detailed information of these genes/DMRs is listed in Additional file 6: Table S5. In total, the qRT-PCR results showed that the sequencing data were reliable.

Fig. 9

The mRNA expression level of DMGs as determined by qRT-PCR. The relative expressions of DMGs in ovaries was detected by qRT-PCR. The experiment was performed using 3 biological repeats and 3 technical repeats. The relative expression levels were normalized to the expression amount of GAPDH. The results are expressed relative to the LP group as the mean ± SEM, and the ordinate represents log10-transformed values. *, P < 0.05. **, P < 0.01

Database enrichment analysis for DMGs: COG, GO and KEGG

To probe changes in the methylation status of gene functions under prolificacy traits, the COG, GO and KEGG pathway databases were analysed to characterize the 12,832 DMGs that were detected in the DMRs. The COG analyses revealed that DMGs were enriched on general function prediction mostly for CG (Fig. 10a). The GO analysis revealed that for the CG type, DMGs were significantly enriched in the categories of cell migration, anatomical structure formation involved in morphogenesis, cell projection, and intracellular membrane-bounded organelle (Fig. 10b). The KEGG analysis revealed that for CG type, DMGs were significantly enriched in the categories mucin type O-Glycanbiosynthesis, long-term depression and nicotine addiction (Fig. 10c). Importantly, we also found that some DMGs were involved in biological processes important for female gonad development, including ovarian follicle development (GO: 0001541), ovulation from ovarian follicle (GO:0001542), antral ovarian follicle growth (GO:0001547), luteinization (GO:0001553), ovulation cycle process (GO:0022602), negative regulation of female gonad development (GO:2,000,195), which suggested that these specific genes, which are influenced by DNA methylation could affect the development of ovarian follicles, subsequently impacting ewes’ prolificacy. For more detailed results of the COG, GO, KEGG analyses of CG, CHG and CHH (see Additional files 7, 8, 9 and 10).

Fig. 10

COG, top GO and top KEGG pathway analysis of CG type DMGs. a, COG analysis. b, top GO analysis. c, top KEGG analysis. In graph a, the abscissa represents the COG function classification; the ordinate represents the number of genes that were enriched in this classification. In graph b, the ordinate represents the GO terms that were the most enriched; the abscissa represents the P-value that was calculated using -log10-transformed values; the green and orange colors indicate biological process and cellular component. The size of the circles represent the number of genes contained in the particular class in the graph c, the larger the circle is, the more genes there are; Differently colored circles represent the enrichment degree of false positives, the redder the circle is, the lower is the false positive rate

Correlation analysis of DMGs and sheep prolificacy

To further understand the relationship between DNA methylation and different levels of prolificacy, we set two limiting factors to perform an association analysis. First, the DMGs of the two groups should be enriched in female reproduction related pathways in the GO analysis. Second, the pathway in KEGG (except for disease and cancer pathways) that was enriched for the selected DMGs should significantly differ between the two groups (P < 0.05). As a result, 28 genes meeting these two criteria were detected (more detailed information on these genes is listed in Additional file 11: Table S7). Subsequently, these genes were analyzed using the STRING database.

As Fig. 11 shows, within the network analysis, we focused on the DMGs that interacted with 5 or more other genes. BMP7, BMPR1B, CTNNB1, FST, FSHR, LHCGR, TGFB2, and TGFB3 are hub genes in the network, related to the female reproduction pathway. More detailed results on the abovementioned genes are listed in Table 2.

Fig. 11

STRING analysis of DMGs associated with prolificacy. The detected DMGs were analyzed using the STRING database, with setting as follows: organism: Bos Taurus, interaction score: medium confidence, 1st shell: none/query proteins only, 2nd shell: none

Table 2 DMGs most likely involved in female prolificacy pathways

Discussion

DNA methylation is the main feature of the epigenetic regulatory mechanism that plays an important role in the regulation of gene expression. Recently, studies have been conducted to identify the genome-wide methylation profiles of farm animals [15,16,17,18]. Previously, some studies have been conducted to describe DNA methylation for sheep ovary [19,20,21,22], but few reports from the ovarian genome-wide methylation pattern [23]. WGBS, which allows unbiased genome wide DNA methylation profiling, has allowed us to investigate prolificacy related DNA methylation in unprecedented detail [9]. In this study, we used WGBS to investigate the DNA methylation profiles of the genome in ovarian tissues of high prolificacy and low prolificacy sheep to discover the relationship between DNA methylation and different levels of prolificacy. Further correlation analysis indicated that several DMR-related genes were most likely involved in Hu sheep prolificacy.

DNMTs are the writers of the epigenome. DNMTs constitute a highly conserved family of proteins in mammals, and there are 3 major DNMTs: DNMT1, DNMT3a and DNMT3b. DNMT1 is a maintenance DNMT, while DNMT3a and 3b are de novo DNMTs. Lynch et al. (2016) indicated that DNMT3a has broad downstream effects on the timing of the genomic control of reproductive function [24]. Our results showed that the mRNA expressions of DNMT1 and DNMT3a in ovary tissue were significantly lower in the HP group than that in the LP group, which indicated that DNMTs may regulate the transcription of genes associated with sheep prolificacy. In our study, approximately 3.5% of cytosine sites were methylated, and the CG methylation type was present in the highest proportion and at the highest level in the genome. These results were similar to those found in other species, including humans and pigs [17, 25]. Non-CpG information was also obtained in the present study, CAG was the most common sequence motif in the CHG mC sites, as also found in other studies [26], and the sequence motif in CHH differed between the two groups. The methylation levels near the TSS were the lowest of all the gene functional regions, which was consistent with the results found in pig ovaries by RRBS [27]. Most hypermethylation CGI (over 68%) were located in distal intergenic regions, while only 1.5% were in UTR.

As the DNA methylation status of promoter and gene body regions could affect gene expression thorough changes in chromatin structure or transcription efficiency [28, 29], we compared the genome-wide methylation patterns of the HP and LP sheep to identify DMGs that may affect prolificacy. We identified 71,271 DMRs and 12, 831 genes related to these DMRs were predicted, and 68.60% of the DMRs were located in distal intergenic regions, but only 0.27% of the DMRs were located in UTRs, a finding that was also similar to the previous research in pig ovary tissues by RRBS [27]. To further validate the sequencing results, qRT-PCR was performed to detect the mRNA expression of 10 randomly selected DMGs, the expression patterns of which were consistent with the sequencing data. Because DMGs such as ABCG2, mTOR, STK3, ACVR1 and PSMD7 may contain two or more DMRs, these did not show significant differences between the HP and LP groups in our study.

In our study, 11,520 of the 12,831 DMGs were enriched in three categories, as determined by GO analysis: biological processes, molecular function and cellular components. Strict conditions were followed to select the most likely DMR related genes involved in the regulation of ovarian functions. Eventually, we identified 10 eligible DMGs, compared with LP group, CTNNB1, FST, LHCGR, and TGFB3 were hypomethylated, and BMP7, BMPR1B, FSHR, TGFB2, INHBA and JUP were hypermethylated in the HP group. Moreover, the DMRs of these DMGs were all located in intron and distal intergenic regions in the genome. Recent evidence has shown that intragenic DNA methylation plays a role in the regulation of alternative splicing [30].

Key DMGs in the TGF-β superfamily

Bone morphogenetic proteins, which belong to the transforming growth factor β (TGFβ) superfamily, are known to have effects on reproduction. A mutation in the BMPR1B gene, called FecB, was the first major gene associated with prolificacy in sheep [31, 32]. Previously, it was also reported that highly prolific Booroola sheep have a mutation in the intracellular kinase domain of BMPR1B (ALK-6) which is expressed in both oocytes and granulosa cells, and is associated with hyper prolificacy of these ewes [31]. BMPR1B has an additive effect on ovulation rates and litter size in several sheep breeds [4]. BMP7 has been reported to have a significant role in ovarian folliculogenesis due to its expression from the time of committed follicles onward in rat thecal cells [33, 34], and BMP7 was shown to have the function of down-regulating StAR and progesterone production in human granulosa-lutein cells [35]. However, the information regarding how the TGFβ family alters mammalian reproduction through DNA methylation is limited. Our findings in Hu sheep also support the earlier reports, in which significant differences were found in the levels of hypermethylation BMPR1B and BMP7 genes in HP ovaries, and the pathways involving these genes were enriched in ovarian cumulus expansion and the menstrual cycle phase respectively. Moreover, using STRING analysis, we also found that BMPR1B was the hub of these reproduction related DMGs and correlated with BMP7, FSHR and LHCGR.

Key DMGs in Gonadotropin hormone

In our study, several DMGs related to hormone function, such as FSHR and FST, showed significant differences in methylation levels (FSHR and FST are up-regulated, and LHCGR is down-regulated) in ovary tissue from the HP group compared with that from the LP group, which may influence mRNA expressions of these genes. Earlier studies provided strong evidence that ovarian follicles lacking FSH or FSH receptors fail to progress to a preovulatory stage, resulting in infertility [36], and expression of the LHCGR at relatively high levels in granulosa cells was required for preovulatory follicles to respond to the midcycle surge of LH that promotes ovulation, oocyte maturation, and corpus luteum formation [37]. The relatively higher levels of expression were found for the transcripts of FSHR and LHCGR across ovaries and ovarian follicles in FecB carrier ewes [33]. FST has been considered to play an important role in ovarian development in species such as mice and pigs [38, 39]. FST secreted by granulosa cells specifically inhibits FSH biosynthesis and secretion. However, FST expression patterns show significant divergence among species; in our study, only one DMR (chr16:25,647,868–25,647,877 on the antisense strand, strong hypomethylation in the HP group) was related to FST and located in the distal intergenic region after the TSS, which indicated that the expression of FST may be influenced by this DMRs, just as intragenic DNA methylation status can down-regulate IGF2 gene expression in bovines [40]. Above all, changes in these gonadotropin receptor mRNA expression levels may determine follicular responses to gonadotropins thereby inducing the release of ovum.

Key DMGs in folliculogenesis and ovulation

In our study, several DMGs relating to folliculogenesis and ovulation were identified, including CTNNB1, INHBA and JUP. Previous studies have shown that CTNNB1 can facilitate FSH induced follicular growth and decreases follicle atresia, but that it negatively affects LH induced ovulation and luteinization in mice [41]. Moreover, it plays important roles in regulating patterning and morphogenesis that are related to adherent junctions and are required for gonadogenesis [42, 43], and similar results were also found in our study, moreover, the DMRs related to CTNNB1 were located in distal intergenic and intron regions (2nd, 4th, 5th of 36). INHBA expression was stimulated by BMP15 in granulosa cells from wild type ewes, and may play roles in the increased ovulation rates [44], and in our study, INHBA was correlated with FST, LHCGR, FSHR, BMPR1B, TGFB2 and TGFB3. Up to now, several studies have reported on the function of JUP in ovarian cancers [45,46,47,48,49], but information regarding JUP function in mammal reproduction has been limited, our study supports JUP as being related to oocyte development, as identified GO analysis, and two DMRs (chr11:41,441,802–41,442,083 on the antisense strand, strong hypomethylation in the HP group; chr11:41,458,174–41,458,384 on the antisense strand, strong hypermethylation in the HP group) were related to JUP and located in distal intergenic and intron regions respectively.

In summary, this study provides a comprehensive analysis of the DNA methylation profiles of Hu sheep ovaries for HP and LP ewes. We identified DMRs and genes associated with these regions. Pathway and network analyses of these DMRs revealed several candidate genes that may affect ovarian function including gonadotropin, folliculogenesis and ovulation. We will validate those DMR-related genes from this study in different stages of follicles development in the future. The results of this study might therefore provide novel clues for deciphering the epigenetic mechanisms of sheep ovarian function and will likely contribute to improving reproductive capacity.

Conclusion

This study revealed the global DNA methylation patterns of sheep ovaries associated with high and low prolificacy. We explained the differences in genomic DNA methylation between HP and LP sheep, and we observed that several DMRs/DMGs were most likely related to changes in Hu sheep prolificacy. Our results demonstrate that DNA methylation may contribute to a better understanding of epigenetic regulation in sheep reproductive capacity.