DNA methylation dynamic in male rat germ cells during gametogenesis

Mechanisms of epigenetic regulation include DNA methylation (DNAme) at the carbon-5 of cytosine in CpG dinucleotides, post-translational modification of histones (PTMs) and localization of non-coding RNA to the chromatin (Wu and Morris, 2001). Each of these has been investigated in relation to the regulation of gene expression ultimately driving cell fate. More specifically for DNAme, its presence in gene-regulating sequences, such as promoter and transcription start sites, can determine whether the gene is expressed (hypomethylation) or silenced (hypermethylation) (Wu and Morris, 2001). Hypermethylation within the gene body may also be associated with gene expression (Aran et al., 2011). In addition, some data have suggested a causal link between DNA hypomethylation and genetic instability (Costello and Plass, 2001; Feinberg and Tycko, 2004). In mammals, zygote formation and germline establishment are the two main developmental windows during which the resetting of DNAme (i.e. epigenetic reprogramming) occurs (Cantone and Fisher, 2013). In the germline, there is a near complete erasure of DNAme, followed by a sex-specific de novo DNAme. This establishes germline-specific gene expression profiles but also, in the longer term, a key component of the sperm epigenome necessary for embryo development (Edwards et al., 2017; Wu et al., 2015). Epigenetic reprogramming is sensitive to in utero environment and long-term DNAme changes in spermatozoa have been demonstrated depending on maternal environment (Wu et al., 2015). However, the overall understanding of these effects remains limited by the absence of a mapping of these methylation dynamics throughout spermatogenesis.

Initial efforts to describe the ontogeny of DNAme in male germ cells came from studies in mice in which the erasure of DNAme was observed in migrating primordial germ cells (PGCs) and shown to occur in two temporally distinct waves (Hajkova et al., 2002; Ly et al., 2015; Seisenberger et al., 2012). The first wave of demethylation, starting at gestational day (GD) 8, is observed genome-wide (CpG islands, introns, exons, promoters and intergenic regions), whereas the second wave, starting at GD10.5, occurs in differentially methylated regions (DMRs) of imprinted genes and CpG-rich promoters of male germ cell-specific genes. In male mice germ cells, DNAme reaches its lowest level at GD13.5 (Seisenberger et al., 2012), yet some sequences resist demethylation. These include intracisternal A particles (IAPs) neighboring CpG islands and some non CpG-rich promoters (Guibert et al., 2012; Lane et al., 2003; Popp et al., 2010; Niles et al., 2011). Thereafter, a genome-wide de novo DNAme occurs between GD14 and GD19 in the mitotically arrested gonocytes affecting all DMRs in paternal imprinted genes (Kato et al., 2007; Lees-Murdock et al., 2003). Most repeats also reach a high methylation level before birth; however, some minor and major satellite sequences only reach an intermediate DNAme level of ∼40% (Kato et al., 2007). In addition, some retrotransposons are known to only start de novo DNAme at GD17.5 and be completed after birth by postnatal day (PND) 2, during a second minor phase of de novo DNAme by a mechanism implicating the piRNA pathway (Aravin et al., 2008). DNAme is considered completed in postnatal spermatogonia (Edwards et al., 2017), yet Oakes et al. described both gain and loss of DNAme in a subset of loci from spermatogonia to early meiotic prophase I (Oakes et al., 2007). In addition, a transient genome-wide reduction of DNAme in spermatocytes, at the beginning of the S phase, was described as the source of retrotransposons reactivations (Gaysinskaya et al., 2018).

There are still significant gaps in the knowledge of DNAme dynamics during spermatogenesis and spermiogenesis, in particular. Yet, the detailed mapping of the DNA methylome in mammalian male germ cells is urgently needed to help understand the genotoxic effects of environmental chemicals on the male germline. It is now recognized that in utero exposure to environmental pollutants such as endocrine disrupting chemicals (EDCs) can impair the sperm epigenome, ultimately affecting the resulting embryo development and health of the progeny (Wu et al., 2015). Accumulating data show that EDCs can affect the epigenetic reprogramming (Feil and Fraga, 2011; Pacchierotti and Spanò, 2015). This has been suggested as the link between early-life exposure and outcomes in later life, and even transgenerational inheritance (Vaiserman, 2014; Wu et al., 2015; Ben Maamar et al., 2021). In addition, evidence is emerging that the environment can influence the sperm epigenome during the period of pre-conception in adults. These include diet, environmental chemicals, alcohol and cigarette smoking, amongst others (Bedi et al., 2022; Jenkins et al., 2017). Yet, the remaining important gaps in knowledge of the ontogeny of the male germline epigenome significantly limit our comprehensive understanding of EDC-induced epigenetic changes. More importantly, this fundamental knowledge mostly comes from mouse models, whereas the majority of EDC experiments have been carried out in rats. This difference in model is crucial, considering the differences in metabolism and sensitivity to toxicants (that is more representative of large mammals) in the rat. Only a few studies using immunofluorescence, pyrosequencing of parental imprints or Methylated DNA Immunopreciptation (MEDIP)-based genome-wide study have given a broad picture of the wave of epigenetic reprogramming and later changes in DNAme in rat male germ cells, from PGC to spermatozoa (Ben Maamar et al., 2022; Rose et al., 2014; Rwigemera et al., 2017), but a detailed map is still lacking.

Here, we contribute a comprehensive reference map of DNAme dynamics during spermatogenesis in rat, that will be instrumental for future assessment of epigenetic perturbations of the male germline throughout life. We combined our abilities to sort male germ cells by fluorescence-activated cell sorting (FACS), particularly during the perinatal developmental stages and spermiogenesis, with methyl-seq capture to investigate DNAme dynamics from fetal life to the production of spermatozoa and, for the first time, the extent of DMRs in the differentiating spermatids. Sensitive developmental stages were identified as well as classes of genomic regions seemingly more dynamic than others and potentially more sensitive to perturbations.

Using transgenic rats expressing GFP specifically in germ cells, we purified male gonocytes at GD16, GD18 and GD20 and spermatogonia at PND5 using FACS. We had previously shown the effectiveness of this method by expression microarray with the enrichment of the germline genes (Rwigemera et al., 2021). Based on ploidy and DNA compaction, we also purified four populations of haploid spermatids. Enrichments of step-specific spermatids were evaluated by visual count under a fluorescent microscope (Table S1) and further confirmed, by immunoblots, based on the presence of the phosphorylated form of histone H2AX (γH2AX) in round spermatids R1-8, transition protein 2 (Tnp2) in late spermatids R13-14, and protamine 1 (Prm1) in condensed spermatids R15-17 (Fig. S1). Finally, we collected spermatozoa from the epididymis. In total, we obtained nine populations of germ cells from fetal stages to spermatozoa, each replicate of the adult stages coming from the same individual. Rat methyl-seq analysis was used to evaluate the DNA methylome at each stage of development (n=4/stage).

Sequencing data were aligned against the Rattus_norvegicusRnor_6.0 assembly (Krueger and Andrews, 2011). The average number of paired-end reads uniquely aligned were 41.1 million reads per sample. The total number of CpG sites successfully sequenced were 6,772,950, including 1,301,041 CpGs that were sequenced with enough coverage over all samples. These CpGs were further studied in genomic segments generated by SMART2, an entropy-based framework, initially developed to compare DNAme marks in 42 different human cells (Liu et al., 2016). The 87,290 segments generated had an average length of 445 bp with ≥5 CpGs and displayed a random chromosomal distribution except for an absence of Y chromosome data (Fig. 1A). Based on average methylation level of all segments per sample, the distribution of each experimental replicate was evaluated using a principal component analysis (PCA) that showed clear clustering of replicates per stage and significant separation between each perinatal stage, away from the spermatids/spermatozoa cluster (Fig. 1B). This high variation between perinatal stages was expected as they were sampled during the epigenetic reprogramming phase, while the fact that spermatids/spermatozoa formed a unique cluster suggested high similarity among those samples. Further analysis of the distribution of DNAme levels between developmental stages showed dynamic profiles during perinatal stages up to PND5 followed by relatively stable levels during spermiogenesis up to spermatozoa (Fig. 1C). More specifically, DNAme reached its lowest level at GD18 and then gradually increased until PND5, with the pattern of DNA methylation being almost similar between PND5 and spermatozoa. It is interesting to note that the distribution of DNAme showed two distinct groups, one having low levels of DNAme remaining stable over time, the other showing dynamic levels during the perinatal stages and a stabilization of DNAme levels in the later stages.

To control the quality of our dataset, we investigated DNA features with known specific DNAme profiles, such as imprinted genes, gene body and CpG islands. First, 23 imprinted genes (germline differentially methylated regions; gDMR) described in the mouse genome (Legault et al., 2021) were located and validated for proximity to the expected genes in the rat genome. Among these rat homologs, 13 were found in our dataset and exhibited the expected dynamics of DNAme levels (Table S2). More specifically, all exhibited decreased levels between GD16 and GD18 followed by de novo DNAme until PND5 for paternal imprinted genes such as H19 and Meg3, whereas the maternal imprinted genes such as Snrpn and Igf2r remained hypomethylated at all stages (Fig. 1D). DNAme profiles were also inspected in gene bodies and CpG islands where all expectations were met. Briefly, the drop at the transcription start site (TSS) followed by an increase along the gene (Fig. S2A) or the low level of DNAme within CpG islands (Fig. S2B) were conserved between stages even though the global level of DNAme varied.

Out of the 87,290 segments, 30,103 were stable and 57,187 were dynamic, i.e. 30% difference in DNAme levels and a false discovery rate (FDR)≤0.05. Interestingly, although these two subgroups included DMRs distributed over the entire genome (Fig. 2A), their chromosomal distributions were significantly non-uniform compared with the all-segment dataset (stable: D=0.38657, P-value<2.2e-16; dynamic: D=0.14636, P-value<2.2e-16). Further annotations of the three datasets by HOMER indeed illustrated different proportions of genomic features (Fig. 2B). When compared with all segments, stable segments showed enrichment in promoters and 5′ untranslated regions (UTR), whereas the dynamic segments displayed high proportion of introns (Fig. 2B). Finally, gene ontology (GO) enrichment analysis showed specific and distinct biological process (BP) terms for each group (Fig. 2C and data deposited in GEO under accession number GSE220440). The stable segments showed enrichment for 23 BP terms, including seven implicated in embryogenesis and morphogenesis and 33 enriched in the dynamic segments, including 14 related to immune response and immune cell differentiation and activation. To follow the dynamic progression of DNAme specifically from one stage to the other, further chronological case-control analysis of DMR (≥30%) was carried out, where the control was the stage before the current case (Table 1).

Analysis of specific genomic features showed the DNAme dynamics over time in CpG islands (Fig. 3A), non-coding (nc) RNA (Fig. 3B) and transposable elements (Fig. 3C). Hierarchical clustering heatmaps confirmed that most CpG islands remain hypomethylated at all stages, although we identified a subset of 17 genic CpG islands and 28 intergenic CpG islands showing reprogramming (Fig. 3A, top cluster) (see data deposited in GEO under accession number GSE220440). ncRNA present in our dataset were mostly microRNA (miRNA) going through reprogramming during perinatal stages (Fig. 3B, two bottom clusters). Interestingly, we identified a cluster that remained hypomethylated and another showing late reprogramming only in spermatids (Fig. 3B, two top clusters). Finally, hierarchical clustering heatmaps of transposable elements in subgroups of LINEs, SINEs and LTRs showed that most of them are reprogrammed by PND5 (Fig. 3C, top cluster). Nevertheless, one cluster remained hypomethylated over time in each of these three subgroups (Fig. 3C, bottom cluster). Word clouds were generated to identify the frequencies of retrotransposon families revealing similar distributions in both clusters (Fig. S3).

In order to better describe the changes occurring during perinatal life, we took advantage of our chronological case-control analysis of DMRs, which evidenced that all DMRs were hypomethylated between GD16 and GD18, whereas most, if not all, DMRs were hypermethylated from GD18 to GD20 to PND5 (Table 1).

We identified 282 hypomethylated DMRs between GD16 and GD18 (≥30%) representing the last sites of DNA demethylation (Table 1). Annotations of these regions revealed enrichment in intergenic regions (Fig. 4A). Furthermore, GO terms enrichment analysis identified four BP, all related to cell movement (Fig. 4B). DNA methylation levels are further detailed for the 11 genes related to the top enriched GO term at GD16 and GD18 for all replicates (Fig. 4C).

We identified 44,477 and 46,661 hypermethylated DMRs between GD18 and GD20 or GD20 and PND5, respectively, (≥30%) representing sites of de novo DNAme (Table 1). Comparing these two lists revealed three groups of segments displaying different de novo DNAme kinetics (Fig. 5A). The first was named ‘early’ and consisted of 8,083 segments for which de novo DNAme started at GD18 and was completed by GD20. The second group, named ‘continuous’, contained 36,401 segments for which de novo DNAme started at GD18 and was completed by PND5. The third group was named ‘late’ and consisted of 10,275 segments for which de novo DNAme mostly started after GD20 and was completed at PND5. Interestingly, annotations analysis showed that the first two groups were highly enriched in introns, whereas the ‘late’ group was highly enriched with intergenic regions (Fig. 5B), suggesting that de novo DNAme did not occur randomly but that the kinetics of de novo DNAme vary according to specific genomic features. Furthermore, analysis of GO terms BP showed both group-specific enrichment and BPs common to all three groups (Fig. 5C and data deposited in GEO under accession number GSE220440). It is interesting to note that the late group has more unique BPs than the other groups, suggesting that de novo DNAme kinetics also vary depending on BP.

Because we observed shared BPs among groups, we questioned the similarities between the genes involved in a similar BP (Fig. 5D). This revealed that some genes are unique to each group but other are shared among groups. Interestingly, although the average percentage of CG sites was similar in all these genes, their sequence length and number of exons differed. Indeed, genes common to the three groups were four times longer (157,630 bp±7486) compared with the ones unique to the early (36,314 bp±2341), continuous (48,860 bp±1418) or late (38,270 bp±3586) groups; they also contained significantly more exons (17.8±0.50, 9.94±0.32, 10.8±0.15 and 8.2±0.21 in the common, early, continuous and late groups, respectively: P<1.10⁻³⁶). The fact that some genes are common to the three groups suggested that a different de novo DNAme kinetic could apply on the same gene depending on where the DMRs is located. This was further confirmed by using the 777 genes common to all groups, as the analysis of DMRs associated with these genes do not overlap with each other (Fig. 5E). Furthermore, these DMRs exhibited the expected methylation progression specific to their group, as presented in the heatmaps of DNA methylation (Fig. S4). Using these lists of DMRs located on the same genes, we attempted to verify whether de novo DNAme kinetics varied logically according to location within genes and the reading frame of the gene. Analysis of location compared with the TSS or along the genes failed to identify a specific pattern.

To further deepen our understanding of the DNAme dynamic during spermiogenesis, we also used the chronological case-control analysis of DMRs but relaxed our threshold to 10% change in DNAme between steps (Table 2). This made our analysis more permissive, revealing finer changes in the genome, showing that spermiogenesis still displays DNAme dynamics, but with a lower amplitude compared with perinatal stages. Although thousands of DMRs were identified, the highest number was found when comparing elongating spermatids R9-12 to late spermatids R13-14 (Table 2). Hierarchical clustering heatmap analysis of all DMR illustrated this dynamics, particularly between R9-12 and R13-14 (Fig. 6A). Interestingly this analysis revealed that the two most similar groups were the round spermatids R1-8 and epididymal spermatozoa (data not shown). Clustering analysis using each stage in a chronological order revealed eight clusters of DMRs (Fig. 6A). These include clusters 3 and 4 for which we observed hypermethylation in R9-12 elongating spermatids and clusters 7 and 8 displaying specific hypomethylation in R9-12 elongating spermatids and R13-14 late spermatids respectively. Annotation of these four clusters indicated that clusters 3 and 4 were enriched in intergenic regions while clusters 7 and 8 were highly enriched in introns (Fig. 6B). Further, GO term analysis showed cluster-specific enrichment for different BP with a large proportion related to neurodevelopment (Fig. 6C).

DNAme is the first epigenetic mechanism to have been demonstrated in epigenetic inheritance by male germ cells (Roemer et al., 1997; Surani et al., 1993; Holliday, 1990). However, there are still many gaps in the understanding of the establishment of DNAme in male germ cells, both during epigenetic reprogramming occurring during perinatal life and spermiogenesis (Edwards et al., 2017). This study contributes a comprehensive reference map of DNAme dynamics during spermatogenesis. Using FACS approaches and a DNA capture strategy, we detailed the precise timing and genomic locations of DNAme variations during two key developmental windows of spermatogenesis. Data obtained confirmed and refined the understanding of the global dynamics of perinatal epigenetic reprogramming and described both losses and gains of DNAme during spermiogenesis, at very specific steps of chromatin remodeling. These waves of DNAme variations could represent specific window of sensitivity to environmental influence in restricted dynamic genomic regions leading to specific outcomes.

With the advancement in cell purification approaches and low-input next-generation sequencing, the study of the DNA methylome in germ cells has been facilitated, yet it suffers from the lack of a standardized platform for the species of interest. The Illumina platforms have been the workhorses used to investigate DNAme in human studies, but their designs are still pending for other species. As a consequence, other methods have been developed (McGraw et al., 2013) including DNA capture. Such a platform has been developed for the rat, targeting CG-rich sequences mostly in RefSeq genes to study the impact of chronic stress exposure on blood sample DNAme (Carey et al., 2018). More specifically, this technique allows the targeting of 2.3 million randomly distributed CpGs, with 50% of the generated data representing gene regions, even though they represent only 2% of the whole genome. In our hands, this assay was highly reproducible, with low DNA input generating high quality sequencing data targeted towards a ‘functional genome’ centered around coding genes. Here, we studied proliferative and quiescent gonocytes, spermatogonia, spermatids at different steps of chromatin compaction and spermatozoa obtained with high purity. Absence of somatic contamination was validated by the expected profiles of DNAme in imprinted genes (Reik and Walter, 2001). In parallel, the quality of the sequencing data could be confirmed by observations of the expected DNAme patterns in CpG islands (Antequera and Bird, 1993) and gene bodies (Aran et al., 2011) at all developmental stages studied.

The overall DNAme dynamics along developmental stages were studied using a threshold of 30% variation in DNAme levels, revealing that 34.5% of the genomic segments analyzed had stable methylation level across developmental stages, whereas 65.5% of the segments were dynamic. Strikingly, these two groups represent distinct hotspots in the genome and cover genes involved in completely different BP. Stable segments were not considered to be DMRs; however, they could be those with methylation changes falling below the threshold. They could also have lost DNAme before GD16 and remain hypomethylated over time or be hypomethylated throughout gametogenesis since germ cell specification. The latter would include CpG-rich promoters and maternal imprinted genes that are not subjected to epigenetic reprogramming because they are unmethylated at all stages of spermatogenesis (Edwards et al., 2017). Our data are consistent with this, as the stable group of segments was highly enriched in promoters. Interestingly, we also identified a fraction of unexpected ncRNA that remained hypomethylated and a subset of CpG islands with DNAme reprogramming profiles. In parallel, among the dynamic segments, a high enrichment in intronic regions was observed suggesting that gene expression control could originate from these genomic regions. Moreover, we identified most of the targeted ncRNA and repetitive sequences to be reprogrammed. In mice germ cells, there is a peak of expression of IAP, LINE and SINE at GD13.5 concomitant with DNAme erasure (Jansz, 2019). In this study, most of the retrotransposons, SINE, LINE and LTR, studied were intronic, which could significantly impact gene expression.

Epigenetic reprogramming occurs in the germline during perinatal life, starting in migrating PGCs when DNAme is erased followed by a resetting until around birth in males (Edwards et al., 2017). It is generally considered that these global changes are conserved across mammalian species (Chen and Riggs, 2011), but the extent to which the exact timing and genomic features are preserved remains largely unknown. Here, to detail the ontogeny of this wave in the male rat, we collected four key, post-migratory stages of germ cell development: GD16, during the fetal proliferation phase; GD18, upon entry into quiescence; GD20, during quiescence; and PND5, after mitosis resumes in differentiating spermatogonia (Culty, 2009). Previous studies by us and others had suggested that de novo DNAme started only at GD18 in the rat (Rose et al., 2014; Rwigemera et al., 2021, 2017). The present data confirmed that the lowest DNAme levels were observed at GD18, identifying the last coding regions to be demethylated between GD16 and GD18 to be highly enriched in genes involved in the negative regulation of cellular movement. PGC migration to the gonad has been described to end at around GD13 in rats (Kemper and Peters, 1987; Culty, 2009), therefore, one could hypothesize that the genes involved in inhibiting the migratory characteristics of gonocytes are the last to be demethylated because their regulation was still needed up to that stage. Of note, our data showing the lowest value at GD18.5 also suggest the demethylation phase lasts longer in male rats than mice, as the lowest level of DNAme in mice has been described at GD13.5, which would correspond to GD15.5 in rat (Culty, 2009). This is of importance when considering windows of sensitivity and species-differences depending on the timing of in utero exposure experiments.

De novo DNAme driven by the specific timely expression of DNA methyltransferase 3L (DNMT3L) (La Salle et al., 2004; Rwigemera et al., 2021) before birth in gonocytes allows genomic imprinting and proper cell differentiation. The onset of de novo DNA methylation occurs at GD14 in mouse, immediately before gonocytes enter quiescence, and is mostly completed before birth at GD19 (Ly et al., 2015). Methylation of some retrotransposons has been shown to have a delayed onset, occurring at GD17.5 and completing after birth. Our data are consistent with such patterns starting in rat quiescent gonocytes at GD18 or GD20 with differences in the dynamics of remethylation, the final pattern being acquired in PND5 spermatogonia. Interestingly, we showed that the genomic regions that gain methylation early are mostly introns. Because DNAme in introns has been suggested to control gene expression (Aran et al., 2011), it would be interesting to investigate whether the genes included in this early group are crucial for germ cell survival and differentiation. On the contrary, the genomic regions to start de novo DNAme only after GD20 are mostly intergenic. Together, these data strongly suggest that gain of methylation is a non-random process, first starting within genes. This also implies that, depending on the timing of exposure, xenobiotics could impair DNAme in distinct gene features leading to different effects. In parallel, we have identified DMRs following the three different kinetics in the same genes. This applies to long genes, containing many exons. Yet, we could not identify specific kinetics of de novo DNAme along the genes based on their location. Further analysis will be required to identify what guide the different dynamics of de novo DNAme based on specific sequences or other epigenetic marks.

DNAme levels in post-pubertal germ cells during spermatogenesis are considered to be globally stable from spermatogonia to spermatozoa (Edwards et al., 2017; Ly et al., 2015). However, although several studies suggested variations during meiosis, none has investigated spermiogenesis. Yet, it is the time for a near complete change in chromatin structure followed by compaction (De Vries et al., 2012). This striking change in chromatin structure and DNA topology involves histone variants and post-translational modifications of histones leading to their almost complete eviction from the chromatin before protamination (Goudarzi et al., 2014). Whether those changes are related to DNAme dynamics in the spermatid haploid genome is still unknown. In this study we purified four populations of spermatids according to their differentiation, namely round spermatids (R1-8), when the chromatin is decondensed, elongating spermatids (R9-12), when histones hyperacetylation occurs to promote histone replacement by transition proteins, late spermatids (R13-14), during chromatin condensation, and condensed spermatids (R15-17), during protamination. The identification of DMRs during spermiogenesis became evident when using a 10% DNAme level change threshold. The peculiar dynamics observed between R9-12 and R13-14 could be linked to chromatin condensation changes and the replacement of histones, especially given that DNAme is known to induce or facilitate chromatin condensation. Overall, these results suggest that DNAme could play a role in chromatin remodeling during spermiogenesis, potentially by its interplay with histones (Prakash and Fournier, 2017). Moreover, these dynamics highlight the possibility of alteration by exogenous factors during the pre-conception period. The transient hypomethylation observed between R9-12 and R13-14 mostly within intronic regions was found to involve synaptic genes with high significance. This may represent a sensitive window through which the mature sperm DNAme landscape is established. Any alteration in this sensitive transition by xenobiotics would therefore be expected to first impact the methylation pattern of synaptic genes, with potential cognitive consequences to the offspring (McCarthy et al., 2018). Because transient DNA strand breaks in spermatids were also found to preferentially arise within synaptic genes (Gregoire et al., 2018) the potential interplay between these mechanisms deserves further investigation.

This study generated a genome-wide map of DNA methylation in male rat germ cells from fetal gonocytes to spermatozoa. The selected capture tool generated a solid and reproducible dataset of the rat ‘functional genome’. It includes most of the coding genes, imprinted genes, repetitive elements and ncRNA. As expected, the overall dynamics of epigenetic reprogramming showed patterns very similar to those described in mice. We did, however, identify some species-specific timing, with the demethylation phase lasting longer in male rats. In addition, our data showing different dynamics of de novo DNAme timing according to classes of genomic features could imply windows of sensitivity and potential hotspots more sensitive to perturbations. Finally, the description of DNAme dynamic during spermiogenesis, which had never been done before, strongly suggested that chromatin compaction is also a time of DNAme variation and revealed a new window of sensitivity to environmental influence. This baseline will be crucial for future toxicology studies in rats to identify pre- and peri-conception DNAme alterations and genomic regions susceptible to environmental influence.

Transgenic Sprague-Dawley rats expressing germ cell-specific GFP (GCS-EGFP) (Cronkhite et al., 2005) were used as previously described in Rwigemera et al. (2017). Rats were housed on a 12 h light/12 h dark cycle and were fed with commercial food (Teklad global 18% protein, Envigo) and tap water ad libitum. All animal studies were conducted in accordance with the guidelines set out by the Canadian Council of Animal Care (CCAC) and as approved by the Institutional Animal Care and Use Committee at the INRS (Protocols 1808-01 and 2109-03). Two females were caged with one male overnight and vaginal smears were carried out the following day to identify sperm-positive females. That day was counted as GD0. To acquire pre-birth germ cells, pregnant rats were euthanized by CO₂ asphyxiation and subsequent cervical dislocation at GD16, 18 and 20. Fetuses were removed from the uterus, and gonads were dissected under a binocular microscope; their sex was determined by the morphology of the gonads. To acquire post-birth germ cells, natural birth occurred at GD22.5 which was counted as PND0. Neonates were euthanized by decapitation at PND5 and their gonads were immediately removed for cell dissociation and sorting.

Ten males from the GCS-GFP colony maintained in the Delbes lab were sent for cell sorting by the Boissonneault team. There, the animal study protocol was approved by the Institutional Ethics Committee of Université de Sherbrooke (protocol 2020-2835). Rats were euthanized between PND90 and PND120 by CO₂ asphyxiation and subsequent cervical dislocation. Their gonads and epididymides were immediately harvested and used to collect testicular spermatids and mature epididymal spermatozoa, as detailed below.

At GD16, 18, 20 and PND5, explanted testes were pooled per litter and processed to obtain a testicular cell suspension using a double enzymatic digestion as previously described (Rwigemera et al., 2017). Briefly, testes were slightly cut and immersed in a solution of 1 mg/ml collagenase (Sigma-Aldrich, C9891) and 0.02 mg/ml DNase (Sigma-Aldrich, D4527) and incubated for 15 min at 37°C with manual shaking every 5 min. After a 5 min centrifugation at 500 g, the pellet was re-suspended in 0.25% trypsin–EDTA (Life Technologies, 25200-056) and incubated for 10 min at 37°C, after which 10% fetal bovine serum (FBS) (Life Technologies, 10099-133) was added to neutralize trypsin. After filtration (70 μm strainer; Thermo Fisher Scientific, 352350) and centrifugation (at 500 g, 5 min, 4°C), cells were resuspended in presort buffer (1× Ca²⁺/Mg²⁺-free PBS supplemented with EDTA 1 mM, HEPES 25 mM and 1% FBS) to obtain 5.10⁶ cells/ml and used for cell sorting with the BD FACSJazz (BD Biosciences) at an event rate of ∼2500/s. The use of GCS-EGFP rats enabled purification of germ cells by FACS with sort gates set on forward and side scatter, trigger pulse width and the GFP signal detected (488 nm laser for excitation and emission at 530/40 nm) using the BD FACS software. After sorting, the GFP-positive cell fractions were washed with Hank's Balance Salt Solution (HBSS) and counted. The number of cells collected was calculated using a hemocytometer under an Eclipse Ti-S fluorescent microscope (Nikon), and purity per fraction (percentage of GFP-positive cells; Table S1) was calculated by FACS over the analysis of 1000 purified cells. Sorted cells were centrifuged (at 500 g, 5 min, 4°C), flash frozen and stored at −80°C until DNA extraction.

FACS of GCS-EGFP rat (Cronkhite et al., 2005) spermatids at specific differentiation steps was performed as previously described (Simard et al., 2014, 2015). Briefly, testes were excised, decapsulated, minced, homogenized and filtered on a 40 µm nylon mesh. Testicular cells were fixed on ice in sorting buffer [1× PBS+1 mM EDTA (pH 8.0)+25 mM HEPES (pH 7.0)+1% heat-inactivated FBS] for 15 min with 75% ethanol+0.8 mM EDTA. Fixed germ cells were stained with SYTO 16 Green Fluorescent Nucleic Acid Stain (Thermo Fisher Scientific) and DAPI (Sigma-Aldrich) and sorted on a BD FACSAria™ III cell sorter (BD Biosciences) using DNA staining intensity, forward-scattered light (size, FSC) and side-scattered light (granularity, SSC) to define and gate populations on the FACSDiva software. Four populations of spermatids were isolated: round spermatids (R1-8), elongating spermatids (R9-12), late spermatids (R13-14) and condensed spermatids (R15-17). To collect epididymal spermatozoa, epididymides were minced and spermatozoa were flushed in sorting buffer. Sperm were collected from the supernatant by centrifugation at 1000 g for 10 min and fixed on ice in 75% ethanol+0.8 mM EDTA for 15 min. Spermatids and spermatozoa purity was validated by microscopy (Table S1) and immunoblotting for specific markers of spermiogenesis (Fig. S1).

Proteins were extracted using the acid extraction protocol previously described (Torregrosa et al., 2006) and separated on a 15% acid–urea polyacrylamide gel. Protein detection was carried out using the following primary antibodies: rabbit antihyperacetylated Histone H4 (Penta) antibody (Millipore, 06-946, 1/1000), mouse anti-H2AX-Phosphorylated (Ser139) antibody [3F2] (BioLegend, 613402, 1/3000), mouse antiprotamine P1 (Briar Patch Biosciences, Mab-001 Hup 1N, 1/3000) and goat anti-TNP2 antibody (K-18) (Santa Cruz Biotechnology, sc-21106, 1/500). The following secondary antibodies were used: Alexa Fluor 680 goat anti-mouse IgG (H+L) (Invitrogen, A21057, 1/10,000), IRDye 800CW goat anti-rabbit IgG (H+L) (Li-Cor Biosciences, 926-32211, 1/10,000) and Alexa Fluor 680 rabbit anti-goat IgG (H+L) (Invitrogen, A21088, 1/10,000). Fluorescence was visualized using the Odyssey Infrared Imaging System (model 9120; Li-Cor Biosciences).

Genomic DNA (gDNA) was extracted from FACS-sorted germ cells using the QIAamp DNA micro kit (Qiagen) and the DNeasy Blood & Tissue kit (Qiagen) for spermatozoa. The purification of DNA from sperm was carried out using the DNeasy Blood and Tissue kit; protocol 2, DY03 (Qiagen). DNA were extracted from perinatal stages (125,276±10,915 cells) and adult stages (1.8×10⁶±4.9×10⁵ cells) including a RNase A (Qiagen) treatment after the cell lysis step. Quantity and quality were assessed using a 4200 TapeStation system with Genomic DNA ScreenTape (Agilent). We selected four biological replicates (n=4/stage; n being one pooled litter) from perinatal stages with the best purity and highest amount of gDNA extracted. For adult stages, cell extracts were purified from four different males; each male represented a replicate for which we had obtained all the stages studied from the same individual and adequate quality and quantity of gDNA.

Methyl-seq capture libraries were prepared using the Rat SureSelectXT Methyl-Seq Target Enrichment Panels kit (Agilent) as described in the SureSelectXT Methyl-Seq target Enrichment System for Illumina Multiplexed Sequencing protocol (Agilent, Version E0, April 2018). The SureSelect Target Capture System was designed for the rat genome to capture CpG islands, island shores (±1 kb flanking CpG islands) and GC-rich regions from all RefSeq genes (Carey et al., 2018). Preliminary control studies had confirmed the robustness of library preparation and genome coverage for a range of DNA inputs (1 µg-100 ng). All libraries were therefore prepared with a DNA input ranging from 300 ng to 191.8 ng (Table S1) with minor modifications in accordance with the Agilent SureSelectXT Methyl-Seq Applications with Low-Input DNA and Smaller Capture Libraries protocol (PR7000-0620, January 2017, 5991-7838EN): the post-capture PCR cycles were increased. In brief, gDNA was diluted in 1× Low TE buffer (Thermo Fisher Scientific) in a total volume of 50 µl and sheared with a M220 Focused-Ultrasonicator (Covaris). DNA ends were repaired and samples purified using SparQ PureMag Beads (QuantaBio). Quality analysis of end-repaired DNA was carried out using High Sensitivity D1000 Screen Tape and 4200 TapeStation (Agilent). DNA fragments were 3′ end dA-tailed and purified with SparQ PureMag Beads. Immediately after, the methylated adapters were ligated, purified with SparQ PureMag beads and analyzed with High Sensitivity D1000 Screen Tape on 4200 TapeStation (Agilent). The volume of each adaptor-ligated DNA sample was concentrated using an Eppendorf Vacufuge Concentrator (Thermo Fisher Scientific) to a final volume of 3.4 µl and the sample was hybridized with SureSelect rat Methyl-Seq Capture library mixture for 16 h at 65°C in a C1000 Touch Thermal Cycler (Bio-Rad) with a heated lid at 105°C. Dynabeads Myone Streptavidin T1 (Thermo Fisher Scientific) was prepared and used to capture hybrids during 30 min at room temperature. The beads were then successively washed in wash buffers 1 and 2 (Agilent Technologies). Capture DNA was submitted to a bisulfite conversion using reagents from the EZ DNA Methylation-Gold kit (Zymo Research). The bisulfite-treated DNA was PCR amplified and the numbers of PCR cycles were increased in accordance with the alternative method for low-input DNA amounts (Agilent). More precisely, the number of PCR cycles was modified depending on the concentration of the DNA at the methylated adapter ligation step: 11, 12, 14 and 16 PCR cycles were performed for 2000-3000 pg/µl, 1000-2000 pg/µl, 500-1000 pg/µl and <500 pg/µl input, respectively. Amplified bisulfite-converted DNA was purified using SparQ PureMag Beads (QuantaBio), followed by an amplification reaction to add their unique index-tag as described in the SureSelect Methyl-Seq protocol (Agilent) and purified with beads. Indexed libraries were assessed for their quality using High Sensitivity D1000 Screen Tape and 4200 TapeStation (Agilent) and pooled at the same equimolar amount of 5 ng for each library. The 36 methyl-seq capture libraries were sequenced in paired-end on NovaSeq6000 S4 system (Illumina). Sequencing was performed at Centre d'expertise et de services Génome Québec in Montreal (Canada).

All sequencing data used in this publication have been deposited in the NCBI Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession number GSE220440. Quality of raw reads was assessed using fastqc (Wingett and Andrews, 2018). Raw reads were then cleaned from remaining adaptor sequences and low-quality nucleotides using Trimmomatics 0.36 (Bolger et al., 2014) and a minimum sequencing quality of 30. Reads longer than 30 bp were then aligned to the converted and unconverted reference genome (http://ftp.ensembl.org/pub/release-104/fasta/rattus_norvegicus/dna/) using Bowtie2 (Langmead et al., 2018) as implemented in the Bismark 0.19.0 pipeline (Krueger and Andrews, 2011). Based on read alignments with converted and unconverted genomes, DNA methylation was extracted using Bismark methylation extractor, followed by a filtration step to obtain data with a minimum of 10× coverage. Then the SMART2 Python package was used for segmentation and determination of DMRs (Liu et al., 2016). The parameter settings for data segmentation were: Euclidean_Distance, 0.2; minimum CpG number per segment, 5; minimum segment_length, 20; CpG distance, 500; similarity entropy: 0.6; Methylation_Specificity: 0.5. The parameter settings for DMR identification were: AbsMeanMethDiffer, 0.3 or 0.1; pDMR, 0.05; and pDMR CaseControl, 0.05. The case-control comparison DMRs were corrected for multiple comparisons FDR≤0.05 by P.adjust R function. The bean plot and the PCA plot were created using SeqMonk software (v1.48.0, Babraham Institute).

The University of California, Santa Cruz (UCSC) Blast tool in the rn6 rat genome assembly was used with default parameters to determine the homolog of 23 imprinted genes (gDMR) described in the mouse genome (Legault et al., 2021). These mouse gDMR were located and validated for proximity to the expected genes in the rat genome using a combination of three criteria: highest score (mean score 1295.5), highest percentage of identity (mean identity 89.6%) and the similarity between the imprinted gene names in mice and rats. The start and stop were those given by the blast tool, and the length was the same as the mouse sequence.

Sequence annotation was carried out using Homer 4.11 (Heinz et al., 2010). The gene identifications linked with gDMRs (including DMRs annotated as 3′ UTR, transcription termination site, introns, exons, promoter, 5′ UTR) in the column ‘annotation’ were then used for GO enrichment analysis using clusterProfiler 4.0 (Wu et al., 2021).

The density distributions of DMRs over the genome were obtained using the circlize R library (Gu et al., 2014). The density distributions of the various DMR subgroups (all, stable, dynamic) were compared using a Kolmogorov–Smirnov test implemented in the ks.test R function. The ComplexHeatmap 2.13.1 package was used for heatmap and clustering analysis (Gu et al., 2016).

This research was enabled in part by support provided by the Digital Research Alliance of Canada (https://alliancecan.ca/en). We thank Hongbo Liu (Institute for Diabetes, Obesity, and Metabolism, University of Pennsylvania, Perelman School of Medicine, Philadelphia, PA, USA) for providing support using the SMART2 tool and performing DMR search analysis. We also thank all the staff of the animal facilities involved in this project.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201606.reviewer-comments.pdf