Molecular Basis of Sperm Methylome Response to Aging and Stress

TL;DR

The paper explores how aging and stress affect sperm DNA methylation in regions linked to embryonic development, proposing that these changes may help diversify offspring traits to improve survival.

Contribution

The study introduces a novel hypothesis that stochastic methylation variation in variable methylation regions (VMRs) explains sperm methylome responses to stress and aging.

Findings

Changes in sperm DNA methylation due to aging and stress occur predominantly in VMRs associated with developmental genes.

VMRs are enriched for a binding motif for ZFP42, an epigenetic remodeler involved in genomic imprinting.

The nature of DNA regions, not the stressor type, determines the sperm methylome's response to aging and stress.

Abstract

Many human and animal studies have shown that epigenetic programs in sperm are affected by aging and various stressful factors, and that these changes may have adverse effects on offspring health. Surprisingly, different stressors often affect similar sets of differentially methylated regions (DMRs) associated with developmental genes regardless of the nature of the stressor. We hypothesized that these unexplained patterns of sperm epigenome response may originate from stochastic epigenetic variation—a hypothesis suggesting that, in response to stressors, naturally variable methylation regions (VMRs) associated with morphogenic genes exhibit increased methylation variation to diversify offspring phenotypes and improve the chances of survival of the genetic lineage. We used data on sperm DNA methylation from mouse and rat experiments and demonstrated that, indeed, changes in DNA methylation in response to aging and various stressors predominantly occur in VMRs associated with major developmental pathways. Next, we propose a model that explains how a stochastic increase in methylation variation in VMRs results in either increased or decreased methylation in these regions. Further, we show that VMRs are enriched for a binding motif for ZFP42, an epigenetic remodeler involved in genomic imprinting. This knowledge may open opportunities to develop interventions to control epigenetic information transfer via germ cells.

Aging and stress-related factors affect sperm DNA methylation in regions associated with genes responsible for embryonic development. The stochastic epigenetic variation hypothesis holds potential to explain these patterns, proposing that, in response to stressors, naturally variable methylation regions (VMRs) associated with morphogenetic genes exhibit increased methylation variation to diversify phenotypes and improve the chances of survival of the genetic lineage. Here, we test predictions from this hypothesis using mouse and rat sperm DNA methylation data from publicly available sources. Specifically, we identify VMRs and analyze their overlap with regions differentially methylated (DMRs) in response to aging, stressors, and with various genomic elements. We demonstrate that the nature of the DNA regions, rather than the nature of the stressor, determines the response of the sperm methylome to aging and stress, and propose a model that explains shifts in methylation within VMRs through stochastic changes, whereby initially hypermethylated regions lose methylation and initially hypomethylated regions gain methylation. VMRs are depleted of open chromatin regions and histones in male germ cells and are enriched for a binding motif for ZFP42, an epigenetic remodeler. This knowledge may open opportunities for the development of interventions to control epigenetic information transfer via germ cells.

Figures

Full text

1. Introduction

Research conducted over the last two decades demonstrates that sperm DNA methylation is sensitive to reprogramming by a broad range of factors including paternal diet [1], childhood abuse [2], psychological stress [3], metabolic status [4,5], heat stress [6], environmental exposures during perinatal development or the preconception period [7,8,9,10,11], and aging [7,8,12,13] among others. After fertilization, parental-specific epigenetic marks of gametes undergo reprogramming to establish totipotency in the developing embryo. However, imprinted loci and certain classes of repetitive sequences escape this reprogramming event and contribute to the non-Mendelian form of inheritance [14]. Sperm DNA methylation in other genomic loci may also resist reprogramming in humans [15,16] and laboratory rodents [3,13,17,18], suggesting that the final sperm methylome is a significant channel for the transfer of inheritable information to the next generation [19].

Final DNA methylation profiles of spermatozoa are a result of epigenetic events that occur during spermatogenesis [20]. Specifically, reduction in global methylation of around 12–13% occurs in preleptotene spermatocytes and methylation is gradually reestablished during leptotene-pachytene stages [20,21,22]. Additionally, transition from elongating spermatids to late spermatids is associated with changes in methylation of more than 5000 DNA regions [23].

Published research demonstrates that changes in sperm epigenome mediate the effects of paternal factors on adverse outcomes in offspring. Some examples include the effects of advanced paternal age on life span reduction [17], reduced exploratory and startle behavior [13] in mouse offspring, and reduced fertilization rates in humans [24]; metabolic dysfunction in mice induced by paternal psychological stress [3], paternal diet [1], and paternal prediabetes [4], and increased risk-taking behavior induced by paternal psychological stress in mouse offspring [3]. The role of sperm epigenetics in offspring health is also suggested by population studies linking advanced paternal age with poor pregnancy outcomes [25], lower odds of live birth [26], and adverse health of the offspring at later life [27,28,29], including increased susceptibility to the early development of cancer [30]; and neurodevelopmental disorders such as schizophrenia [31] and autism [28,29,32,33]. These changes in the sperm epigenome induced by age or various stressors are often explained by the accumulation of random epigenetic errors or epimutations. However, as more data on sperm epigenetics are becoming available, common patterns of epigenetic change are emerging and challenge the random nature of sperm epigenetic changes [34].

First, two epigenetic mechanisms, DNA methylation and small non-coding RNA (sncRNA), undergo concordant changes in sperm. For example, Xie et al. [17] demonstrated that in mice both age-dependent changes in DNA methylation and sncRNA in sperm were associated with major developmental and aging pathways. In our experiments with rats, age-dependent changes in sncRNA [8] and DNA methylation [7] were enriched in the same set of developmental pathways, including those identified by Xie et al. In another study, increased adiposity and reduced glucose tolerance in mouse offspring were equally induced by sperm and seminal plasma of fathers on a low protein diet [1], suggesting that sncRNA cargo of extracellular vesicles in seminal plasma [35,36,37] and spermatozoa-specific epigenetic mechanisms transfer similar epigenetic information to the next generation.

Second, in a broad range of studies of the effects of various factors on sperm DNA methylation regardless of the factor analyzed (e.g., age, diet, chemical exposure, heat stress, psychological stress) and biological species (human, rat, mouse) genes associated with differentially methylated regions (DMRs) in sperm are almost always enriched for the same set of biological categories, including embryonic development, neurodevelopment, behavior, and metabolism [2,4,5,6,7,9,12,16,17,18,38,39,40,41]. These conserved patterns of sperm epigenetic changes suggest that a purely stochastic error model cannot explain the shared changes in sperm epigenome and strongly suggest an evolutionary basis of this observed phenomenon [42].

According to the evolutionary stochastic epigenetic variation (SEV) hypothesis, stress may increase variation in naturally variable methylation regions (VMRs) associated with development and morphogenesis to ensure increased phenotypic variability in offspring and improve chances of survival of the corresponding genetic lineage in conditions of a changing environment [43]. This hypothesis predicts that sperm DNA contains VMRs associated with genes involved in embryonic development and that VMRs non-specifically respond to stressors by an increase in methylation variation. This hypothesis also predicts that differential methylation in response to aging or stressors may result from increased variation in VMRs, although the mechanism of conversion of VMRs to DMRs remains unclear.

To our knowledge, these predictions as well as the SEV hypothesis in general were never tested experimentally, and it is not clear if sperm epigenetic changes observed in many published studies may be explained by such a mechanism. Thus, the goal of the current study consists of (1) the testing of predictions from the SEV hypothesis using sperm DNA methylation data and (2) the development of evolutionary and molecular understanding of the nature of DNA methylation changes induced by aging and stressors. This study demonstrates concordance of the sperm DNA methylation responses to aging and stressors to the SEV hypothesis and explains some molecular events behind these responses.

2. Materials and Methods

2.1. Datasets

This study is based on sperm DNA methylation datasets previously published by our and other groups (Table S1, Supplementary Materials). To achieve uniformity of data analysis we used only studies using reduced representation bisulfite sequencing (RRBS) as the major method for DNA methylation analysis and reanalyzed all data using a pipeline described below. The following datasets were used to analyze changes in sperm DNA methylation in response to aging and other stressful factors. To identify VMRs in young, untreated wildtype mice we used RRBS data from 8 individual 56-day-old animals from our recent study [41,44]. DNA methylation profiles of 56- and 154-day-old mice from the same study [41,44] were used to identify age-dependent changes in mouse sperm epigenome. The suitability of this dataset for the identification of changes in sperm methylome associated with aging was provided in the original study [41], which demonstrated that the absolute majority of DNA methylation regions undergo a linear increase or decrease in DNA methylation during the studied period from postnatal day 56 to 334. Similarly, a linear association between sperm DNA methylation and age across a broad age range was shown in humans [24,45]. Therefore, any two time points that are sufficiently spaced to detect methylation differences can be used to identify age-related changes. Additionally, the study by Amir et al. was the source of RRBS datasets representing wildtype (WT) mice and mice with Sertoli-specific knockout of Rictor—a genetic manipulation resulting in a disruption of the blood-testis barrier (BTB) [41,44]. To analyze age-dependent changes in sperm methylome in rats we used our RRBS data from 65- and 120-day-old animals from previously published research [7,46]. The same dataset was used to analyze the effects of perinatal exposure to 0.2 mg/kg/day of a brominated flame retardant 2,2’,4,4’-tetrabromodiphenyl ether (BDE-47) which was conducted via pipette feeding of pregnant/lactating dams to a solution of BDE-47 in corn oil or vehicle starting from gestational day 8 till post-partum day 21. Changes in mouse sperm methylome in response to phthalate exposure were analyzed using data from our previous study [18]. In this study, sperm was obtained from 17- to 18-week-old animals following continuous exposure to 25 mg/kg/day di(2-ethylhexyl) phthalate (DEHP), a compound found in plastics and personal care products, or vehicle control (50% dimethyl sulfoxide and 50% ethanol) via subcutaneous osmotic pumps for the duration of 2 spermatogenesis cycles. Changes in mouse sperm methylome in response to cadmium exposure were analyzed using data from a published study [9] in which animals were exposed to approximately 0.69 mg of CdCl_2_ per kg of body weight for 9 weeks via drinking water. Sperm DNA methylation data used in the current study [7,18,41] were produced following an optimized DNA extraction protocol that ensured DNA yield >90% [47]. The bisulfite conversion rate for these datasets was >99%, calculated as the number of genomic cytosines outside a CpG context that are unconverted, divided by the total number of cytosines outside a CpG context. This conversion rate is comparable with conversion rates for DNA extracted from fresh or frozen somatic tissues [48].

Data on mouse sperm ATAC-seq and CHIP-seq with antibodies for various histones were taken from ChIP-Atlas a data-mining suite powered by full integration of public sequencing data. RepeatMasker was used to identify DNA sequences for interspersed repeats and low-complexity DNA sequences. The lists of imprinted genes were taken from the Geneimprint database [49]. Finally, the list of genomic regions that survive epigenetic reprogramming in mice was taken from the study which provides a comparison of DNA methylation values in sperm, oocytes, blastocyst, and embryonic stem cells [50].

2.2. Identification of DMRs

Raw RRBS reads were trimmed using TrimGalore (v 0.6.6) and a NuGEN-specific adaptor trimming script (trimRRBSdiversityAdaptCustomers.py) available from GitHub (nugentechnologies/NuMetRRBS). Trimmed reads were mapped using Bismark-Bowtie2 with no mismatch allowed. GRCm38/mm10 and mRatBN7.2/rn6 genome assemblies were used to map mouse and rat data respectively. Methylation counts were called using Bismark extract. Differentially methylated regions were identified using the Methyl kit (v 1.24.0) pipeline. In brief, the genome was tiled into sliding windows and a weighted methylation level was calculated for each window. To minimize error in base calling, we filtered out bases with less than 10× coverage and read counts more than the 99.9th percentile, and coverage values were normalized using default settings. We used a logistic regression model for p-value calculation subsequently adjusted for multiple comparisons (FDR) using the SLIM method for final DMR identification. Individual DMRs were identified for a 100 bp sliding window with q < 0.05.

2.3. Identification of VMRs

Methylation variance was calculated for each 100 bp region for mouse (56 and 154 days of age, n = 8 per age) and rat (65 and 120 days of age, n = 6 per age) sperm DNA (Table S1). Top VMRs were then identified for each age/species as regions with the highest values of interindividual variance, using a method for the identification of cut-off points in descriptive high-throughput omics studies [51]. This approach assumes that the variance of DNA regions follows a generally biphasic distribution with the majority of regions having low variance and a smaller number of regions having high variance. The method identified the threshold between low-variance and high-variance as a maximum inflection point in a ranked distribution of all regions.

2.4. Enrichment Analyses

Fisher’s exact test was used to test the significance of overlaps of different genomic regions (VMRs, DMRs, ATAC-seq, ChIP-seq, repetitive elements, regions escaping methylation reprogramming, and genes associated with VMRs in different species). RRBS method used for DNA methylation analysis may introduce some biases as it selects only for regions enriched with CCGG sequence. To avoid potential biases introduced by RRBS, we conducted all enrichment analyses only for regions that were identified by RRBS. To identify fold enrichment/depletion of VMRs with other genomic regions (VMRs, DMRs, ATAC-seq, ChIP-seq, repetitive elements, regions escaping methylation reprogramming, and genes associated with VMRs in different species) we compared the observed number of overlapping elements with the calculated number of overlapping elements based on the probability of overlap if no interaction between the two elements exists (null hypothesis).

Additionally, gene set enrichment analysis (GSEA v 4.4.0) was used to visualize overlaps between DMRs and VMRs and establish if DMRs are enriched with VMRs. For functional enrichment analysis DMRs or VMRs were assigned to the closest gene (<5 kb upstream transcription start site, promoter, 5′UTR, exon, or intron) or intergenic region. Spatial genomic annotation was conducted using annotatr package (v 1.24.0) and annotated to the genomic features of the Ensembl genome. Genomic features were compared using the GenomicRanges package (v 1.50.2). Metascape (v 3.5) analysis was used to analyze biological categories associated with genes annotated to DMR or VMRs.

2.5. Motif Analysis

The lists of VMRs and non-VMRs (control) were submitted to STREME tool of the MEME Suite 5.5.5 to identify DNA motifs significantly overrepresented in VMRs. These motifs were submitted to the Tomtom tool of the MEME Suite to identify matching annotated motifs in the HOCOMOCO Mouse (v11 CORE) database. Given that DNA motifs are not yet well annotated in the rat genome, this analysis was done for mice only.

3. Results

3.1. VMRs Are Stable Across Ages and Conserved in Rats and Mice

To identify VMRs in sperm of rodent models, we leveraged our previously published RRBS datasets [7,41] (Table S1). We identified 48,948 overlapping DNA methylation regions in 56- and 154-day-old mice. In this overlapping list, 2566 and 3248 regions were VMRs in 56- and 154-day-old mice, respectively, with 1345 overlapping VMRs between age groups (Figure 1A). In the rat dataset, 279,439 DNA methylation regions were identified in both 65- and 120-day-old animals. In this list, 9872 and 17,418 regions were VMRs in younger and older rats, respectively, with 3589 overlapping VMRs between age groups (Figure 1B). In both species, VMRs identified in younger and older groups overlapped highly significantly with Fisher’s exact test p < 0.00001. These overlaps were 8-fold and 6-fold higher in mice and rats, respectively, than overlaps predicted based on an assumption of independence of methylation variance in young and old animals. In both species, the amount of VMRs increased in older animals, suggesting age-induced gain in regions with high methylation variation. Notably, genes associated with VMRs overlapped between mice and rats (p < 0.00001), indicating their conservation across species. Overall, these results suggest that, concordant with the stochastic epigenetic variation hypothesis, DNA regions prone to methylation variation in sperm are not random and aging is a significant factor increasing the number of such regions.

3.2. DMRs Induced by Age and Stressors Are Enriched for Developmental Pathways

Similar sets of categories were enriched in DMR-associated genes in sperm in response to a broad range of factors (age, diet, chemical exposure, heat stress, psychological stress) in different species (human, rat, mouse) [2,4,5,6,7,9,12,16,17,18,38,39,40,41], suggesting the common mechanisms of DNA methylation response to various stressors. We analyzed enriched categories of genes associated with DMRs induced by age in mice and rats [7,41]; environmental pollutants 2,2’,4,4’-tetrabromodiphenyl ether (BDE-47) (rats, perinatal exposure) [7], phthalates (mice, preconception exposure) [18], and cadmium chloride (mice, preconception exposure) [9]; and blood-testis barrier (BTB) disruption by genetic manipulations in mice [41] (see Table S1 for the details of datasets used in this analysis). In all cases, DMR-associated genes were highly enriched with developmental categories, including, for example, cell morphogenesis, neuron projection morphogenesis, synaptic signaling, axon development, and behavior (Figure 1C). To test if RRBS approach selects for DNA methylation regions associated with developmental and neurodevelopmental genes, we conducted several enrichment analyses with 500 or 2000 random DNA methylation regions identified by RRBS in mouse sperm. Random regions did not demonstrate enrichment for developmental and neurodevelopmental categories observed for DMRs responding to stressors and aging—see example in Figure 1C. These later patterns suggest a possibility that the context of DNA regions determines the specific response of the sperm DNA methylation rather than the nature or timing of factors causing DNA methylation changes.

3.3. VMRs Are Enriched for a Similar Set of Developmental Pathways

We used the lists of genes associated with significant VMRs in young mice and rats for enrichment analysis using Metascape (v 3.5). In both species, VMR-associated genes were highly enriched with developmental categories, including, for example, cell morphogenesis, neuron projection morphogenesis, and axon development (Figure 1C). Similar categories enriched in DMR- and VMR-associated genes in sperm suggest a functional link between DMRs and VMRs.

3.4. Age-Dependent DMRs Overlap with VMRs

To test a hypothesis that age-dependent DMRs are a result of increased variation in VMRs, we analyzed the overlap between VMRs and DMRs in mouse and rat sperm using GSEA. First, age-induced methylation changes were identified across all DNA methylation regions. Next, VMRs in young animals were plotted against DNA methylation regions ranked by their percentage of age-dependent methylation change. The results of these analyses are shown in Figure 2A,B. Indeed, all significant VMRs are grouped in the leftmost and rightmost regions of the differential methylation distribution, indicating that all DNA regions that undergo hypo or hypermethylation with age correspond to VMRs in both species, while DNA regions not changing in methylation with age were depleted with VMRs.

3.5. Stress-Induced DMRs Overlap with VMRs

We further hypothesized that an increase in methylation variation in VMRs may also be a cause of differential methylation induced by stressors, such as chemical exposures and genetic perturbation affecting testis physiology. To test this hypothesis, we used the same GSEA approach as described in the previous paragraph. VMRs in young mice or rats were plotted against DNA methylation regions ranked by their percent methylation change induced by perinatal exposure to BDE-47 in rats, preconception exposure to phthalates or cadmium in mice, and BTB disruption via Sertoli-specific KO of Rictor in mice. Despite the very different nature and timing of these interventions, changes in DNA methylation occurred in VMRs for all analyzed datasets (Figure 2C–F). These results indicate that the genomic context of these regions determines non-specific response to various stressors, regardless of the type of stressor.

3.6. DMRs Are a Result of the Increased Variability of VMRs

The observed link between VMRs and DMRs is unclear at first glance as VMRs respond to stressors by an increase in variation, while DMRs represent regions with increased or decreased methylation values, a seemingly different phenomenon. To understand the functional link between VMRs and DMRs, we developed the following hypothesis. Given that DNA methylation has a binary nature (presence or absence of methyl group on a CpG site), the majority of methylation values averaged for a DNA region, a population of cells, and/or organisms have values close to extreme methylation—namely approximately 0 or 100%. An increase in methylation variance in such regions will result in an asymmetric shift in averaged methylation values (Figure 2G) as more cells will gain methylation in an initially hypomethylated site and more cells will lose methylation in an initially hypermethylated site. This model explains how an increase in the stochasticity of VMRs may result in DMRs. Indeed, the distribution of methylation values in age-dependent DMRs in young and mature mice and rats demonstrate that predominant age-dependent changes consist of an increase in methylation of sites hypomethylated in young animals and a decrease in methylation of sites hypermethylated in young animals (Figure 2H,I).

3.7. Overlap with Repetitive Elements Does Not Characterize VMRs

Given that one of the most important functions of DNA methylation consists of suppression of repetitive elements, we analyzed whether VMRs differ from other methylation regions in their association with these elements. We analyzed the overlap of VMRs and all methylation regions in DNA methylation datasets in young WT untreated mice (Table S2) and rats (Table S3). Most repetitive elements were present in less than 1% of all DNA methylation regions and/or VMRs in both species. Thus, association with these elements cannot explain higher methylation variance in VMRs in comparison with other regions. Two types of repetitive elements were present in higher than 1% of regions in mice: low complexity DNA was found in 1.15% of all regions and 0.31% of VMRs, while simple repeats overlapped with 6.14% of all regions and 2.3% of VMRs. Thus, VMRs were 3.4-fold depleted with low complexity DNA and 2.7-fold for simple repeats with Fisher’s exact test p < 0.00001 for both element types. In rats, four types of repetitive elements were present in higher than 1% regions in all regions and/or VMRs: LINE were found in 1.21% of all regions and in 2.86% of VMRs, LTR were found in 2.56% of all regions and in 4.51% of VMRs, SINE were found in 2.96% of all regions and in 3.14% of VMRs, and simple repeats overlapped with 2.80% of all regions and 2.59% of VMRs. Out of these elements, VMRs were 2.35-fold enriched for LINE and 1.76-fold for LTR with Fisher’s exact test p < 0.00001 for both repetitive elements. The low percent of co-occurrence between methylation regions and repetitive elements in both VMRs and non-VMRs suggests that differences in methylation variance are likely not linked causatively to the presence or absence of repetitive elements.

3.8. VMRs Are Depleted in Open Chromatin Regions

We then hypothesized that higher variance in methylation observed in VMRs may be due to higher accessibility of corresponding regions for enzymes involved in methylation and demethylation of DNA. Given the DNA methylation profile of spermatozoa results from methylation/demethylation events that occur earlier in spermatogenesis, we tested whether VMRs overlap with open DNA regions in spermatogonia, spermatocytes, spermatids, and mature sperm. For this analysis, ATAC-seq data for each cell type were extracted from ChIP-Atlas [52,53,54] and compared with 48,948 methylation regions in sperm, and enrichment/depletion of open chromatin regions were calculated for 2566 VMRs as compared with 46,382 non-VMRs (Table 1). Unexpectedly, we observed that sperm VMRs are significantly (Fisher’s exact test p < 0.00001) depleted with open chromatin regions in all male germ cells. Interestingly, the number of VMRs associated with open chromatin is decreasing in the course of spermatogenesis. ATAC-seq is known to capture only hyper-accessible regions [55]; thus, depletion of VMRs in ATAC-seq regions suggests that their differential features may be due to binding with some biological molecules.

3.9. VMRs Are Differentially Enriched with Histones

Given that VMRs are depleted of open chromatin regions, we hypothesized next that VMRs may be associated with histones or specific modifications of histones. To test this assumption, we extracted ChIP-seq data for spermatogonia, spermatocytes, and spermatids from ChIP-Atlas [52,53,54]. These datasets included data on the following histones: H3K27ac, H3K27me3, H3K4me1, H3K4me2, H3K4me3, and H3K9me3. Enrichment/depletion of histone regions was calculated for 2566 VMRs as compared with 57,988 non-VMRs for germ-cells (Table 2). Given that around 99% of histones are replaced by protamines in mature spermatozoa, spermatozoa were excluded from this overlap analysis. Overall, these results demonstrate that VMRs are significantly differentially enriched with most histone modifications in germ cells and are generally depleted of histones throughout spermatogenesis.

Given that the most significant DNA methylation changes in the course of spermatogenesis occur at spermatocyte [20,21,22] and spermatid stages [23], association of VMRs with histones during these stages may be more important. In spermatocytes, VMRs are significantly depleted with H3K27ac, H3K27me3, H3K4me1, H3k4me2, and H3K4me3. Spermatid VMRs were depleted of H3K27ac, H3K4me2, and H3k4me3. It is reasonable to assume that depletion of every specific histone modification in VMRs in comparison with non-VMRs results from the overall smaller number of nucleosomes associated with VMRs. From that point of view any positive enrichment with histone marks deserves special attention as going against the overall trend for histone depletion. VMRs were only enriched with H3K27me3 and H3K4me1 in spermatogonia and with H3K9me3 in all germ cells.

H3K4me1 is a well-recognized mark of poised or active enhancers, while H3K27me3 and H3K9me3 are major repressive histone modifications usually attributed to facultative and constitutive heterochromatin respectively [56,57,58,59].

These results demonstrate that VMRs are associated with distinctly different profiles of modified histones than non-VMRs. However, given only a small fraction of VMRs is associated with differentially enriched/depleted histone variants, it is unlikely that these associations causally determine the variable nature of VMRs.

3.10. Age-Dependent DMRs Are Differentially Enriched with Histones

Methylation and demethylation of DNA are conducted by different molecular mechanisms which may have different links with the chromatin environment. For example, it was shown that abnormal de novo methylation takes place preferentially in CpG islands associated with H3K27me3 mark and polycom repressive complex 2 [60], suggesting that DMRs undergoing age-dependent increase in methylation may be enriched with H3K27me3 mark as compared with DMRs undergoing age-dependent decrease in methylation. To test if age-dependent hypo and hypermethylated age-dependent DMRs have different patterns of association with histone regions we conducted an overlap analysis for DMRs and histones in the same way as described in the previous section for VMRs.

Similarly to VMRs, age-dependent DMRs are depleted with histones across spermatogenesis (Table S4) and both VMRs and DMRs are enriched only with H3K27me3 and H3K4me1 in spermatogonia and H3K9me3 in all germ cells.

Within DMRs, regions hypomethylated with age are depleted with all types of histone modifications as compared with hypermethylated regions. This finding is concordant with the known dependence of de novo methylation on histones and the known independence of DNA demethylation machinery on histones. Indeed, it was demonstrated previously that the N-terminus of histone H3 tail is required for de novo methylation by Dnmt3a and Dnmt3b [61], while both passive and TET-mediated demethylation do not require histones [62].

3.11. VMRs Distribution Across Genomic Elements and Motifs

To analyze if VMRs are differentially associated with different genomic features, mouse VMRs and non-VMRs were annotated to the genomic features of the mouse genome. The results of this analysis are shown in Figure 2J. Non-VMR DNA methylation regions captured by RRBS are enriched in promoters, while the highest number of VMR overlap is with intergenic regions.

To test if VMRs are characterized by specific binding sites, we submitted the lists of VMR and non-VMR sequences in the sperm of young untreated mice to the STREME tool [63]. This analysis identified 14 differentially enriched motifs with Fisher’s exact test p < 0.05 (Table S5). Only five of these motifs had an E-value (Fisher’s exact test p-value adjusted for multiple comparisons) < 0.05. Given that the goal of this analysis is to identify discriminative features of VMRs, we focused further on the only significant (E-value < 0.05) motif (AARATGGMWW, Figure 2K) that was present in the majority (90.5%) of VMRs. Subsequently, we submitted the selected motif to the Tomtom tool of the MEME Suite and identified 21 matching motifs in the HOCOMOCO Mouse (v11 CORE) database, corresponding to binding sites of the following proteins: ZFP42, NR4A1, MEF2A, MEF2C, MEF2D, TAF1, YY1, PRDM9, SRF, FOXL2, OTX2, ATOH1, IRF3, NKX6-1, NEUROD2, NEUROD1, ISL1, ALX1, SOX10, ELF5, MSGN1, listed by descending significance. ZFP42/REX1 is an epigenetic remodeler and a marker of pluripotency [64]. PRDM9 is also an epigenetic remodeler. Other proteins are transcription factors essential for muscle development (MEF2A, MEF2C, MEF2D), nervous system development (NEUROD1, NEUROD2, SOX10, OTX2, ATOH1), mesoderm differentiation (MSGN1), keratinocyte differentiation (ELF5), development of craniofacial structures (ALX1), cell morphogenesis (SRF), cell cycle regulation (Taf1), gonad differentiation (FOXL2), type I interferon response (IRF3), development of insulin-producing beta cells and insulin secretion (NKX6-1, ISL1), and stress response cell differentiation, apoptosis, proliferation, inflammation, and metabolism (NR4A1). Although binding sites of these proteins matched the significant VMR motif with p < 0.05, both q-value and E-value calculated to adjust for multiple comparisons were >0.05 for every protein. ZFP42/REX1 binding site had the highest significance of match (p = 4.21 × 10^−4^) with highest enrichment of VMRs (90.5%).

3.12. Overlap of VMRs with Regions Escaping Methylation Reprogramming

Given that VMRs are enriched with ZFP42/REX1 binding motifs and that ZFP42/REX1 plays a role in imprinting, we hypothesized that VMRs may represent regions related to imprinting phenomena. To test this hypothesis, we first compared the overlap of genomic coordinates of VMRs with coordinates of imprinted genes in mice and rats. The lists of imprinted genes were taken from the Geneimprint database and included 159 mouse and 55 rat genes. No overlap with sperm VMRs in either species was identified. Next, we compared the coordinates of mouse sperm VMRs with regions that escape epigenetic reprogramming during early embryogenesis. Regions that escape demethylation were identified as regions with >15% methylation in blastocyst and regions that escape remethylation were identified as regions with <15% methylation in both spermatozoa and embryonic stem cells. Our analysis demonstrated that VMRs are not enriched with either type of escapees.

4. Discussion

Thus, the goals of the current study were two-fold: to test if changes in sperm DNA methylation induced by aging and stressful factors occur in VMRs associated with developmental genes as predicted by SEV hypothesis and to develop an evolutionary and molecular understanding of the nature of DNA methylation changes induced by aging and stressors. Indeed, we identified sperm VMRs as the material basis for DNA methylation changes in response to a broad range of stressors and aging. We demonstrate that the genomic context of DNA regions determines the response of sperm methylome to various factors rather than the type of stressor and/or timing of these factors. Because methylation variation in VMRs increases with age, sperm methylome response to stressful factors may be characterized as the acceleration of epigenetic aging. These findings are supported by a recent report that environmental exposures (heat stress and cadmium) result in mouse sperm methylome changes concordant with changes induced by aging [65].

These findings provide an evolutionary framework for the understanding of epigenetic inheritance via male germline. VMRs in mice and rats are conserved and are associated with genes enriched for developmental categories. These characteristics match predictions from the SEV hypothesis, suggesting that VMRs represent an evolutionary adaptation that increases the chances of survival of the genetic lineage in conditions of a changing environment [43]. VMRs may represent stochastic switches in branching points of developmental trajectories of Waddington’s epigenetic landscape [66]. That view of epigenetic changes in sperm suggests that phenotypic changes in offspring conceived by older fathers or fathers exposed to stressful environments are characterized by increased phenotypic variation, rather than by directional shift in phenotypic traits.

At the molecular level, age-dependent changes in sperm DNA methylation may have a similar nature in somatic cells. Indeed, it was demonstrated that many cytosines increase in methylation variation in human somatic cells with age [67,68,69,70,71] and VMRs were identified in human whole blood samples [72] and mouse and human livers [43], mouse blood [73], and human brain [43]. It was also shown that human blood VMRs are enriched in enhancer regions [68]. Recent single-cell methylation analysis using mouse stem cells demonstrated that major age-dependent changes in DNA methylation are characterized by a stochastic increase in methylation of initially hypomethylated regions and a decrease in hypermethylated regions [73]. It was also shown using purely simulated data that accumulation of stochastic variation is sufficient to build epigenetic clocks that capture age-dependent shifts in DNA methylation [74]. Additionally, changes in bulk blood DNA methylation throughout functional aging (3–35 months of age) in mice were small, with the largest change around 30% per CpG [73]. This observation is concordant with our model (Figure 2G), suggesting that age or stress-dependent changes in DNA methylation are a result of increasing asymmetric variation across a population of cells, rather than the result of directional and coordinated change. Indeed, if the initially hypomethylated region undergoes a coordinated directional shift in methylation, the maximum shift in methylation of approximately 100% may be achieved when all cells gain methylation. In the case of methylation shift due to increasing variation, the maximum methylation shift may only reach 50%—the point when variance is maxed out (max entropy, min information).

We tried to elucidate if VMRs may be involved in the intergenerational transfer of epigenetic information as suggested by the stochastic epigenetic variation hypothesis. We observed that mouse sperm VMRs are enriched with ZFP42/REX1 motif (Table S5). This gene is highly expressed only during spermatogenesis (in spermatogonia and spermatocytes) and in early embryo, specifically in placenta [64]—extra-embryonic tissue mostly controlled by paternally expressed imprinted genes [75]. Although ZFP42/REX1 is not expressed in the late embryo, its gene dosage during early embryogenesis is critical for the survival of the late-stage embryos and neonates, suggesting that during spermatogenesis ZFP42/REX1 creates some epigenetic signature that instructs embryo development at later stages [64]. This assumption is supported by the fact that paternal inheritance of deactivated ZFP42/REX1 causes higher levels of lethality in embryos than maternal inheritance [64]. It is also the central negative regulator of X chromosome inactivation during early female embryogenesis [76] and a central positive regulator of X chromosome reactivation [77]. ZFP42/REX1 also plays a role in genomic imprinting and is hypothesized to bind unmethylated allele of imprinted genes and protect them from DNA methylation removal during embryogenesis [64].

The association of VMRs with ZFP42/REX1 motif suggests that these regions may survive methylation reprogramming in the early embryo and affect developmental trajectories. However, we did not find any significant overlap between VMRs and methylation reprogramming escapees. Thus, it remains unclear if DNA methylation in sperm VMRs has a higher propensity to contribute to the developmental trajectories in offspring than methylation in non-VMRs. It is important to mention here that the ability of sperm methylation patterns to affect initial stages of embryo development before being reprogrammed during blastocyst is poorly understood. According to the classical view, the first transcription in newly formed embryos starts in the eight-cell stage, around 3 days after fertilization [78]. Recent research demonstrated, however, that embryonic transcription initiates as early as the one-cell stage in humans [79] and mice [80], and that zygotic transcription activation starts from the paternal genome [81]. These reports taken together suggest that paternal epigenome may play an important role in determining the developmental trajectories of the embryo before the global wave of methylome reprogramming.

We also attempted to identify specific features of VMRs that may be causatively linked to their high methylation variation. Although VMRs were differentially enriched with some repetitive elements and histone modifications, the fact that both modified histones and repetitive elements overlap with only a small portion of VMRs suggests that these features cannot explain causatively the stochastic nature of DNA methylation in these regions. Surprisingly, VMRs were depleted with open chromatin regions in male germ cells, indicating that these regions are bound to some biological molecules (Table 1). However, they were mostly depleted with histones as well (Table 2), suggesting that molecules other than histones protect VMRs from hyperactive transposase that is the central element of ATAC-seq technology. While the analyses of DNA regions overlap with histones provide a meaningful comparison of VMR and non-VMR associations with various histone marks in spermatogenic cells, caution is warranted when comparing these associations across different cell types. This is due to the fact that ChIP-seq data for different cell types were sourced from distinct studies, each employing varying sequencing depths.

We hypothesize that association with ZFP42/REX1 may play a role in the stochastic instability of methylation in VMRs. Indeed ZFP42/REX1 gene has been duplicated from YY1 in the Eutherian lineage [82] and both transcription factors still share similar binding motifs. The two proteins competing for the same binding site may have different effects on DNA methylation of corresponding regions. For example, ZFP42/REX1 was shown to protect from methylation the paternal allele of imprinted Peg3 gene [64], while YY1 promotes methylation of Peg3 during oogenesis [64]. Additionally, ZFP42/REX1 itself may have a dual effect on DNA methylation of regions to which it is bound. On one hand, it may shield them from methylation; on the other hand, some evidence indicates that ZFP42/REX1 can recruit de novo DNA-methyltransferase DNMT3b [83] and increase Dnmt1 expression [64]. Thus, these two mechanisms—competition between ZFP42/REX1 and YY1 for the same binding sites and their opposite effects on methylation, as well as the potentially dual role of ZFP42/REX1 itself—may be responsible for high methylation variation in VMRs. This hypothesis is to be tested experimentally in future research.

This study has several limitations to be addressed in future research. For example, in our study VMRs were identified based on the analysis of intra-individual methylation variability, where methylation value for each individual organism was averaged for millions of spermatozoa. Future single-cell DNA methylation analysis is needed to test if the same regions are identified as VMRs at the inter-individual level. Functional enrichment analyses in this study were done for genes closest to VMRs or DMRs, and thus they do not account for distant epigenetic regulatory elements.

Future analysis of sperm DNA methylation across a broad age range in rodent models and humans using whole genome bisulfite sequencing may produce a better understanding of age-dependent methylation dynamics and provide a detailed view of the overlap of VMRs, DMRs and various genomic elements, including repetitive elements, open chromatin regions, binding motifs, and others. Additional research is needed to evaluate conservation of VMRs across different cell types and dissect experimentally the molecular mechanisms responsible for unstable methylation in VMRs.

5. Conclusions

This study demonstrates that a broad range of stressors and aging induce DNA methylation changes in sperm in specific regions characterized by high stochastic variability of methylation levels. These findings are concordant with the stochastic epigenetic variation hypothesis, which proposes that VMRs result from adaptive evolution by increasing the chances of genetic lineage survival in changing environmental conditions. This hypothesis explains observed association of VMRs with morphogenetic genes. Our data suggest that high level of stochastic variation in DNA methylation within VMRs may be causally linked to ZFP42/REX1 protein, an epigenetic remodeler involved in genomic imprinting. The nature of stochastic instability of VMRs, as well as their ability to deliver epigenetic information to offspring, requires additional research.