Comparative analysis of sperm DNA methylation supports evolutionary acquired epigenetic plasticity for organ speciation

Methylation of DNA is a well-studied epigenetic mark that is evolutionary conserved among most eukaryotic organisms such as fungi, plants and animals [1,2]. It involves adding a methyl group to a cytosine base within the DNA molecule and is associated with a number of key processes, including regulation of gene expression [3], genomic imprinting [4], X-chromosome inactivation [5], preservation of chromosome stability [6], repression of transposable elements [7,8], aging [9,10] and gametogenesis [11]. Aberrant DNA methylation is observed in a plethora of pathological conditions, such as carcinogenesis [12], kidney disease [13] and metabolic disease [14,15], supporting the notion that DNA methylation also participates in the risk to develop disease.

The role of DNA methylation in gametogenesis is of particular importance: methylation marks in the gametes establish the developmental programming, thereby influencing cell fate, organ formation and the nature of the response to environmental factors later in life [11]. DNA methylation patterns are substantially lost during gametogenesis and subsequently restored at the fully differentiated gamete stage [16]. While such epigenetic reprogramming is conserved in mammalian species [16,17], alterations may influence gene expression and the predisposition to disease [18]. Epigenetic modifications of sperm cells in different phases of spermatogenesis take part in different aspects of embryonic development upon fertilization [19]. Environmental factors like dietary habits, physical activity and pollutants modulate the DNA methylation signature of spermatozoa and the phenotype of the offspring [20–22], indicating that epigenetic variation in sperm may contribute to environmentally driven phenotypic plasticity in the next generation of offspring.

While several studies have mapped the epigenetic patterns of spermatozoa from humans and rodents separately [23–25], little is known about the comparative differences across these species. In the present study we compared the DNA methylomes of spermatozoa from humans, mice, rats and mini-pigs, with a specific interest in conserved genomic regions. Species-specific epigenetic signatures were analyzed for their potential role during early embryo development by investigating the expression of genes near epigenetically variable regions in human and mouse embryos at the preimplantation stage. Our analysis identified epigenetic variation near genes referenced as participating in brain function and the perception of smell, supporting species-specific means to develop certain physiological functions.

Materials & methods

Biological samples & bisulfite sequencing

We performed our analyses on a panel of publicly available data sets, generated for different studies. Human (Gene Expression Omnibus [GEO] accession: GSE49624, GSE109478), mouse (GSE49624, GSE174093) and rat (GSE79566, GSE173898) sperm [20,25–28] are included in whole-genome bisulfite sequencing (WGBS) and reduced-representation bisulfite sequencing (RRBS) data sets, respectively. Samples were only taken from the control groups of the studies in which an exposure/intervention was conducted.

For RRBS on mini-pig sperm, semen samples were collected from four 8-month-old genetically inbred 100% miniature pigs from Specipig^® (Pioneer breeding and biomedical research centre, Barcelona, Spain). Genomic DNA was extracted from the sperm cells using an AllPrep DNA/RNA/miRNA Universal Kit (Qiagen, Hilden, Germany; cat. no. 802224). The procedure was changed for processing of sperm according to the recommendations of the manufacturer. RRBS libraries were constructed from 100 ng genomic input DNA using the Ovation^® RRBS Methyl-Seq library preparation kit (Tecan, Männedorf, Switzerland; cat no. 0553-32). Briefly, genomic DNA was digested with MspI enzyme, followed by adaptor ligation to Ovation RRBS Methyl-Seq adaptors containing unique molecular barcodes, followed by a final end-repair step. The ligated DNA was subjected to bisulfite conversion and left in the thermal cycler overnight at room temperature. The following day, the bisulfite-converted DNA was amplified by PCR followed by a final clean-up step. Sequencing of RRBS libraries was performed on a NovaSeq 6000 platform (Illumina, CA, USA) as 101-bp single-end reads. The results are available in the NCBI GEO database, with the accession number GSE213548.

For human WGBS data, samples were obtained from men attending the University of Utah Andrology laboratory consented for research (Institutional Review Board approved protocol # 7790). For mouse WGBS data, the study was done under an approved Institutional Animal Care and Use Committee protocol at University of Utah School of Medicine, UT, USA. For rat WGBS data, part of the study was approved by the Zoos scientific review committee and also by the Institutional Animal Care and Use Committee at University of Southern California, CA, USA. For human RRBS data, the study was approved by the Ethics Committee from the Capital Region of Denmark (reference H-1-2013-064) and informed consent was obtained from all participants. For mouse RRBS data, all procedures followed the guidelines of the National Institutes of Health Guide for the Care and Use of Laboratory Animals and the approval for the study was received from the Institutional Animal Care and Use Committee at University of Massachusetts, Amherst, USA. For rat RRBS data, the study was performed following the National Institutes of Health guidelines for the care and use of laboratory animals and approved by the Institutional Animal Care and Use Committee at The Lundquist Institute for Biomedical Innovation at Harbor–UCLA Medical Center, CA, USA. For mini-pig data, all animal procedures were performed following the National Institutes of Health guidelines for the care and use of laboratory animals at Bionea Lab and approved by the French Ministry of Education and Research (reference APAFIS#22853-201911212455784v1).

Processing bisulfite sequencing data sets

Trim Galore! v.0.4.3 (www.bioinformatics.babraham.ac.uk/projects/trim_galore/) was used to trim adapter sequences incorporated in the bisulfite sequencing reads with default parameters, except that we applied a stringency of 3 for both WGBS and RRBS data sets with the –rrbs flag only for the RRBS data set. We eliminated reads shorter than 15 bp after trimming. MultiQC v.0.8 [29] was used for quality control. The trimmed reads were aligned to the human (hg38), mouse (mm10), rat (rn6) or pig (susScr11) genomes, using the Bismark aligner v.0.18.1 [30]. We applied the Bismark deduplication script to the WGBS data set to remove reads mapped to identical starting positions [31]. RRBS reads from mini-pigs were deduplicated using the deduplicate_bismark tool with the –barcode flag because they contained unique molecular identifiers. The bismark_methylation_extractor tool was used to estimate the methylation level at each particular cytosine site. SNPs from the dbSNP database [32] (hg38 [build 151]; mm10 [build 142]; rn6 [build 149]; susScr11[build 150]) and CpGs with coverage higher than ten-times the 95th percentile of coverage were filtered out.

Orthologous CpGs were assigned using the liftOver tool from the University of California Santa Cruz (UCSC) genome browser [33] using the recently published Perl scripts [34] and the following procedure. The mouse, rat and mini-pig CpG sites were mapped to human coordinates (hg38) using chain files available in ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/liftOver/ and then mapped back to the first genome to ensure that it had a unique correspondent in both genomes. CpGs with a minimum of six counts were considered for further analysis. Multidimensional scaling (MDS) analysis was performed using the ratio of methylated to unmethylated CpGs (M value). The methylation percentages for different genomic regions, including gene bodies, promoters, exons and introns as well as global DNA methylation were obtained by aggregating the individual CpG methylation levels using the R package genomation v.1.20 [35]. Given that the assumptions of parametric testing were violated, the Scheirer–Ray–Hare test (similar to a two-way analysis of variance setup) was performed as previously described [28]. Subsequently, the Dunn post hoc test was performed to identify pairwise comparisons. The R package VennDiagram v.1.7.3 [36] was used to plot Venn diagrams, and Fisher’s exact test was carried out to determine the significance of the overlaps in a pairwise fashion. The density profiles of genes and gene features methylation were performed using the sm.density.compare function of the R package sm v.2.2-5.7.1 [37], applying 10,000 bootstrap for testing significant difference.

Differentially methylated gene components

The comparisons between the human, mouse, rat and mini-pig data sets were carried out using the R package GeneDMRs v.1 [38]. RefGene database hg38 and CpG island annotation from the UCSC genome browser (http://genome.ucsc.edu) were used. We obtained weighted methylation levels for each group, Q-values with the false discovery rate method [39] and methylation differences between species for all promoters, exons and introns. The differentially methylated regions (DMRs), including promoters (DMPs), exons (DMEs) and introns (DMIs), were defined as those with Q-values of less than 0.05.

Functional annotations of cluster networks of all DMRs & biological enrichment of DMP-related genes

To decipher the mechanism of action of DMPs, DMEs and DMIs identified from WGBS and RRBS data in each comparison, we made use of the ClueGO plug-in of Cytoscape v.3.9 [40], which illustrates functionally grouped networks of the nonredundant biological terms for large gene clusters. Enrichment analyses were performed using 20 genes for WGBS and ten genes for RRBS data sets in the cluster with gene ontology (GO) tree interval range of 6–14, a κ-score of 0.4 for network connectivity, two-sided hypergeometric test with p < 0.05 and Bonferroni correction. GO biological processes (GO:BP) of DMP related genes were determined with the Database for Annotation, Visualization and Integrated Discovery website [41] and visualized by the networkD3 v.0.4 [42] and GOplot v.1.0.2 R packages [43].

Identification of hypomethylated regions & visualization of methylation profiles

We additionally used MethPipe v.3.4.3 [44] to identify hypomethylated regions (HMRs) in each species using the WGBS data sets. Briefly, the alignment BAM files from Bismark were converted to the standardized SAM format using the MethPipe ‘format reads’ program. The MethPipe ‘methcounts’ function was then used to count the methylation level at each cytosine site. The ‘merge-methcounts’ program was used to merge the individual methcounts files from multiple replicates to produce a single estimate with higher coverage for each species. The ‘hmr’ program was then used to identify HMRs in each species. This program defines the methylation levels at single sites while considering the number of reads affecting those levels using a hidden Markov model method with a β-binomial distribution. We used MethBase [45] from the UCSC Genome Browser track hub to visualize the methylation profiles and HMRs.

Expression analysis

To investigate the association between DNA methylation level in DMP-related genes and early gene expression, we used stage-specific modules and looked at gene expression levels at different time points of human and mouse early embryos. Given that the zygotic genome activation (ZGA) starts during the one-to-two-cell stage in mice [46,47] and the four-to-eight-cell stage in humans [47–49], we obtained single-cell RNA sequencing reads for the oocyte and two-cell stages in mouse and for the oocyte and eight-cell stages in human from the GEO database (accession: GSE44183) [50].

We aligned the reads to the human (hg38) and mouse (mm10) genomes, using STAR v.2.5.3a [51]. FeatureCounts v.1.5.3 [52] was used to conduct read summation onto genes. Differential expression testing between the oocyte and ZGA stages was performed using edgeR v.3.14.0 [53] with the quasi-likelihood test, then normalized values of transcripts per million (TPM; adjusted with gene length) were calculated for those genes that were overexpressed by the ZGA in at least one species. OrthoFinder [54] was used to identify orthologous gene groups using all available protein sequences (www.ensembl.org/) for human and mouse. Average TPM values of the cells or embryos from the same stage for the orthologous DMP-related genes were used for visualization.

Results

DNA methylation patterns at CpG sites orthologous to humans, mice, rats & mini-pigs

To compare the DNA methylation pattern of spermatozoa from human, mouse and rat, we extracted WGBS and RRBS public data from the GEO database. For mini-pigs, we performed RRBS on spermatozoa isolated from mini-pig semen samples (Figure 1A). A summary describing the basic features of the data is shown in Supplementary Table 1. Sequencing reads were mapped to in silico-generated bisulfite-converted genomes corresponding to each species. For WGBS, we measured a DNA methylation signal from both strands in 28.8–44.4 million CpG sites per sample, while for RRBS we measured 2.1–7.3 million CpG sites per sample (Supplementary Table 2).

Next, we mapped nonhuman DNA methylomes to the human reference genome at single-CpG resolution (see Methods section). CpG sites from mice and rats were mapped to the human genome, even when their orthologous positions did not correspond to a CpG site in humans; for example, in the context of CHG, CHH (where ‘H’ represents any base except G). This allowed us to map 12.5–17.6 million CpG sites per sample for WGBS data, and 1.6–4.2 million CpG sites per sample for RRBS data (Supplementary Table 2). We annotated 0.9–1.9 million fully conserved orthologous CpGs in the human and mouse genomes (i.e., CpG in both genomes) and 0.6–2.5 million fully conserved CpGs in the human and rat genomes using the RRBS and WGBS data sets, respectively, as well as 1.4 million fully conserved CpGs in the human and mini-pig genomes using the RRBS data set (Supplementary Table 2). For the subsequent analyses, we opted to only use fully conserved orthologous CpG sites.

The methylation levels of orthologous CpGs were mapped against densities of genes, as shown in Supplementary Figure 1A. Comparing human and mouse genomes using WGBS data, we identified a respective distribution of orthologous CpG sites of 14.8, 37.3, 34.4 and 13.5% within promoters, exons, introns and intergenic regions, respectively, and 16.1, 14.1 and 69.8% within CpG islands, CpG island shores and other regions, respectively. Using the RRBS data, the distribution was 54.3, 21.8, 20.9 and 3.0% within promoters, exons, introns and intergenic regions, respectively, and 83.8, 8.1 and 8.1% within CpG islands, CpG island shores and other regions, respectively (Supplementary Figure 1B). Comparing the human and rat genomes and using WGBS data, the respective distribution within promoter, exon, intron and intergenic regions was 14.4, 42.7, 32.1 and 10.8%, and 20.4, 14.3 and 65.3% for CpG islands, CpG island shores and other regions. Using the RRBS data, the distribution was respectively 51.8, 23.5, 21.1 and 3.6% for promoter, exon, intron and intergenic regions and 78.8, 10.4 and 10.8% for CpG islands, CpG island shores and other regions (Supplementary Figure 1B). Comparing human and mini-pig genomes using the RRBS data, the distribution was respectively 50.0, 22.2, 23.8 and 4.0% within promoters, exons, introns and intergenic regions, and 84.8, 6.4 and 8.8% within CpG islands, CpG island shores and other regions (Supplementary Figure 1B). Thus our analysis showed a trend for a higher CpG density in promoter and CpG islands when using RRBS compared with WGBS data.

Multidimensional scaling analysis of methylation levels for orthologous CpGs in humans, mice, rats and mini-pigs showed a very distinct separation of data (Figure 1B). While we found no significant differences between species in CpG methylation across the entire genomes (Scheirer–Ray–Hare test, p_groups = 0.26 and 0.31 based on WGBS and RRBS data sets, respectively), differences were highly significant for specific regions (Scheirer–Ray–Hare test, p_region = 1.06E-89 and 1.86E-24 based on WGBS and RRBS data sets, respectively). Most genomic regions displayed different DNA methylation levels, according to the Dunn test for pairwise comparisons (Table 1). Based on both WGBS and RRBS data sets and across species, we found that the annotation with the highest methylation levels were the exons, particularly the exons beyond the first exons (annotated ‘Rest of exons’ in Figure 1C & D). On the other side of the spectrum, promoter regions exhibited the lowest methylation levels at orthologous CpGs (Figure 1C & D).

Human and mouse orthologous CpG sites identified in WGBS data were found in 13,134 genes, 4839 promoters, 15,147 exons and 10,769 introns. From these CpG sites, we found sites in common with the RRBS data for 6778 genes, 2207 promoters, 3575 exons and 2337 introns (Fisher’s exact test, all p < 0.05) (Figure 1E). When analyzing human and rat orthologous CpG sites based on the WGBS data, we found orthologous CpG sites in 7844 genes with 2224 promoters, 7052 exons and 4679 introns, of which we found sites in common within the RRBS data for 3214 genes, 773 promoters, 1439 exons and 970 introns (Fisher’s exact test, all p < 0.05) (Figure 1F). Human and mini-pig orthologous CpG sites identified in RRBS data were found in 16,028 genes, 14,513 promoters, 12,872 exons and 10,497 introns.

Using WGBS data for each genomic region of interest, the DNA methylation patterns had a bimodal distribution, with important densities of hypermethylated or hypomethylated regions but a low density of intermediate methylation levels in both human–mouse (Figure 1G) and human–rat (Figure 1H) comparisons. For each annotation, the bimodal distribution was differently balanced between hypermethylation and hypomethylation: while introns and exons showed a comparable bimodal DNA methylation pattern with a higher density of hypermethylation, in promoters, CpG sites were largely unmethylated (Figure 1G & H). For each annotation, the methylation density of orthologous CpG methylation showed significant hypomethylation in humans compared with both mice and rats (all p < 0.001). Methylation patterns of all genes/cytosine sites across the different species comparisons and gene bodies in different CpG island regions confirmed that methylation levels within CpG islands are lower compared with CpG island shores (all p < 0.0001) (Supplementary Figure 2).

We next aimed to study the distribution of DNA methylation across gene bodies. As the density of exons and introns containing orthologous CpGs decreased within a gene as the ordinal positions of the exons or introns increased, we opted to compare the methylation levels at promoter regions versus the first 20 exons/introns (Supplementary Figure 3). We observed that the methylation levels of promoters were lower than those of the first 20 exons (all p < 0.0001) (Figure 2A). We found that exons/introns at different ordinal positions have different methylation levels, with the first exon’s methylation levels being lower than those of the remaining exons in all species (all p < 0.0001) (Figure 2A). Likewise, the first introns had lower levels of methylation than the other introns (all p < 0.0001) (Figure 2B).

Differential methylation identifies pathways related to the function of the central nervous system

We next aimed to identify differential DNA methylation for the genomic regions carrying the orthologous CpG sites characterized above. We focused on regions with methylation differences higher than 20% across all comparisons, because we found that most methylation differences were between 0 and 20% in all species (Supplementary Figure 4). The detailed information regarding the DMPs, DMEs and DMIs is listed in Supplementary Tables 3–7.

Comparing humans and mice using WGBS and RRBS, we identified, respectively, 558 and 195 DMPs, 1769 and 578 DMEs, and 1602 and 266 DMIs. We found that the numbers of genes in common between DMPs, DMEs and DMIs were 51 and 18 using WGBS and RRBS, respectively (Fisher’s exact test; all p < 0.05) (Figure 3A). Methylation differences within promoters, exons and introns were similar using WGBS and RRBS data sets, which supports that these two techniques are analogous to detect DNA methylation differences (Figure 3B). Accordingly, our results revealed a good relationship in DNA methylation differences from the same DMRs (Pearson’s r = 0.87) between WGBS and RRBS (Figure 3C). To get insight into the functional consequence of species-specific methylation at conserved regions, we analyzed the GO:BP networks of the genes with proximity to the DMRs that we identified. We found that DMRs were concentrated within 20 and ten GO:BP groups based on the WGBS and RRBS data sets, respectively (Figure 3D & E & Supplementary Tables 8 & 9). A substantial number of GO:BP terms relating to the central nervous system were shared between the data sets, including ‘chemical synaptic transmission’ and ‘nervous system development’.

Comparing humans and rats, we identified 184 and 163 DMPs, 643 and 419 DMEs, and 398 and 175 DMIs using WGBS and RRBS, respectively. We identified 15 and 20 common genes related to DMPs, DMEs and DMIs (Fisher’s exact test; all p < 0.05) (Figure 3F). Similar to the human–mouse comparison, methylation differences within promoters, exons and introns showed a similar distribution using the WGBS and RRBS data sets (Figure 3G). Here also, methylation differences of DMRs identified in the WGBS and RRBS data sets were highly correlated (Pearson’s r = 0.91) (Figure 3H). The biological process network of all identified DMRs mainly consisted of 12 groups based on both WGBS and RRBS data sets (Figure 3I & J & Supplementary Tables 10 & 11). As in the human–rat comparison, many of the biological processes shared were related to the central nervous system, including ‘neurogenesis’ and ‘modulation of chemical synaptic transmission’. Comparing human and mini-pig RRBS data, we identified 609 DMPs, 2223 DMEs and 1460 DMIs, of which 96 were common (Fisher’s exact test; p < 0.05) (Figure 3K). DMRs were enriched in 40 GO:BP groups, several of which were related to signal transduction and central nervous system development (Figure 3L & Supplementary Table 12), similar to human and rodent comparisons. We identified some specific terms including ‘regulation of hair follicle development’, ‘sensory perception of sound’, ‘eye development’ and ‘ear morphogenesis’ (Figure 3L & Supplementary Table 12).

Biological enrichment of DMP & DME1-related genes in human–mouse comparison

Our analysis of WGBS data from human and mouse sperm identified the promoters of MIR127 (Q-value: 1.6 × 10^-162; methylation difference: - 93%), MBD1 (Q-value: 2.3 × 10^-146; methylation difference: - 92%), LRRC30 (Q-value: 8.6 × 10^-101; methylation difference: - 97%) and AMIGO3 (Q-value: 1.3 × 10^-86; methylation difference: - 89%) as the most differentially methylated between species. The genes AMIGO3 and LRRC30 were common between the most DMPs and first exon (Figure 4A & B). GO enrichment analysis of the promoters harboring lower methylation in humans compared with mice returned terms related to muscle development, RNA processing and the development and function of the nervous system, while the DMPs with lower methylation levels in mice than in humans were largely related to genes controlling muscle function, intracellular signaling pathways and the perception of smell (Figure 4C & Supplementary Table 13). To get further insight into the cell functions possibly affected by distinct epigenetic signatures in different species, we plotted links between the genes nearest to the DMRs common to DMPs and DME1 and their corresponding GO terms. Out of the 171 DMRs common to DMPs and DME1, we noticed an important fraction of DMRs near genes related to signal transduction (35 genes) and linked to significant GO terms like ‘olfactory receptor activity’ (GO: 0004984) and ‘detection of chemical stimulus involved in sensory perception of smell’ (GO: 0050911), while four DMRs were near genes linked to GO terms related to neuronal function (e.g., ‘synapse assembly’ [GO: 0007416]) (Figure 4D). Thus, our analysis supports that regions carrying orthologous CpG sites have distinct epigenetic patterns in human and mouse.

Mapping of HMRs in sperm

Taking advantage of the higher resolution provided by WGBS compared with RRBS, we used a hidden Markov model approach to identify HMRs in each species individually. We identified 89,761, 72,917 and 73,185 HMRs with an average size of 2.4, 1.8 and 1.3 kbp for humans, mice and rats, respectively (Supplementary Table 14). Promoters were overlapping more with the HMRs compared with exons and introns, and the density of HMRs decreased toward the 3′ end of the genes (Supplementary Figure 5 & Supplementary Table 15).

We next aimed to study the cross-species differences in the size of HMRs that overlapped with the identified DMPs in the human–mouse comparison. Our results showed that species have distinct HMR boundaries. HMRs that overlapped with the promoters of genes harboring lower methylation in humans compared with mice covered a larger region in humans than in rodents. Five genes showed methylation gain (loss of HMRs) in rodents (Supplementary Figure 6A). Conversely, the promoters of genes with lower methylation levels in mice than in humans largely corresponded to lost HMRs in humans (i.e., HMR gain in rodents) (Supplementary Figure 6B). Thus, our analysis reveals a general evolutionary extension of the HMR sizes in humans compared with rodents for the genes with lower methylation levels in human promoters than in mouse promoters (Supplementary Figure 6C). HMR gain in rodents occurred more frequently in the genes with lower methylation levels in mice than in humans, and the HMRs were not significantly larger in size for these genes in rodents (Supplementary Figure 6D).

Four examples of methylation levels and identified HMRs across the AMIGO3, TAS2R41, SHANK2 and OR6S1 genes in humans and their orthologs in rodents are shown in Figure 5. AMIGO3, identified as a common gene between the top 20 DMP and DME1 genes (Figure 4A & B) is a gene well characterized for its role in brain development [55]. AMIGO3 overlapped with an HMR in humans (light blue line) but had a hypermethylated pattern in mice and rats (Figure 5A). We identified TAS2R41, a gene that encodes a bitter taste receptor, as linked to G-protein-coupled receptor (GPCR) activity and as harboring both DMP and DME1 (Figure 4D). In contrast to AMIGO3, TAS2R41 harbors a hypomethylated promoter and exon in mice and rats, and a hypermethylated promoter and exon in humans (Figure 5B). SHANK2, which was found in three GO terms – ‘memory’, ‘synapse assembly’ and ‘long-term synaptic potentiation’ (Figure 4C–D) – contains an HMR at the start of the gene in humans, but hypermethylated regions in mice and rats (Figure 5C). OR6S1, which showed a higher methylation level in humans than in mice and was involved in three significant GO terms – ‘G-protein-coupled receptor activity’, ‘olfactory receptor activity’ and ‘detection of chemical stimulus involved in sensory perception of smell’ (Figure 4C & D) – contains an HMR in mice and rats but is mostly hypermethylated in the same regions in humans (Figure 5D).

Gene expression in the early embryo is moderately associated with DNA methylation in sperm

Given that methylation of promoters controls gene expression programming [56], we investigated whether the DMPs that we identified in spermatozoa are linked to the expression levels of their nearest genes in early embryos. We studied the expression of genes near human–mouse DMPs in the oocyte and eight-cell stages in humans and in the oocyte and two-cell stages in mice. We used the OrthoFinder program to identify orthologous gene groups in the human and mouse genomes. OrthoFinder returned 16,552 orthologs with one or more genes in mouse and human, of which 14,045 were one-to-one (∼85%). OrthoFinder found only 355 and 270 unique genes in human and mouse, respectively (Supplementary Table 16). Out of 30 genes upregulated by ZGA in humans at the eight-cell stage and with lower methylation levels in humans compared with mice, 19 were similarly upregulated in mice by ZGA at the two-cell stage. Of these genes, eight did not show a significant pattern in mouse by ZGA, and three were significantly downregulated between the oocyte and two-cell stages in mice, showing an expression pattern loaded from the maternal side (Figure 6A & Table 2). Out of 15 genes with lower methylation levels in mice compared with humans, 12 were upregulated in mice by ZGA but did not show a significant pattern in humans, and three were similarly upregulated in humans by ZGA (Figure 6B & Table 2). Promoters of genes with lower methylation levels in mice compared with humans consistently exhibited higher transcription abundance (Figure 6C), thereby validating our approach.

Discussion

Here we performed an epigenomic comparison of DNA methylation using WGBS and RRBS data from spermatozoa of humans, mice and rats and RRBS data from spermatozoa of mini-pigs. We identified that, while a substantial fraction of epigenome methylation marks are conserved across evolution, some orthologous CpG sites are differentially methylated. Strikingly, the pathways of genes differentially methylated between species were related to distinct functions. Our analysis supports an evolutionary selection of species-specific DNA methylation footprint in orthologous genomic sequences that allows speciation of organ function.

Our study confirmed a bimodal distribution in sperm-cell DNA methylation, consistent with previous reports [23,25]. This bimodal distribution pattern may be responsible for maintaining the basal transcription profile of the preimplantation embryo [57]. Moreover, our study showed that the first exons and first introns had lower levels of methylation than the other exons and introns. This is in accordance with the results of previous findings in pig [58] and fish [59] that reported a lower methylation level on the first exons and introns. Expression levels of the exon gene product are negatively associated with the DNA methylation levels of the first exons [60] and the first introns [59], thus supporting a functional role of sperm DNA methylation on gene expression post-fertilization in the early phase of embryogenesis.

GO terms of orthologous genomic regions found to be differentially methylated when comparing humans with mice, rats and mini-pigs were largely found near genes associated with the development of the central nervous system and signal transduction. Thus genomic regions under epigenetic evolution may be involved in species-specific phenotypic variation, with humans using this means for adapting cognitive function.

We found that promoters with lower methylation in humans compared with rodents are near genes related to the development and function of the nervous system. Conversely, promoters with hypomethylation in rodents compared with humans appeared to be enriched in the proximity of genes related to intracellular signaling and the perception of smell. These results are in line with a study showing that genes with human-specific hypomethylated promoters are related to the development of the central nervous system and enriched for genome-wide association study signals of brain-related traits [23]. In that study, genes with cattle-specific hypomethylated promoters were found to be primarily involved in lipid storage and metabolism [23]. Similarly, a study comparing human and chimpanzee sperm methylomes identified that genes linked to human-specific HMRs are linked to neuronal functions [61]. Previously, we discovered that the differences in DNA methylation in spermatozoa from obese individuals after gastric bypass-induced weight loss are located near genes linked to brain function [62]. Genes referenced as participating in brain function and development appear to carry hot spots of epigenetic variation, as endurance exercise training also altered the DNA methylation signature of human spermatozoa at gene hot spots related to brain function and development [20]. Given that DNA methylation at promoter and CpG-rich regions plays a role in gene repression, our results suggest that promoter hypomethylation in sperm is associated with transcriptional activation of genes after fertilization. A function of DNA methylation in spermatozoa on developmental programming and brain function is further supported by the observation that age-related DNA methylation changes in sperm are passed on to the next generation and contribute to the children of older fathers having an elevated incidence of neuropsychiatric illnesses [63,64]. Thus, we speculate that the hypomethylation of genes controlling the central nervous system in humans contributes to the speciation of brain function in humans compared with rodents.

We identified that genes hypomethylated in rodents are enriched for pathways related to the perception of smell. It was previously posited that genes with promoters that lack HMRs in both human sperm and human embryonic stem cells have a strong enrichment for GPCR and sensory receptors that detect and respond to smell [61]. Odor molecules are detected by olfactory receptors (ORs), which constitute the largest group of GPCR genes. ORs include approximately 350 intact members in humans and approximately 1000 intact members in mice and rats [65]. These genes are largely expressed in sensory neurons found in the main olfactory epithelia of the nasal cavity. The difference in the number of genes between these species appears to be the consequence of both a significant increase in OR genes in the rodent lineage and a large deactivation of OR genes in the human lineage following the primate–rodent divergence [66]. It has been suggested that inactivation of OR genes in higher primates may be the consequence of the development of a full trichromatic vision which aids higher primates in finding food, mates and territory, potentially explaining a decreased need for olfaction [67]. Here we identified TAS2R41, a gene that encodes a bitter taste receptor and which is linked to GPCR activity, which harbors a hypomethylated promoter and exon in mice and rats and a hypermethylated promoter in humans. These findings support that the TAS2R41 gene is selectively silenced in humans, which is in line with previous studies showing a downregulation of taste receptor genes in humans compared with rodents [68,69].

Our results identified a global trend for size extension in HMRs of humans which overlapped with promoters of the genes with lower methylation levels in humans compared with mice. We found rodent-specific HMR gains for genes with lower methylation levels in mice compared with humans. These observations are consistent with a previous study showing that rodent genomes exhibit a stronger tendency toward the birth of novel HMRs, whereas HMR boundaries of promoters have evolved to be significantly longer in primates than in rodents [25].

It is fascinating to hypothesize that HMR size extension in human promoters might be associated with the evolution of new regulatory elements to enhance specificity of the transcriptional regulation of the highly specialized genes involved in the development of the nervous system. A particularly intricate regulation of gene expression is required for the highly specialized functions of neurons [70], which may be accomplished, at least partially, by utilizing the adaptability of transcription factor binding sites. A previous theoretical study on evolution of these sites showed that gain of transcription factor binding sites can be facilitated by the availability of longer regulatory sequences in which the sites can evolve simultaneously and by the biophysical cooperation between transcription factors, which can also match the theoretical time-scale computations derived from comparative genomics [71]. Lineage-specific promoter HMR extensions revealed enhanced RNA polymerase II and EP300 binding sites, and enrichment of numerous transcriptional activators [25].

Our findings that humans and rodents have different types of genes that are hypomethylated supports a mechanism by which orthologous genes are selectively used between species to allow organ specialization. Whether species-specific epigenetic reprogramming of orthologous genes has occurred specifically in response to environmental factors, or may have been selected from stochastic epigenetic variation, remains to be established.

Previous studies suggest that the epigenomic status of embryonic cells and the early activation of developmental genes are influenced by the methylation of germ-cell DNA [72]. Based on our comparison of the gene expression profiles during early embryo development for genes near DMPs, we identified that there are slightly more gene expression differences between the oocyte and eight-cell stages in human than between the oocyte and two-cell stages in mouse when only considering genes that display lower promoter methylation in human compared with mouse. Additionally, as expected, we confirmed that promoters of genes with a lower methylation level in mouse than in human consistently exhibited more gene expression differences between the oocyte and two-cell stages in mouse than between the oocyte and eight-cell stages in human. Together, interspecies differences in gene expression levels during early embryogenesis are modestly associated with promoter DNA methylation profiles in spermatozoa. It is worth mentioning that because the RNA expression and DNA methylation data sets were from independent studies, and therefore distinct tissues or cells, we were unable to make a paired analysis where we could more directly link changes in methylation to changes in gene expression. On the other hand, by using independent samples, our analysis was less affected by internal experimental biases.

Here we found that, when overlapping CpGs annotated by both WGBS and RRBS within gene features, the proportion of cytosines annotated to non-promoter regions decreased, whereas the percentage of cytosines annotated to the promoter section of the genome increased. Additionally, the proportion of cytosines annotated among CpG island shores decreased, whereas the percentage of cytosines annotated within CpG islands increased. This discrepancy in mapping using the two sequencing techniques may be attributed to low coverage of CpGs within CpG-poor regions by the RRBS technique due to the restriction enzyme cut sites used in RRBS. Indeed, RRBS mainly targets promoters and CpG islands because the regions sequenced are flanked by the CCGG sequence, a CG-rich site [73]. It should be noted that the RRBS method superficially covers CpG shores and other intergenic regions and that RRBS data sets have a lower average methylation level compared with that of WGBS data sets, as large stretches of repeated elements in noncoding DNA regions are mainly methylated [74]. Given that the WGBS approach can generate a lot of reads in poorly assembled noncoding DNA regions, this could explain why WGBS data have a lower mapping yield than RRBS data [75]. For instance, even though all CpG sites should, in theory, be covered, some of them have low coverage (one to ten counts) or are not sequenced using the WGBS method [73]. Therefore the average read depth of the RRBS method was greater compared with the WGBS method in our study. Our results highlight the relevance of RRBS as a sequencing method for measuring CpG methylation in comparative studies among species, when sequencing costs, read coverage and the availability of enough methylation information are taken into account.

Conclusion

In conclusion, our detailed comparative analysis identified a divergence in DNA methylation patterns in spermatozoa from rodents, mini-pigs and humans, although the underlying DNA sequence is highly conserved. Our results revealed that the differences in the methylation footprint between human, rodents and mini-pigs are associated with species-specific phenotypic traits. Divergent epigenetic marks on orthologous sequences suggest a mechanism driving phenotypic plasticity and organ speciation, allowing long-term and stable adaptation to the environment.

Future perspective

The growing availability of epigenetic editing techniques may provide a way to test the function of hypomethylated orthologous sequences in sperm on the phenotype of the offspring. Future studies comparing the sperm DNA methylome of additional species will be useful to identify the universal factors driving species adaptation by epigenetic means.