A multi-scale analysis of bull sperm methylome revealed both species peculiarities and conserved tissue-specific features

Animals and cell/tissue collection

All study methods were implemented in accordance with EU guidelines and regulations (Directive 2010/63/UE). For animals maintained in INRA facilities, the experimental protocols were approved by the INRA local Ethics Committee (COMETHEA, authorization numbers 12/160 and Méjusseaumes Animal Care committee 0162503). The bull samples originated from bulls selected for artificial insemination and were provided by commercial companies: Montbéliarde breed by GEN’IATEST (France) and UMOTEST (France), Holstein and Normande breeds by EVOLUTION (France) and Belgian White Blue breed by AWE (Belgium). Other bovine tissues were collected from Holstein cows maintained at the INRA experimental farm (UCEA, INRA, France). The ram, goat and boar semen and blood samples were supplied by commercial companies (OSON, Capgenes and LNCR, respectively, France). Mice semen samples were collected from the caudal epididymis of 7-week-old male C57Bl/6JOlaHsd mice supplied by Harlan Laboratory (Netherlands) and euthanized by cervical dislocation. Stallion semen was supplied by Dr. M. Magistrini (UMR INRA 0085 PRC, France). Human sperm samples originated from patients included in a PHRC METASPERME study, coordinated by Dr. R. Levy (Laboratoire d’Histologie Embryologie Cytogénétique CECOS, Hôpital Jean Verdier, France); this study received the approval from French ethics board (Conseil d’évaluation éthique pour les recherches en santé, CERES) and all the patients gave their informed written consent to participate.

In Fig. 1c-d, semen and blood from the same individuals were collected using standard procedures on bulls and rams maintained in semen production centers. Total blood was used for DNA extraction. For both bulls and rams, collected semen was extended with Optidyl (Cryo-Vet) and either underwent direct DNA extraction (fresh semen) or was subjected to standard techniques for semen processing (straw conditioning, freezing and storage in liquid nitrogen; frozen semen). Other bull semen samples were in the form of frozen straws stored in liquid nitrogen.

For all semen samples except those from humans and stallions, the possible contamination of spermatozoa by somatic cells was checked systematically under the microscope and confirmed to be below detectable levels. The human and stallion semen samples contained observable somatic cells and were therefore processed as previously described to ensure the absence of potential contamination: the human semen samples were subjected to a stringent somatic cell lysis protocol [18] and the sperm from stallions were purified by single layer centrifugation using Androcoll-E-Large (SLU, Uppsala, Sweden) [31].

Primary cultures of fibroblasts were derived from ear skin biopsies from three separate adult heifers and cultured until passage 11 in Dulbecco’s modified Eagle medium supplemented with 10% fetal calf serum and 1% penicillin-streptavidin (Life Technologies) at 38 °C with 5% CO₂. Livers were obtained from adult cows slaughtered at the INRA experimental facilities. Peripheral blood mononuclear cells (PBMCs) were isolated from blood collected from the jugular vein and centrifuged using a Ficoll gradient. To obtain the monocyte fraction, PBMCs were incubated in the presence of microbeads conjugated to monoclonal anti-CD14 antibodies (mouse IgG2a; Miltenyi Biotec) in MACS BSA buffer, for 15 min. at 4 °C under gentle agitation. Magnetic separation was then performed using MS Columns following the manufacturer’s instructions. The tissues and cells were snap-frozen in liquid nitrogen and stored at − 80 °C until DNA extraction.

Genomic DNA extraction and genotyping

One straw of bull semen was used for DNA extraction (about 20 million spermatozoa). After thawing at 37 °C, the semen was washed with phosphate buffer saline (PBS) to remove the extender, and incubated overnight at 55 °C in 200 μl lysis buffer (10 mM Tris-HCl pH 7.5, 25 mM EDTA, 1% SDS, 75 mM NaCl, 50 mM dithiothreitol (DTT) and 0.5 μg glycogen) in the presence of 0.2 mg/ml proteinase K. After incubation with 25 μg/ml RNAse A for 1 h at 37 °C, genomic DNA was extracted twice using phenol and phenol:chloroform (1:1), then ethanol precipitated and washed. The dried pellet was re-suspended in TE buffer (10 mM Tris HCl pH 7.5, 2 mM EDTA) and the DNA concentration was measured using a Qubit 2.0 Fluorometer (Invitrogen). Fresh bull semen, total blood from bulls and rams, fresh or frozen semen from rams, boars, mice, goats, processed semen from stallions and men were treated in an identical manner as the frozen straws from bull semen.

DNA extraction from liver samples was performed as described elsewhere [32]. The same procedure was used for fibroblasts and monocytes, except that the cells were lysed by the direct addition of lysis buffer and proteinase K to the cell pellet.

The genotyping of two Montbéliarde bulls was performed by LABOGENA (France) on semen and blood DNA from the same individuals using the commercially available BovineSNP50 v2 BeadChip (Illumina). Genotypes were determined using the Genotyping Module of GenomeStudio software (Illumina). For each animal, copy number variations (CNVs) were searched for and compared between tissues. The Log R Ratio (LRR, normalized measurement of total signal intensity) and B Allele Frequency (BAF, measurement of the allelic intensity ratio) were used to infer copy number changes in the genome. For example, in the presence of a deletion, LRR values increase and BAF values cluster around 0 or 1 but are absent at around 0.5, due to a lack of heterozygotes. LRR and BAF were then plotted along the genome and compared between tissues.

In silico analyses

A script developed in house was used to extract the coordinates of all the CCGG sites present in the genomes of different species (cattle, sheep, horse, pig, mouse and human; Fig. 1e). The CCGG sites were then annotated relative to different gene features, CpG density and repeats by means of a pipeline developed in house (https://github.com/FAANG/faang-methylation/tree/master/RRBS-toolkit/Annotation) and using the genome annotation features indicated in Additional file 1: Table S1. The following criteria were applied: promoter-TSS, − 1000 to + 100 bp relative to the transcription start site (TSS); TTS: -100 to + 1000 bp relative to the transcription termination site (TTS); shore, up to 2000 bp from a CpG island (CGI); and shelve up to 2000 bp from a shore. A site/fragment was considered to belong to a CGI (respective shore and shelve) if an overlap of at least 75% was observed between the site/fragment and the CGI (respective shore and shelve). A site/fragment was considered as being overlapped by a repetitive element whatever the extent of this overlapping.

In Additional file 1: Tables S2 and S3, in silico reduced representation (RR) genomes digested by MspI restriction enzyme were produced for different species and using different size selection criteria by means of another pipeline developed in house (https://github.com/FAANG/faang-methylation/tree/master/RRBS-toolkit/RR_genome). The RR genome fragments were then annotated as explained above.

Luminometric methylation assay (LUMA)

The conditions were compared using non-parametric tests suited to small samples (permutation tests for two independent samples or for two paired samples according to the situation, with Monte-Carlo sampling of 100,000 permutations) or using t-test and analysis of variance when appropriated (larger samples with a normal distribution).

Methylated DNA immunoprecipitation (MeDIP), microarray hybridization and data analysis

MeDIP and quality controls by PCR were performed as described elsewhere [32]. The antibody used for immunoprecipitation was BI-MECY-1000 5-methylcytidine antibody (Eurogentec). To prevent any technical bias, the products of five independent MeDIP experiments were pooled for each sample. After moderate genome amplification, the pooled MeDIP reactions and corresponding input DNA were labelled with Cy3 and Cy5 and hybridized on a Roche-NimbleGen 3x720K microarray, with technical dye-swaps for every sample. The microarray targeted the promoter region (− 2000 to + 1360 bp relative to the gene start) of 21,296 bovine genes, according to an annotation file downloaded from the Johns Hopkins University Center for Computational Biology FTP website (ftp://ftp.ccb.jhu.edu/pub/data/assembly/Bos_taurus/Bos_taurus_UMD_3.0/annotation/; accessed Aug. 2010). The microarray design and hybridization protocol, as well as more details on the data analysis, can be found in [35]. The identification of regions of interest containing clusters of probes enriched in at least one tissue, the identification of differentially methylated regions (DMRs) among these regions of interest, as well as the calculation of mean percentages of enriched probes (Pr) in each tissue for each region r (DMR or region of interest), are detailed in the Additional file 1: Supplementary methods. Three Pr were obtained per region: Pr_sperm, Pr_liver and Pr_fibroblasts. The scatterplot shown in Fig. 3d illustrates the Pr_liver - Pr_sperm and Pr_fibroblasts - Pr_sperm differences for the regions of interest (in black) and for the DMRs specific to the comparison between sperm and somatic cells (in red). Positive values for both Pr_liver - Pr_sperm and Pr_fibroblasts - Pr_sperm indicated that the two somatic cell types were more methylated than sperm in the region considered, while negative values for both Pr_liver - Pr_sperm and Pr_fibroblasts - Pr_sperm indicated the contrary. Similarly, an opposite sign for the values of Pr_liver - Pr_sperm and Pr_fibroblasts - Pr_sperm reflected the contrasting behaviors of liver and fibroblasts. Among the DMRs, those undermethylated in sperm were selected based on a positive value for both Pr_liver - Pr_sperm and Pr_fibroblasts - Pr_sperm and were annotated as explained regarding the in silico analyses, together with the 27,684 regions of interest. Genes containing DMRs were subjected to DAVID analysis (Database for Annotation, Visualization and Integrated Discovery; [36]) using genes containing the 27,684 regions of interest as the background.

Reduced representation bisulfite sequencing (RRBS) and data analysis

RRBS libraries were prepared as described elsewhere [37, 38]. Briefly, 200 ng of genomic DNA were digested by MspI (Thermo Scientific), end-repaired and ligated to 55 bp Illumina adapters for paired-end sequencing. Size selection by gel excision was performed in order to select fragments ranging from 150 to 400 bp (genomic fragments of 40-290 bp + adapters). The DNA was then purified using the MinElute gel extraction kit (Qiagen) and then bisulfite-converted twice consecutively with the EpiTect bisulfite kit (Qiagen), following the manufacturer’s instructions for DNA extracted from FFPE tissues. Converted DNA was amplified with Pfu Turbo Cx hotstart DNA polymerase (Agilent) using 14 PCR cycles for sperm and fibroblasts and 12 cycles for monocytes. The libraries were then purified using Agencourt Ampure beads (Beckman-Coulter) and sequenced on an Illumina HiSeq2500 sequencer to produce 75 bp paired-end reads (Integragen SA, France).

RRBS sequences were analyzed using an integrated pipeline combining scripts developed in house in Python, R and Shell, together with external tools (https://github.com/FAANG/faang-methylation/tree/master/RRBS-toolkit/). Details about the analysis and the identification of differentially methylated cytosines (DMCs) are provided in the Additional file 1: Supplementary methods.

For each tissue, the mean methylation percentage was calculated (mean of the methylation percentages obtained in the two biological replicates) as well as the difference between two tissues (Additional file 2, column L). The scatterplot shown in Fig. 3e illustrates the differences between monocytes and sperm and between fibroblasts and sperm for the 1,580,644 CpGs covered by 5 to 500 uniquely mapped reads (CpGs 5-500) in all six samples (in black) and for the DMCs specific to the comparison between sperm and somatic cells (in red). Among these DMCs, those undermethylated in sperm were selected based on a positive value for both differences. Together with the CpGs 5-500, they were then annotated relative to gene features, CGIs and repeats as explained for the in silico analyses. Genes containing DMCs were subjected to DAVID analysis using genes containing the 1,580,644 CpGs 5-500 as the background.

To better characterize repetitive elements, an artificial genome containing the consensus sequence of each bovine repeat was constituted from the Repbase database [39]. Reads were aligned on this artificial genome as explained above, and the average methylation percentage was calculated for each repeat and each sample (average methylation percentage for all CpGs included in one genomic repeat and covered by either 5-500 reads or by > 500 reads).

Bisulfite-pyrosequencing

Bisulfite conversion was performed on 1 μg genomic DNA as described elsewhere [40]. After ethanol precipitation, the DNA pellet was re-suspended in 20 μl H₂O.

For LSM4 and BTSAT4, primers were designed using the MethPrimer program [41] and amplifications were carried out from 1 μl treated DNA with Platinum Taq DNA polymerase (Invitrogen), according to the manufacturer’s instructions with variable MgCl₂ concentrations. The following program was used: 3 min. at 94 °C followed by 50 cycles of 30 s. at 94 °C, 1 min. at variable hybridization temperatures, 1 min. at 72 °C, and finally 10 min. at 72 °C. For DDX4 and SYCP3, primers were designed using the Pyromark assay design software (Qiagen) and amplifications were performed using the Pyromark PCR kit (Qiagen) according to the manufacturer’s instructions. The primers used to amplify each region are listed in Additional file 1: Table S4, together with the hybridization temperatures and MgCl₂ concentrations. The reverse primers were 5′-biotinylated.

After denaturation and purification, the biotinylated antisense strand of PCR product was used as a template for pyrosequencing with 0.3 μM pyrosequencing primer, using the Pyromark Q24 device and Pyromark Gold Q96 reagents (Qiagen). The pyrosequencing primers are listed in Additional file 1: Table S5. Each CpG was assayed in duplicate, and inconsistent duplicates (more than 5% difference) were repeated. The methylation percentage per CpG was then obtained by calculating the mean of all replicates that passed quality control by the Pyromark Q24 software. The statistical analysis was performed on the mean percentage per CpG using permutation tests as explained for LUMA.

Global DNA methylation level is low in bull sperm

We first assessed the global level of DNA methylation in bull sperm relative to somatic cells (PBMCs) using LUMA [33, 34] on a large sample size. The average methylation at CCGG sites was dramatically lower in sperm (45.5%, n = 185) than in PBMCs (74.8%, n = 73; Fig. 1a). The standard deviations in this large sample were weak (2.6% for sperm and 1.6% for PBMCs), demonstrating the limited inter-individual variability within each cell type and the reliability of the technique to assess global DNA methylation. Because the sperm samples were collected from Holstein, Montbéliarde, Normande and Belgian White Blue bulls, we investigated the effect of the breed on global sperm DNA methylation. Methylation was significantly lower in Belgian White Blue than in any other breed and significantly lower in Holstein than in Normande and Montbéliarde (Fig. 1b). A bootstrap analysis confirmed that these breed-related differences were not due to an unbalanced number of bulls from each breed (Additional file 1: Figure S1). These results demonstrated that global methylation varied across different bovine breeds, suggesting that the presence of DNA polymorphism could influence the global CCGG content and methylation. However, the range of variation (from 42.6% in Belgian White Blue to 47.1% in Normande) was weak relative to the 30% difference we observed between sperm and PBMCs, suggesting that sperm weaker methylation was not breed-dependent. To determine whether the global DNA methylation of sperm was comparably low in another species, we collected paired semen and blood samples from bulls and rams. While sperm was less methylated than blood in both species (Fig. 1c), the difference between the two cell types was much greater for bulls (30% less methylation in sperm than in blood) than for rams (only 10% less methylation in sperm than in blood). Differences between species were observed in both cell types, but were broader for sperm (20% difference) than for blood (only 1.2% difference). In sperm, the difference between species was independent of the cryopreservation process (Fig. 1d).

This lower methylation of sperm may have resulted from a biased representation of the regions present in our sperm genomic DNA. Indeed, the high level of sperm chromatin compaction impeded the complete extraction of genomic DNA using standard procedures, which were therefore optimized by the addition of reducing agents such as DTT (see for instance [42]). The protocol we used during this study performed well for DNA extraction from human sperm (50 mM DTT, [7]). To investigate whether these conditions were also suited to bovine sperm, we used genomic DNA from two paired sperm and blood samples (the same samples as in Fig. 1c), both of which were purified in the presence of 50 mM DTT, as a template for genotyping. As expected, the same genotype was obtained from blood and sperm DNA in each bull. Likewise, for each animal, CNV profiles were compared between tissues because any difference might be indicative of preferential extraction. The plots were similar and no gross discrepancies could be observed between tissues (Additional file 1: Figure S2). This result therefore ruled out the possibility that some specific regions of the bovine genome failed to be extracted from sperm chromatin under our experimental conditions.

Because the genomic distribution of CCGG sites may display species-specific variations that might influence the methylation results, we compared sperm CCGG methylation in several species relative to the genomic features available. In the horse, pig, mouse, human and goat, sperm DNA methylation level was much higher than that measured in bovine and closer to that determined in sheep (Fig. 1e). The CCGG distribution in the bovine genome was within the range observed in other species for most of the features examined, except for CGIs that accounted for 14.2% of total CCGG sites in bovine vs. 8.0 to 12.2% in other species. The important representation of bovine CCGGs in CGIs was balanced by a low representation in open sea (regions with a low CpG density). Most mammalian CGIs are unmethylated in somatic cells and can become methylated during development and disease [43]. Because blood cells (in which most CGIs are supposed to be unmethylated) displayed a roughly similar methylation level in bovine and sheep, it is unlikely that the lower level of DNA methylation in bull sperm resulted solely from the higher percentage of bovine CCGG sites in CGIs. Another remarkable feature of the bovine CCGG sites was their strong enrichment in satellites, which represented 1.6% of all CCGGs vs. 0.1 to 0.6% in other species. However, whatever the species examined, the < 2% difference in CCGG sites present in satellites could not account for the > 10% difference in sperm methylation (Discussion).

Taken together, these data demonstrated that compared to somatic cells, bull sperm displayed a dramatically lower level of CCGG methylation which seemed to be specific to the bovine species. This lower methylation was neither related to the process of semen cryopreservation nor to a technical artefact, and could not be fully explained by the genomic distribution of CCGG sites in the bovine species.

The tissue/cell type is a major determinant of DNA methylation landscapes in cattle

Because high-throughput analyses were necessary to identify regions that were hypomethylated in bovine sperm, we decided to assess two cost-effective approaches which are widely used to study DNA methylation: MeDIP-chip [44] and RRBS [37]. MeDIP-chip enables precise targeting of specific regions in the genome through custom design of the microarray (in our case, 3360 bp spanning the promoter and upstream region of each of the 21,296 bovine genes [35]), while RRBS offers a base-resolution analysis of CpG-rich regions through the combined use of enzymatic digestion (MspI) and the size selection of restriction fragments. Because no data were available regarding the optimal size window for RRBS in cattle, we conducted an in silico prediction of the genome coverage by RRBS and compared the results with those of MeDIP-chip (Additional file 1: Supplementary methods, Tables S2-S3). These in silico analyses suggested that the RRBS procedure could be successfully adapted to the bovine genome using MspI with a size window of 40-290 bp, and would lead to a base-resolution map of the methylome with coverage that would complement that of the MeDIP-chip.

Independent bovine samples were analyzed using MeDIP-chip (sperm, liver and fibroblasts) and RRBS (sperm, monocytes and fibroblasts) in order to determine cell type-dependent variations of the methylome using these two complementary technologies. For RRBS we selected CpGs covered by 5 to 500 uniquely mapped reads for each sample (CpGs 5-500), from which an average methylation rate was calculated. The average methylation in sperm (51.8%) was higher than in fibroblasts (48%) and lower than in monocytes (57.6%; Additional file 1: Table S6). We next categorized the CpGs 5-500 into hypo- (< 20% methylation), intermediate (20–80% methylation) and hypermethylated CpGs (> 80% methylation) and observed a larger proportion of hypo- and a smaller proportion of intermediate CpGs in sperm than in somatic cells, which was counterbalanced by a large proportion of hypermethylated CpGs (Fig. 2a). This bimodal distribution of methylation in sperm probably explained the intermediate level of average methylation for CpGs 5-500. It was noted that when only CpGs covered by more than 500 uniquely mapped reads were considered (CpGs > 500), average methylation increased in monocytes (81.4%) and in fibroblasts (66.6%), but fell dramatically in sperm (22.5%; Additional file 1: Table S6).

We next conducted descriptive analyses of the MeDIP-chip and RRBS data. For MeDIP-chip, a normalized factor NEpi, representing the number of enriched probes at promoter p for sample i, was calculated for all promoters and samples, and hierarchical clustering and principal component analysis (PCA) were run on the resulting matrix. For RRBS, PCA and hierarchical clustering were computed from the matrix of methylation percentages obtained for each CpG 5-500 and each sample (Fig. 2b). With both types of descriptive analysis, the samples were clearly clustered according to the tissue/cell type using both MeDIP-chip and RRBS. Interestingly, hierarchical clustering revealed that the distance between sperm and other cell types was more important than the distance between liver and fibroblasts or monocytes and fibroblasts, highlighting the methylation specificities of germinal cells compared to somatic cells. This could also be seen in PCA, where dimension 1 opposed sperm to one or both somatic cells/tissue, while the difference between the two somatic cell types was more apparent along dimension 2.

Our results therefore showed that the tissue/cell type represented the main source of variation in methylation, and that sperm-specific methylation profiles could emerge from our MeDIP and RRBS data.

Identification of regions and CpGs hypomethylated in bull sperm

We next ran a differential analysis on each pair of tissue/cell types. For MeDIP, we compared the proportion of enriched probes between two tissues/cell types in 27,684 regions of interest containing clusters of probes enriched in at least one tissue/cell type (see Additional file 1: Supplementary methods). This led to the identification of 4329 DMRs between sperm and liver; 3780 DMRs between sperm and fibroblasts and 2803 DMRs between fibroblasts and liver. The features of each DMR and the corresponding promoter regions are summarized in Additional file 3. For RRBS, we identified 298,901 DMCs between monocytes and sperm; 450,971 DMCs between fibroblasts and sperm, and 239,036 DMCs between monocytes and fibroblasts using stringent criteria (Additional file 1: Supplementary methods, Table S7 and Additional file 2).

Consistent with the results of clustering and PCA, the number of DMRs/DMCs was higher between sperm and any somatic cell type/tissue than between two somatic cell types, for both MeDIP and RRBS. Fig. 3a-b shows the distribution of DMRs and DMCs in different Venn territories corresponding to pairwise comparisons under MeDIP (Fig. 3a) and RRBS (Fig. 3b). DMRs and DMCs specific to the comparison between sperm and somatic cells could be deduced from these territories (in red). These 1678 DMRs and 174,103 DMCs were located at the intersection between the “sperm vs. somatic cell type 1” and “sperm vs. somatic cell type 2” territories (yellow and orange) excluding the DMRs/DMCs also shared by the “somatic cell type 1 vs. somatic cell type 2” comparison (three-color territory). Because the comparison between sperm and fibroblasts was performed using both MeDIP and RRBS, the actual complementarity of the two technologies could be assessed using real data. As shown in Fig. 3c, only a limited subset of CpGs was shared by DMRs identified using MeDIP and DMCs identified using RRBS. This result, together with the in silico analysis indicating that the targeted regions were largely different under MeDIP-chip and RRBS, clearly demonstrated that we were able to identify two distinct subsets of sperm-specific DMRs/DMCs.

We next investigated whether these two subsets of sperm-specific DMRs/DMCs might display similar features in terms of variations in methylation. For both MeDIP (Fig. 3d) and RRBS (Fig. 3e), we plotted the differences in methylation between each somatic cell type and sperm for sperm-specific DMRs/DMCs. The appropriate background was used for each technology: the 27,684 regions of interest subjected to differential analysis in MeDIP and the 1,580,644 CpGs 5-500 covered in all six samples analyzed using RRBS. With both MeDIP and RRBS, background regions/CpGs (in black) were particularly concentrated around the center of the plot, illustrating that most of the regions/CpGs analyzed did not display cell type-dependent variations. By contrast, sperm-specific DMRs/DMCs were located along a diagonal that ran from the lower left-hand corner to the upper right-hand corner, meaning that these DMRs/DMCs behaved similarly in the two somatic cell types when compared to sperm. Most strikingly, a great majority of sperm-specific DMRs/DMCs (68% for MeDIP and 79% for RRBS) were grouped in the upper right-hand corner, indicating that they were hypomethylated in sperm when compared to the two somatic cell types examined.

Taken together, these results demonstrated that using two complementary technologies we were able to identify two distinct subsets of DMRs/DMCs specific to the comparison between sperm and somatic cells, most of them being less methylated in sperm. These hypomethylated sperm-specific DMRs/DMCs (hypo-DMRs/DMCs) could therefore partly explain the lower global DNA methylation of sperm observed in LUMA which was particularly marked in cattle.

Hypomethylation in bull sperm targets specific genomic features and functions

To determine whether specific gene ontology (GO) terms were enriched in the sperm hypomethylated regions, we next annotated the hypo-DMRs identified by MeDIP relative to genes. Consistent with the microarray design, most of the 1144 hypo-DMRs were located in or close to genes according to the criteria described in the Methods, resulting in a list of 701 unique genes which were then subjected to DAVID analysis. Significant enrichments were found for biological processes such as sexual reproduction (36 genes), fertilization (15 genes) and RNA transport (11 genes). Further analysis using a more restrictive list of GO terms led to the identification of functional clusters involved in mRNA processing (Fig. 4a) and meiosis/spermatogenesis (Fig. 4b).

Because CpGs targeted by RRBS are scattered along the genome, we then started to characterize the hypo-DMCs identified by RRBS relative to different genomic features (genes, CpG density and overlapping repeats; Fig. 5a). The most remarkable observation was a dramatic enrichment of hypo-DMCs for repeats (24.5% vs. 13.2% in background; right panel), and particularly for satellites (64.7% vs. 17.8%). In order to get a more precise picture of the methylation status of repetitive elements by rescuing some of the information included in the ambiguous reads, we aligned the reads on a Repbase artificial genome containing the consensus sequence of each bovine repeat (see Additional file 1: Supplementary methods for details, and Additional file 4 and Additional file 1: Table S8 for data). The hypomethylation of sperm was clear in satellites and also in rDNA repeats encoding ribosomal RNAs (Additional file 1: Figure S3). We therefore concluded from both CpGs targeted by uniquely mapped reads (that had been used to identify hypo-DMC) and CpGs targeted by ambiguous reads, that specific families of repetitive elements were massively hypomethylated in bovine sperm.

Regarding gene features with respect to RRBS data, the proportion of genes containing hypo-DMCs was relatively unchanged compared with background (Fig. 5a, left panel), but interestingly, exons were more frequently represented (26.8% vs. 13.8% in background) while promoter-TSS were represented less than in background (4.8% vs. 13.1%). The distribution of hypo-DMCs in CGIs, shores and shelves was identical to that observed with the background (middle panel). We investigated whether specific GO terms were enriched in gene features displaying a different representation in hypo-DMCs and background. For exonic hypo-DMCs, which accounted for 2713 unique genes, significant enrichments were found for biological processes such as cell adhesion (213 genes), the regulation of signaling (382 genes) and cell migration (177 genes), and one main functional cluster related to cell adhesion could also be identified (Fig. 5b). For hypo-DMCs located in promoter-TSS, which accounted for 1200 unique genes, the most enriched biological process was sexual reproduction (44 genes), and two functional clusters were identified (Fig. 5c) as being related to piRNA metabolism (left panel) and to meiosis and spermatogenesis (right panel).

From an analysis of both hypo-DMRs obtained in MeDIP and hypo-DMCs identified through RRBS, we therefore concluded that undermethylation in sperm essentially targeted repeats and the promoters of genes important to spermatogenesis (which is the differentiation process that eventually leads to the mature sperm we analyzed), but also to genes involved in cell communication, signaling and migration that may be essential to both sperm functions and post-fertilization steps.

Hypomethylation of four regions is confirmed by bisulfite-pyrosequencing

Four regions were selected for validation, based on their hypomethylation in sperm, the presence of both DMRs obtained by MeDIP analysis and DMCs obtained by RRBS analysis whenever possible, and their position relative to genes involved in sperm functions. Figure 6a shows the detailed localization of these regions together with their coverage and individual methylation in both MeDIP and RRBS. LSM4, which contained a hypo-DMR, is a gene involved in RNA processing [45]. Genes SYCP3 and DDX4, which contained both hypo-DMRs and hypo-DMCs identified in our study, play a major role in spermatogenesis insofar as either mutation or aberrant methylation of these genes associate to male infertility [46, 47]. BTSAT4 was the most frequently represented bovine satellite in our RRBS data and displayed undermethylation in sperm (Additional file 4). Overall, the four regions represented 42 analyzed CpGs. We used the pyrosequencing of bisulfite-converted DNA [40] to quantify the absolute methylation percentage of individual CpGs in sperm, liver and fibroblasts (Fig. 6b; genomic DNA from monocytes was in limited amounts and was saved for separate investigations). Consistent with the MeDIP data showing enriched probes along the LSM4 promoter in liver and fibroblasts but not in sperm, the CpGs analyzed by pyrosequencing were all hypomethylated in sperm. In this region, the CpGs assessed by pyrosequencing were not covered by RRBS, but those in DDX4, SYCP3 and BTSAT4 were covered by the three techniques and displayed hypomethylation in sperm whatever the technique used. By pyrosequencing four additional CpGs, we also checked that one region specifically hypermethylated in sperm compared to somatic cells validated (Additional file 1: Figure S4).

In conclusion, the results obtained using MeDIP, RRBS and pyrosequencing were in excellent agreement, which validated the high-throughput data and led to the characterization of gene regions with methylation patterns specific to bull sperm.