Potential role of KRAB-ZFP binding and transcriptional states on DNA methylation of retroelements in human male germ cells

  1. Kei Fukuda  Is a corresponding author
  2. Yoshinori Makino
  3. Satoru Kaneko
  4. Chikako Shimura
  5. Yuki Okada
  6. Kenji Ichiyanagi
  7. Yoichi Shinkai  Is a corresponding author
  1. Cellular Memory Laboratory, RIKEN Cluster for Pioneering Research, Japan
  2. Laboratory of Pathology and Development, Institute for Quantitative Biosciences, The University of Tokyo, Japan
  3. Department of Obstetrics and Gynecology, Ichikawa General Hospital, Tokyo Dental College, Japan
  4. Laboratory of Genome and Epigenome Dynamics, Department of Animal Sciences, Graduate School of Bioagricultural Sciences, Nagoya University, Japan

Abstract

DNA methylation, repressive histone modifications, and PIWI-interacting RNAs are essential for controlling retroelement silencing in mammalian germ lines. Dysregulation of retroelement silencing is associated with male sterility. Although retroelement silencing mechanisms have been extensively studied in mouse germ cells, little progress has been made in humans. Here, we show that the Krüppel-associated box domain zinc finger proteins are associated with DNA methylation of retroelements in human primordial germ cells. Further, we show that the hominoid-specific retroelement SINE-VNTR-Alus (SVA) is subjected to transcription-directed de novo DNA methylation during human spermatogenesis. The degree of de novo DNA methylation in SVAs varies among human individuals, which confers significant inter-individual epigenetic variation in sperm. Collectively, our results highlight potential molecular mechanisms for the regulation of retroelements in human male germ cells.

Editor's evaluation

Retrotransposons undergo massive reprogramming of their methylation states during germ cell development, but some elements are immune to this remodeling. This manuscript explores the contribution of binding motifs for KRAB-Zinc Finger Proteins (KZFPs) and position towards genes to explain the variable methylation dynamics of different retrotransposon families, namely L1, SVA and LTR12, as well as potential inter-individual variation during male germ cell development in humans, using an integrative analyses of available sequencing datasets. By bringing insights into the complex regulation of retrotransposons, it could be of particular interest to the epigenetics community.

https://doi.org/10.7554/eLife.76822.sa0

Introduction

Transposable elements comprise more than 40% of most extant mammalian genomes (Lander et al., 2001). Among these, certain types of transposable elements called retroelements, including short/long interspersed elements (SINEs/LINEs) and hominoid-specific retrotransposons SINE-VNTR-Alus (SVA) are active in humans and can be transposed (Huang et al., 2012; Maksakova et al., 2013; Ostertag et al., 2003). As retrotransposons cause genome instability, insertional mutagenesis, or transcriptional perturbation and are often deleterious to host species (O’Donnell and Burns, 2010), multiple defense mechanisms have evolved against transposition. The first line of defense is transcriptional silencing of integrated retroelements by chromatin modifications, such as DNA methylation and histone H3 lysine 9 (H3K9) methylation (Fukuda and Shinkai, 2020; Goodier, 2016). Most retroelement families are bound by Krüppel-associated box domain zinc finger proteins (KRAB-ZFPs), which coevolved to recognize specific retroelement families (Imbeault et al., 2017; Jacobs et al., 2014; Wolf et al., 2015). KRAB-ZFPs repress retroelements by recruiting KAP1/TRIM28 (Sripathy et al., 2006) and other repressive epigenetic modifiers (Schultz et al., 2002; Schultz et al., 2001).

Restricting retroelements is especially important for germ cells, because only germ cells transmit genetic information to the next generation. During embryonic development, primordial germ cells (PGCs) undergo epigenetic reprogramming, characterized by DNA demethylation and global resetting of histone marks in mice and humans (Gkountela et al., 2015; Guo et al., 2015; Kobayashi et al., 2013; Seisenberger et al., 2012; Seki et al., 2007; Tang et al., 2015). A subset of young retroelements resists this global DNA demethylation event in PGCs, which may be required for retroelement silencing (Gkountela et al., 2015; Guo et al., 2015; Kobayashi et al., 2013; Seisenberger et al., 2012; Seki et al., 2007; Tang et al., 2015). H3K9 trimethylation mediated by SETDB1 is enriched in DNA demethylation-resistant retroelements in mouse PGCs (Liu et al., 2014). As SETDB1 regulates DNA methylation of a subset of retroelements (Matsui et al., 2010; Rowe et al., 2013), and it is recruited to the retroelements via interaction with KRAB-ZFPs, it has been hypothesized that SETDB1/KRAB-ZFPs may contribute to DNA demethylation resistance in PGCs.

In contrast to the extensive DNA hypomethylation in PGCs, the genomic DNA of sperm is highly methylated in both humans and mice (Hammoud et al., 2014; Kobayashi et al., 2013; Molaro et al., 2011; Okae et al., 2014). Retroelements are also subjected to de novo DNA methylation during spermatogenesis in mice via the PIWI/piRNA pathway (Aravin et al., 2008; Inoue et al., 2017). Epigenetic alterations in retroelements and dysfunction of retroelement silencing pathways in male germ cells are associated with male sterility linked to azoospermia (Aravin et al., 2007; Bourc’his and Bestor, 2004; Carmell et al., 2007; Heyn et al., 2012; Urdinguio et al., 2015). In addition, epigenetic alterations of retroelements in male germ cells can be potentially transmitted to the next generation with phenotypic consequences (Daxinger et al., 2016; Rakyan et al., 2003). Therefore, deciphering the regulatory mechanisms of retroelements in germ cells contributes to the understanding of sterility and transgenerational epigenetic inheritance. Extensive studies have been conducted to understand DNA methylation mechanisms in mouse spermatogenesis; however, limited progress has been achieved in humans.

In this study, we aimed to clarify the regulatory mechanisms of DNA methylation of retroelements in human germ cells and performed an integrative analysis using three sets of previously reported data, which included whole-genome bisulfite sequencing (WGBS) data for human PGCs (hPGCs) and sperm, the transcriptome of human male germ cells, and comprehensive human KRAB-ZFPs ChIP-exo data.

Results

Transposable elements showing DNA demethylation resistance in hPGCs

To learn more about the factors that contribute to DNA demethylation resistance in hPGCs, we reanalyzed publicly available WGBS data for male hPGCs (Guo et al., 2015). The global erasure of DNA methylation is mostly complete at 19 weeks of gestation (Figure 1A), therefore, we analyzed the DNA methylation status of full-length transposable elements (a copy whose length is 90% or more of the length of the consensus sequence of each subtype, listed in Supplementary file 1) in male hPGCs at 19 weeks of gestation to identify retroelements that showed resistance to demethylation. We generally focused on the retroelement types that had been analyzed for the DNA methylation status more than 30 copies. Among the retroelements we analyzed, the primate-specific retroelement families L1PA, SVA, and LTR12 showed high levels of DNA methylation (Figure 1B). In the SVA family, SVA_A, which emerged 13–14 million years ago (Mya) and is the oldest SVA type, showed the highest DNA methylation levels relative to other SVA types. This includes the currently active SVA_E/F (Figure 1C). In the L1 family, L1PA3–5, which emerged 12–20 Mya and is moderately young, showed higher methylation levels than the older (L1PA5–8) and younger L1 types, including the currently active L1 (L1HS) (Figure 1D). LTR12 (also known as HERV9 LTR) is not currently active and all LTR12 types are highly methylated (Figure 1E). Therefore, it appears that young but inactive L1PA, SVA, and LTR12 types are resistant to DNA demethylation in hPGCs. Because 100 bp short-read NGS data did not map efficiently onto the currently active L1 transposon, L1HS (Figure 1—figure supplement 1A), and DNA methylation of only about 10% of full-length L1HS copies could be analyzed (Figure 1—figure supplement 1B), it is possible that a subset of L1HS is resistant to DNA demethylation. Epigenome analysis using long-read sequence technology, such as nanopore sequencing, may provide an answer to this question (Ewing et al., 2020). Even though some retroelement types showed relatively high DNA methylation levels in hPGCs, the DNA methylation levels of each retroelement type were highly variable among full-length copies (Figure 1C–E), which prompted us to try to identify potential DNA sequences required for DNA demethylation resistance. To this end, we classified each retroelement copy according to DNA methylation levels as follows: low <20%, 20% ≤ medium < 60%, and high ≥60%. Using this classification, we determined that both the ‘high’ and ‘low’ classes of copies exist in highly methylated retroelement types in hPGCs, such as SVA_A, L1PA3, and LTR12C (Figure 1F–H).

Figure 1 with 1 supplement see all
Retroelements showing DNA demethylation resistance.

(A) Violin plots showing DNA methylation levels of each CpG site during human male germ-cell development. DNA demethylation was almost completed at 19 weeks of gestation. (B) Scatter plots showing average DNA methylation level of each retroelement type between somatic cells and male human primordial germ cells (hPGCs) at 19 weeks of gestation. Only full-length copies were used for this analysis, and retroelement types with ≧30 full-length copies were shown. Each plot was colored according to its retroelement family (red: SVA, blue: L1, green: LTR, gray: Alu). (C–E) Violin plots showing DNA methylation level of each retroelement type in hPGCs at 19 weeks of gestation. p-Value was calculated by Tukey’s test and was described in Supplementary file 2. The number in parentheses was analyzed copy number. (F–H) Bar graphs showing the fraction of ‘low’, ‘medium’, and ‘high’ methylated class of each retroelement type in male hPGCs at 19 weeks of gestation. The retroelement copies used in these figures were same as those in (C-E).

The presence of ZNF28- and ZNF257-binding motifs are correlated with demethylation resistance in SVA_A

KRAB-ZFPs are important factors for retroelement silencing. Their activity is mediated by the recruitment of KAP1 and SETDB1, which induces retroelement DNA methylation (Matsui et al., 2010). To investigate whether KRAB-ZFPs could be involved in the DNA demethylation resistance of SVAs, we reanalyzed the binding peak data of 250 KRAB-ZFPs identified by ChIP-exo using exogenously tagged KRAB-ZFPs in human HEK293T cells (Helleboid et al., 2019; Imbeault et al., 2017). We observed that the ZNF257 and ZNF28 peaks overlapped more frequently with highly methylated SVA_A copies than with lowly methylated copies (Figure 2A). Because peaks of ZNF611 and ZNF91, which interact with SVAs in human embryonic stem cells (hESCs) (Haring et al., 2021; Jacobs et al., 2014), were observed in both lowly and highly methylated SVA_A copies (Figure 2A), it is unlikely that these two KRAB-ZNPs contribute to the differences in DNA methylation among SVA_A copies. Of the ‘high’ SVA_A elements, 63.8% and 44.7% were bound by ZNF257 or ZNF28, respectively. However, no ‘low’ SVA_A showed binding (Figure 2A), and both ‘high’ and ‘medium’ SVA_A copies significantly overlapped with the ZNF257- or ZNF28-binding peaks (Figure 2—figure supplement 1A). The frequency of overlap with the ZNF257/28 peaks and the enrichment of ZNF257/28 in SVA_A were positively correlated with DNA methylation (Figure 2B and C), and both ZNF257 and ZNF28 showed the highest enrichment of SVA_A when SVA family members were compared (Figure 2D).

Figure 2 with 2 supplements see all
Identification of Krüppel-associated box domain zinc finger proteins (KRAB-ZFPs) associated with DNA demethylation resistance in SINE-VNTR-Alus (SVAs).

(A) Scatter plots showing the fraction of low-methylated or highly methylated SVA_A copies which overlaps of KRAB-ZFP peaks. ZNF257 and ZNF28 peaks were more frequently overlapped with ‘high’ methylated SVA_A than ‘low’ methylated SVA_A. For this analysis, publicly available ChIP-exo data from 250 human KRAB-ZFPs in HEK293T cells were used. (B) Bar graphs showing the fraction of SVA_A copies with ZNF257 and ZNF28 peaks. SVA_A copies were classified by DNA methylation levels in 19 W human primordial germ cells (hPGCs) (N = 8, 26, 47 in low, medium, and high, respectively). p-Value was calculated by chi-square test. (C) Enrichment of ZNF257 and ZNF28 on SVA_A classified by DNA methylation levels in male hPGCs at 19 weeks of gestation. (D) Enrichment of ZNF257 and ZNF28 on each SVA subtype. (E) Position of ZNF257- and ZNF28-binding motifs along SVA consensus sequences. VNTR1 and VNTR2 are composed of multiple copy number of tandem repeats, and the copy number of these number tandem repeats (VNTRs) is highly variable among SVA copies. Both ZNF257- and ZNF28-binding motifs were found within VNTR1 of SVAs. (F) Sequence logo of ZNF257- and ZNF28-binding motifs. (G) Violin plots showing copy number of VNTR1 of each SVA subtype. (H) Violin plots showing VNTR1 copy number of SVA_A classified by its DNA methylation status in male hPGCs at 19 weeks of gestation. (I) Violin plots showing the number of ZNF257 and ZNF28 motifs in SVA_A classified by DNA methylation status in male hPGCs at 19 weeks of gestation. p-Value was calculated by Tukey’s test.

It is possible that the correlation between the DNA demethylation resistance of SVA_A and the binding potential of specific KRAB-ZNFs based on ChIP-exo mapping data in HEK293T cells could result from differences in read mappability. To determine the likelihood of this, we calculated the mappability of each transposon copy by virtually creating reads from the retroelements and mapping them onto the genome. Although highly methylated SVA_A copies showed greater mappability than those that were lowly methylated (Figure 2—figure supplement 1B), the correlation between SVA_A DNA methylation levels and enrichment for ZNF28/257 was observed even when only SVA_A copies with similar mappability (50–70%) were used for analysis (Figure 2—figure supplement 1C). Therefore, we concluded that the enrichment of ZNF28/257 in SVA_A in HEK293T cells is correlated with SVA_A DNA methylation levels in hPGCs.

The SVA element has a region containing variable-number tandem repeats (VNTRs) in the middle segment. SVA_A contains one type of VNTR (VNTR1), whereas the other SVA classes possess two types of VNTRs (VNTR1 and VNTR2) (Figure 2E). The ZNF257- and ZNF28-binding motifs, which were predicted by HOMER (Heinz et al., 2010; Figure 2F), are in VNTR1 (Figure 2E, Figure 2—figure supplement 1D). The number of ZNF257- and ZNF28-binding motifs within SVAs was the highest in SVA_A (Figure 2E) and was most strongly correlated with the copy number of VNTR1 in SVA_A out of all the SVA classes (Figure 2G). The VNTR1 copy number was also highly variable among SVA_A copies (Figure 2G), and DNA methylation of SVA_A was positively correlated with the VNTR1 copy number (Figure 2H, Figure 2—figure supplement 1E) and the number of ZNF257/28 motifs (Figure 2I, Figure 2—figure supplement 1F). We also confirmed that DNA methylation levels within the VNTR were correlated with ZNF257 or ZNF28 association (Figure 2—figure supplement 1G, H). These results indicate that a high number of ZNF257- and ZNF28-binding motifs within the VNTR increases the enrichment of KRAB-ZFPs. This might contribute to maintaining SVA_A DNA methylation during hPGC development. We confirmed the RNA expression of ZNF257 and ZNF28 in hPGCs by reanalysis of single-cell RNA-seq data from hPGCs and neighboring somatic cells (Guo et al., 2015; Figure 2—figure supplement 2A, B). However, there was no direct evidence for ZNF28/257 protein expression and its binding to SVAs in hPGCs, which warrants further studies.

The presence of the ZNF649-binding motif is correlated with demethylation resistance in L1s

We also analyzed the correlation between KRAB-ZFP-binding motifs and the DNA methylation status of L1s and LTR12s in hPGCs. Consistent with previous reports that ZNF649 and ZNF93 bind L1s (Cosby et al., 2019; Jacobs et al., 2014), ZNF649 and ZNF93 peaks were frequently found in L1PA2–6 and L1PA3–6, respectively (Figure 3A), and these two KRAB-ZFPs were enriched at the 5ʹ terminus of the L1 sequences (Figure 3B). The frequency of L1 copies overlapping with ZNF649 and ZNF93 peaks was correlated with the DNA methylation levels of L1s in hPGCs (Figure 3C, Figure 3—figure supplement 1A). Because read mappability in L1 (L1PA4) was similar across the different DNA methylation groups (Figure 3—figure supplement 1B), ZNF649 and ZNF93 are candidate factors for the DNA demethylation resistance of these L1s. As was the case for SVA_A, the presence of ZNF649- or ZNF93-binding motifs (Figure 3D) was also correlated with DNA methylation levels (Figure 3E).

Figure 3 with 2 supplements see all
Identification of Krüppel-associated box domain zinc finger proteins (KRAB-ZFPs) associated with DNA demethylation resistance in L1.

(A) Bar graphs showing the fraction of full-length L1 copies with ZNF93 and ZNF649 peaks. (B) Heatmaps showing enrichment of ZNF649 and ZNF93 along full-length L1 copies. ZNF649 binds 5’ regions of L1PA2–PA8, while ZNF93 binds the 5’ regions of L1PA3–PA8. (C) Bar graphs showing the fraction of L1 copies with ZNF93 and ZNF649 peaks. L1 copies were classified by their type and DNA methylation levels in male human primordial germ cells (hPGCs) at 19 weeks of gestation. (D) Sequence logo of ZNF93- and ZNF649-binding motifs. (E) Bar graphs showing the fraction of L1 copies with ZNF93- and ZNF649-binding motifs. The presence of ZNF93- and ZNF649-binding motifs was correlated with higher DNA methylation of L1 in male hPGCs at 19 weeks of gestation (L1PA2 and -PA6 were not significant for ZNF93). p-Value was calculated by Hypothesis Testing for the Difference in the Population Proportions using a function of prop.test by R. (F) Comparison of sequences of ZNF649-binding sites among L1 types. L1HS lost the ZNF649 motif by a base substitution. (G) Comparison of sequences at ZNF649-binding sites between low- and high-methylated L1. (H) Representative view of correlation between DNA methylation of L1PA4 in hPGCs and ZNF649-binding peak.

Reanalysis of single-cell RNA-seq data for hPGCs and neighboring somatic cells (Guo et al., 2015) showed that both ZNF649 and ZNF93 were expressed more in hPGCs than in neighboring somatic cells (Figure 2—figure supplement 2A, B). Because the correlation between the presence of binding motifs and DNA methylation levels was stronger in ZNF649 than in ZNF93 (Figure 3E), we investigated ZNF649 in more detail. The ZNF649-binding motif was located at the 5ʹ UTR of L1s (Figure 3F), consistent with the enrichment of ZNF649 in the 5′ UTR (Figure 3B). The enrichment of ZNF649 in L1s was decreased in L1PA2 and abrogated in L1HS (Figure 3B). Along with the decreased ZNF649 enrichment, a base substitution at the fifth position of the ZNF649-binding site was observed in the consensus sequences of L1HS (Figure 3F). Because the fifth position of the ZNF649-binding site (T) is conserved in highly methylated L1 copies (Figure 3G), a T in this position may be required for ZNF649 to bind to L1. Although highly methylated L1 copies had two mismatches within the ZNF649-binding motif, one at the third position (T→G) and one at the sixth position (A→T) (Figure 3G), a minor fraction of the ZNF649-binding motif had the same base composition at these sites (Figure 3D). Thus, these two mismatches may not abrogate ZNF649 binding. We also confirmed high DNA methylation in the ZNF649-binding motifs at individual loci (Figure 3H).

The presence of the ZNF850-binding motif is correlated with demethylation resistance in LTR12C

For the LTR12C family, we found that ZNF850 more frequently overlapped with highly methylated LTR12C/D/E copies than lowly methylated ones when we analyzed the binding peak data for 250 KRAB-ZFPs (Helleboid et al., 2019; Imbeault et al., 2017; Figure 3—figure supplement 2A, B). We focused on LTR12C because it had the highest analyzable copy number (LTR12C: 2054; LTR12D: 130; LTR12E: 46 copies). The ZNF850-binding motif was more frequently found in highly methylated LTR12C copies than in lowly methylated copies (Figure 3—figure supplement 2C). Two high-confidence binding motifs (q-value < 0.01) were identified at the 5′ portion of LTR12C consensus sequences (Figure 3—figure supplement 2D), which was consistent with ZNF850 enrichment in the 5′ portion of LTR12C (Figure 3—figure supplement 2E). Lowly methylated LTR12C copies contained mismatches more frequently at the eighth and tenth positions of the first and second predicted binding sites along LTR12C, respectively (Figure 3—figure supplement 2F). An example of highly methylated LTR12C loci with a ZNF850 peak is shown in Figure 3—figure supplement 2G. These data suggest that KRAB-ZNFs prevent DNA demethylation during male germ-cell development.

The mode of DNA methylation acquisition during spermatogenesis varies depending on retroelement type

To investigate whether lowly methylated retroelements in hPGCs acquire DNA methylation during spermatogenesis, we analyzed publicly available human sperm WGBS data from two donors (Hammoud et al., 2014). The two donors were of similar age (donor 1 – 32 and donor 2 – 37), and both were white Caucasians. The dynamics of DNA methylation in retroelements during spermatogenesis vary depending on retroelement type and individual characteristics. Most L1 copies acquired DNA methylation during spermatogenesis in both individuals, whereas LTR12C copies maintained their DNA methylation status in hPGCs during spermatogenesis (Figure 4A). A substantial difference between individuals was observed in the SVAs. The majority of SVA copies acquired DNA methylation during spermatogenesis in sperm donor 1, but not in donor 2 (Figure 4A). To evaluate these trends more efficiently, we classified retroelement copies based on DNA methylation levels in sperm (common high: > 60% in both donors; high and low: > 60% in donor 1 and < 20% in donor 2; common low: < 20% in both donors). The majority of lowly methylated L1 copies in hPGCs were highly methylated in sperm cells from both donors (Figure 4B). In contrast, most LTR12C/D copies maintained their PGC DNA methylation status during spermatogenesis (Figure 4C). Among the SVA types, SVA_A showed high levels of DNA methylation in both sperm donors, whereas other SVA types showed variable DNA methylation levels when both sperm donors were compared (Figure 4D), especially in SVA copies that had low DNA methylation levels in hPGCs (Figure 4E).

DNA methylation dynamics of retroelements during human spermatogenesis.

(A) Scatter plots showing DNA methylation levels of each retroelement copy in male human primordial germ cells (hPGCs) at 19 weeks of gestation and sperm. Whole-genome bisulfite sequencing (WGBS) data from two sperm donors (Hammoud et al., 2014) were used for this analysis. Donor 1 and donor 2 were colored by orange and cyan, respectively. (B–D) Bar graphs showing the fraction of groups determined by DNA methylation patterns in two sperm donors in L1 (B), LTR12 (C), and SINE-VNTR-Alus (SVA) (D). ‘Other’ indicates groups except for common high, common low, and high and low, such as low methylated in donor 1 and mediumly methylated in donor 2. Bar graphs were also separated by DNA methylation levels (high or low) in male hPGCs at 19 weeks of gestation. (E) Violin plots showing DNA methylation levels of SVA copies in male hPGCs at 19 weeks of gestation, sperm donor 1, and sperm donor 2. The violin plots were also separated by DNA methylation levels of SVA copies in male hPGCs at 19 weeks of gestation. Although hypomethylated SVA copies in male hPGCs at 19 weeks of gestation acquired DNA methylation during spermatogenesis, the degree of DNA methylation increase was significantly different between sperm donors. p-Value was calculated by Dunnett’s test.

The degree of DNA methylation acquisition during spermatogenesis varies among SVA copies

Although the DNA methylation status of SVAs was highly variable between the sperm donors, a subset of SVA copies acquired DNA methylation or maintained a low methylation state during spermatogenesis in both sperm donors (Figure 4D). It is possible to get insight for mechanisms of de novo DNA methylation of SVAs during spermatogenesis by comparing SVAs that acquired DNA methylation (‘low’ in hPGCs and ‘common high’ in sperm) to SVAs that maintained hypomethylation (‘low’ in hPGCs and ‘common low’ in sperm) in both sperm donors. Phylogenetic analysis of ‘common low’ and ‘common high’ SVA copies (SVA_B and _D) showed that these two classes were not genetically separated (Figure 5A), indicating that the acquisition of DNA methylation in SVAs during spermatogenesis is not genetically determined.

Transcription-associated regulation of DNA methylation of SINE-VNTR-Alus (SVA) during spermatogenesis.

(A) Phylogenetic analysis of SVA_B (left) and SVA_D (right) copies low methylated in male human primordial germ cells (hPGCs) at 19 weeks of gestation. SVA copies highly methylated by both sperm donors were colored by red, while those hypomethylated by both sperm donors were colored by blue. (B) Bar graphs showing the fraction of SVA_B–F copies inserted in a gene body. SVA copies were classified by DNA methylation patterns in two sperm donors. Only low-methylated SVA copies in male hPGCs at 19 weeks of gestation were used for this analysis. The number in parentheses represents analyzed copy number. (C) Violin plots showing the expression of genes in adult spermatogonial stem cells 2 (Sohni et al., 2019). Genes were classified according to the DNA methylation status of SVAs inserted in them in the antisense direction. p-Value was calculated by Tukey’s test. (D) Enrichment of RNA-seq reads from undifferentiated spermatogonia (Tan et al., 2020) around non-genic SVAs. Only low-methylated SVA copies male hPGCs at 19 weeks of gestation were used for the analysis, and SVA copies were classified by DNA methylation patterns in two sperm donors (common low, high, and low and common high).

The presence of transcription-directed retroelement silencing mechanisms, such as the PIWI/piRNA pathway (Watanabe et al., 2018), prompted us to investigate the correlation between the genomic distribution of SVA copies and DNA methylation. Approximately half of the SVA_B–F copies were inserted into the gene body, and most of them were in the antisense direction (Figure 5B). ‘Common high’ SVA_B–F copies were enriched in the gene body in the antisense direction, while ‘common low’ SVA_B–F copies were depleted from the gene body (Figure 5B). Reanalysis of publicly available single-cell RNA-seq data in human testes (Sohni et al., 2019) revealed that genes with ‘common high’ SVA_B–F copies in the antisense orientation showed greater expression in spermatogonial stem cells relative to genes with ‘common low’ (Figure 5C). Therefore, SVAs located in actively transcribed regions in the antisense orientation are efficiently subjected to de novo DNA methylation during spermatogenesis. The expression of genes with ‘high and low’ SVA_B–F copies in the antisense direction was higher in spermatogonial stem cells than the expression of other randomly extracted genes and genes with ‘common low’ SVA_B–F copies. However, the expression of these genes was lower than the expression of genes with ‘common high’ SVA_B–F copies (Figure 5C). Approximately half of the ‘high and low’ SVA_B–F copies were located in non-genic regions, but RNA-seq reads from previously reported undifferentiated spermatogonia (Tan et al., 2020) mapped more frequently around the non-genic ‘high and low’ SVA_B–F copies than the ‘common low’ B–F copies (Figure 5D). Therefore, non-genic ‘high and low’ SVA_B–F copies are frequently inserted in transcribed regions during spermatogenesis. These results implicate the possibility that SVAs acquire DNA methylation during DNA methylation via transcription-directed machinery, and that the effectiveness of de novo DNA methylation varies among individuals.

SVAs are a potential source of inter-individual epigenetic variation in sperm

Inter-individual variation in DNA methylation in SVAs was also observed when another set of publicly available sperm WGBS data from three Japanese donors was analyzed (Okae et al., 2014; Figure 6A). For additional validation, we performed amplicon sequencing (amplicon-seq) of bisulfite PCR products for SVAs on sperm from five Japanese donors (Figure 6B). Our amplicon-seq yielded approximately 1.7–2.2 M read pairs and measured the DNA methylation level of over 90% of the full-length SVA_B–F copies (minimum read depth of CpG ≥ 5, analyzed CpG number ≥ 10) (Figure 6C). Again, ‘high and low’ SVA_B–F copies showed variations in DNA methylation among donors (Figure 6D). Thus, inter-individual variation in SVA methylation in sperm is a common phenomenon and is not ethnically specific.

Figure 6 with 1 supplement see all
SINE-VNTR-Alus (SVAs) constitute a major source of inter-individual epigenetic variations in sperm.

(A) Violin plots showing DNA methylation of SVA copies in previously reported three sperm donors (#1–#3) (Okae et al., 2014). Only low-methylated SVA copies in male human primordial germ cells (hPGCs) at 19 weeks of gestation were used for the analysis. SVA copies were classified by DNA methylation levels of two sperm donors from Hammoud et al., 2014. Donor #1 showed significantly higher DNA methylation levels in ‘high and low’ SVA copies than other sperm donors. p-Value was calculated by Tukey’s test. (B) Scheme of amplicon sequencing (amplicon-seq) for analyzing SVA methylation. (C) Bar plots showing the fraction of analyzed full-length SVA copies by amplicon-seq. (D) Violin plots showing DNA methylation levels of SVA copies in five sperm donors from amplicon-seq. Only low-methylated SVA copies in male hPGCs at 19 weeks of gestation were used for the analysis. SVA copies were classified by DNA methylation levels of two sperm donors from Hammoud et al. Donor #5 showed significantly higher DNA methylation levels in ‘high and low’ SVA copies than other sperm donors. p-Value was calculated by Tukey’s test. (E) Scatter plot showing the DNA methylation between sperm donor 1 and sperm donor 2 from Hammoud et al. DNA methylation levels between these two donors were highly correlated. (F) Heatmap showing DNA methylation levels, genomic distribution, and overlap with SVAs of differentially methylated regions (DMRs). (G) Representative view of DMRs overlapping with SVA. Black and green boxes represent SVA and DMR, respectively.

To estimate the impact of SVAs on inter-individual epigenetic variations in sperm, we identified differentially methylated regions (DMRs) in two sperm donors, as shown in Figures 4 and 5; Hammoud et al., 2014. Although the DNA methylation profiles between the two donors were highly correlated (Figure 6E), 2008 regions were identified as DMRs (donor 1 < donor 2: 332, donor 2 < donor 1: 1676). Of the 1676 donor 1-specific methylated DMRs, 772 (46.1%) overlapped with SVAs (Figure 6F). We also observed differential DNA methylation among individuals in SVA-associated DMRs in our amplicon-seq data (Figure 6G). Therefore, SVAs significantly contribute to inter-individual variations in the sperm epigenome. In contrast to the inter-individual epigenetic variation of SVAs in sperm, a reanalysis of WGBS data of adult skeletal muscle from 15 individuals and of helper CD4-positive T cells from 18 individuals, which was deposited in the International Human Epigenome Consortium (IHEC) portal, showed high DNA methylation of SVAs in all individuals (Bujold et al., 2016; Figure 6—figure supplement 1A, B). Thus, inter-individual variations in DNA methylation of SVAs in sperm are essentially canceled during development.

Finally, to investigate whether inter-individual DNA methylation variations are associated with physiological or disease conditions, we reanalyzed publicly available WGBS data from five healthy donors and six oligozoospermic patients (European Nucleotide Archive [ENA] under the accession number PRJEB34432) (Leitão et al., 2020). The disease condition was not associated with inter-individual DNA methylation variations in SVAs, because both healthy donors and oligozoospermic patients showed inter-individual variations of DNA methylation in ‘high and low’ SVA_B–F copies (Figure 6—figure supplement 1C). On the other hand, a comparison of various physiological conditions between highly methylated individuals and lowly methylated ones (median methylation levels of ‘high and low’ > 50% vs. < 50%) revealed that blood testosterone levels were significantly higher in lowly methylated individuals than in highly methylated ones (Figure 6—figure supplement 1D). However, prolactin, follicle stimulating hormone, luteinizing hormone (LH), sex hormone-binding globulin blood levels, and age were not significantly different between the two groups (Figure 6—figure supplement 1D). Although further validation of this correlation is required, DNA methylation of SVAs in sperm may be associated with physiological conditions.

Discussion

In this study, we showed that the binding potential of KRAB-ZFPs correlates with retroelement DNA demethylation resistance in hPGCs. Furthermore, we found that de novo DNA methylation patterns in spermatogenesis vary among the L1, LTR, and SVA retroelements. In addition, we ascertained that the SVAs located in transcription-active regions in the antisense orientation are prone to methylation during spermatogenesis, which implies that the transcription-directed DNA methylation machinery might contribute to de novo DNA methylation of SVAs in male germ cells. Notably, the extent of de novo DNA methylation of SVAs in male germ cells is variable among human individuals, with SVAs being a major source of epigenetic variation in sperms.

We showed that DNA demethylation resistance in hPGCs frequently occurred in moderately young retroelements such as L1PA, SVA_A, and LTR12, but not in currently active retroelements. Because we targeted full-length transposons, our analysis was biased toward relatively young transposons. Thus, it is possible that some fragmented older transposons may also be resistant to DNA demethylation in hPGCs. A subset of LTR transposons, including LTR12, function as enhancers (Deniz et al., 2020). It was recently reported that LTR5s, which are Hominidae-specific LTR-type transposons and hypomethylated in hPGCs (DNA methylation levels < 10%), can function as enhancers to promote hPGC differentiation (Xiang et al., 2022). Therefore, in the case of LTR12C, maintaining DNA methylation might be beneficial for hPGC development because it suppresses inappropriate activation of transposon-embedded enhancer function.

In addition, KRAB-ZFP binding potentially contributed to the DNA demethylation resistance of L1s and SVAs. ZNF257/28, ZNF649/ZNF93, and ZNF850 were associated with the DNA demethylation resistance of SVAs, L1s, and LTR12Cs, respectively (Figure 7). In hPGCs, multiple KRAB-ZNPs were correlated with DNA demethylation resistance in the same retroelements, which may contribute to more robust or cooperative retroelement suppression. Although ZNF91 reportedly binds to the VNTR in SVAs and silences SVA expression in embryonic stem cells (ESCs) (Haring et al., 2021; Jacobs et al., 2014), the DNA demethylation resistance of SVAs did not correlate with ZNF91 binding, indicating that a different KRAB-ZFP set is used to suppress SVAs in human PGCs and ESCs. Both ZNF257 and ZNF28 bound to VNTR1 (Figure 2E), and high copy numbers of VNTR1 were correlated with high ZNF257 and ZNF28 enrichment and DNA methylation (Figure 2H1). The reduction in VNTR1 copy number after SVA_B emergence (Figure 2G) may have been necessary for SVAs to escape silencing mechanisms in hPGCs. The human genome encodes at least 350 KRAB-ZFPs, and not all KRAB-ZFPs were included in the ChIP-seq dataset used in this study (100 copies remained unmapped). Thus, the involvement of other KRAB-ZFPs in DNA demethylation resistance of retroelements in hPGCs is possible. Although we observed a strong correlation between KRAB-ZFPs and DNA demethylation resistance, direct evidence for this correlation remains elusive because of the limited availability of human fetal gonads and of high-specificity antibodies for KRAB-ZFPs. Because they can function as in vitro derivation systems, PGC-like cells (PGCLCs) may be a promising model for investigating the biology of PGCs. Although successful establishment of human PGCLCs has been reported (Sasaki et al., 2015), sufficient DNA demethylation has not been observed in human PGCLCs (von Meyenn et al., 2016). Thus, the currently available human PGCLCs are not suitable models for investigating the mechanisms of DNA demethylation resistance. Optimizing the derivation conditions for human PGCLCs will aid in our understanding of retroelement silencing in PGCs.

Summary of this study.

Our data demonstrated an association between KRAB-ZFP-binding motifs and the DNA demethylation resistance of L1, SVA, and LTR12. ZNF649, ZNF257/28, and ZNF850 were associated with DNA demethylation resistance of L1, SVA, and LTR12C, respectively. The dynamics of DNA methylation during spermatogenesis are largely different among retroelement types. The majority of L1 copies acquired DNA methylation during spermatogenesis, whereas the DNA methylation status of LTR12 in human primordial germ cells (hPGCs) tended to be maintained during spermatogenesis. The mode of DNA methylation changes in SINE-VNTR-Alus (SVAs) during spermatogenesis largely differs between copies and individuals. SVA copies located in highly transcriptionally active regions acquire DNA methylation during spermatogenesis, while those located in transcriptionally inactive regions maintain a hypomethylated state during spermatogenesis. In contrast, the degree of DNA methylation in sperm in SVA copies located in low transcriptionally active regions was highly variable among the individuals. These results suggest that SVAs may be methylated by transcription-directed DNA methylation mechanisms during spermatogenesis, and their activity varies among individuals.

Additionally, we showed that the mode of DNA methylation acquisition during spermatogenesis was very different among retroelement types. The majority of L1 copies acquired DNA methylation during spermatogenesis, whereas LTR12 maintained its DNA methylation status in hPGCs during spermatogenesis (Figure 7). L1HS, which both ZNF93 and ZNF649 were unable to bind, also acquired DNA methylation during spermatogenesis (Figure 4B), suggesting the involvement of other factors in the de novo DNA methylation of L1 during spermatogenesis. The PIWI-piRNA pathway is responsible for the DNA methylation of L1 transposons in mouse male germ cells (Aravin et al., 2007; Carmell et al., 2007; Kojima-Kita et al., 2016; Manakov et al., 2015; Shoji et al., 2009). The PIWI-piRNA pathway may also be functional in humans, because mutations in genes associated with the PIWI-piRNA pathway are linked to human male infertility (Arafat et al., 2017; Gu et al., 2010). Moreover, the majority of putative piRNAs that mapped to transposons at gestational week 20 are derived from L1 (Reznik et al., 2019). Therefore, the PIWI-piRNA pathway is a candidate pathway for L1 silencing in human male germ cells.

Our data showed that the acquisition of DNA methylation of SVAs during spermatogenesis correlated with the inserted regions and not with the nucleotide sequence. SVAs inserted in transcriptionally active regions in the antisense direction are efficiently targeted for de novo DNA methylation during spermatogenesis. In mouse spermatogenesis, MIWI2 binds piRNAs and is recruited to the nascent transcribed regions that are complementary to piRNAs (Watanabe et al., 2018). Subsequently, MIWI2-interacting protein SPOCD1, which forms a complex with DNMT3A and DNMT3L, and potentially with a rodent-specific DNA methyltransferase DNMT3C (Barau et al., 2016), induces DNA methylation on transposons (Zoch et al., 2020). Therefore, one possible mechanism for the de novo DNA methylation of SVAs during spermatogenesis is that the MIWI2/SVA-derived piRNA complex targets nascent transcripts with antisense SVAs and induces DNA methylation. There are also other transcription-directed repetitive element silencing mechanisms, such as those involving the HUSH complex, which repress L1s and SVAs (Fukuda et al., 2018; Liu et al., 2018; Robbez-Masson et al., 2018). The HUSH complex targets young full-length L1s located within the introns of actively transcribed genes (Fukuda and Shinkai, 2020; Liu et al., 2018). In addition to the HUSH complex, efficient pericentromeric heterochromatin formation requires the transcription of pericentromeric satellite repeats, which stabilize SUV39H pericentromeric localization (Johnson et al., 2017; Shirai et al., 2017; Velazquez Camacho et al., 2017). Because SUV39H is also associated with retroelement silencing (Bulut-Karslioglu et al., 2014), both the HUSH complex and SUV39H are candidate factors associated with the transcription-directed DNA methylation of SVAs in human male germ cells. In eukaryotes, gene bodies are the most conserved targets of DNA methylation. Gene body DNA methylation levels are often correlated with transcriptional levels (Teissandier and Bourc’his, 2017). This is because of the interaction between the elongating RNA polymerase II and SETD2, which results in H3K36me3. H3K36me3 participates in the de novo methylation of DNA by recruiting DNMT3 enzymes via their chromatin reading PWWP domains (Baubec et al., 2015; Shirane et al., 2020). Furthermore, antisense RNAs embedded within protein-coding genes are selectively silenced by H3K36 methyltransferase SET2 in Saccharomyces cerevisiae (Venkatesh et al., 2016). Thus, the machinery for gene body DNA methylation regulated by SETD2 is also a candidate for the de novo DNA methylation of SVAs during spermatogenesis.

The mechanism underlying inter-individual epigenetic variations in SVA in human sperm is unknown. In addition to genetic differences among individuals, both intrinsic and extrinsic environmental differences may contribute to inter-individual variations in SVAs. Our data indicate that in sperm, the degree of DNA methylation of SVAs located in genomic regions with low transcriptional activity varies among individuals. Thus, the effectiveness of transcription-directed de novo DNA methylation in male human germ cells may vary among individuals. Previous studies have shown that hypermethylation of the PIWIL2 and TDRD1 promoter regions, which are involved in the PIWI-piRNA pathway, is associated with abnormal DNA methylation and male infertility in humans (Heyn et al., 2012). Therefore, the effectiveness of the PIWI-piRNA pathway may vary among individuals and contribute to the epigenetic variation of SVAs in male germ cells. SVAs function as enhancers (Gianfrancesco et al., 2017), alter the chromatin state near the insertion site (Fukuda et al., 2017), and are associated with Fukuyama-type congenital muscular dystrophy and Lynch syndrome (Ostertag et al., 2003; Payer and Burns, 2019). Therefore, differences in SVA regulation among individuals may induce changes in gene regulation in male germ cells, alter the risk of genome instability, and affect the incidence of disease among individuals.

Materials and methods

Semen collection

Request a detailed protocol

Ejaculates were provided by patients who visited the Reproduction Center of the Ichikawa General Hospital, Tokyo Dental College. All study participants were briefed about the aims of the study and the parameters to be measured, and consent was obtained. The study was approved by the ethics committees of RIKEN, Tokyo University, and Ichikawa General Hospital. Sperm concentration and motility were measured using a computer-assisted image analyzer (C-Men, Compix, Cranberry Township, PA). Human semen was diluted twice with saline, layered on 5.0 mL of 20 mM HEPES buffered 90% isotonic Percoll (Cytiba, Uppsala, Sweden), and centrifuged at 400 × g for 22 min. The sperm in the sediment was recovered to yield 0.2 mL, and then introduced to the bottom of 2.0 mL of Hanks’ solution to facilitate swim-up. The motile sperm in the upper layer were recovered.

Preparation of SVA amplicon-seq

Request a detailed protocol

Genomic DNA was subjected to bisulfite-mediated C to U conversion using the MethylCode Bisulfite Conversion Kit (ThermoFisher Scientific), and then used as a template for PCR for 35 cycles with EpiTaq (Takara) using the following primers: SVA_1_Fw TTATTGTAATTTTTTTGTTTGATTTTTTTGTTTTAG. SVA_1_Rv AAAAAAACTCCTCACATCCCAAAC SVA_2_Fw TTAATGTTGTTTAGGTTGGAGTGTAGTG SVA_2_Rv CAAAAAAACTCCTCACTTCCCAATA. SVA_3_Fw TTTGGGAGGTGTATTTAATAGTTTATTGAGAA SVA_3_Rv TAAACAAAAATCTCTAATTTTCCTAAACAAAAAACC. The PCR products from three sets of primers were combined, purified using a MinElute PCR Purification Kit (QIAGEN), and fragmented using Picoruptor (Diagenode) for 10 cycles of 30 s on and 30 s off. Then, the amplicon-seq library was constructed using KAPA LTP Library Preparation Kits (KAPA BIOSYSTEMS) and SeqCap Adapter Kit A (Roche). The amplicon-seq libraries were sequenced on a HiSeq X platform (Illumina).

WGBS and amplicon-seq analysis

Request a detailed protocol

We used the hg19 version of the human genome for NGS analysis because the predicted KRAB-ZFP peaks were derived from this version. Using the newest version of the human genome (GRCh38) did not significantly affect the conclusions. The following publicly available WGBS data were used in this study: hPGCs (SRP050499), sperm (SRP028572, ERP117337, JGAS00000000006), and adult tissues (IHEC data portal). For the IHEC data, we used processed data for our analysis.

  • Quality control, read mapping, and DNA methylation calculation

Low-quality bases and adaptor sequences were trimmed using Trim Galore version 0.3.7 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). For WGBS data from hPGCs, the first nine bases were further trimmed. The trimmed reads were aligned to the hg19 genome using Bismark v0.14.1, with paired-end and non-directional mapping parameters (--non_directional) (Krueger and Andrews, 2011). The unmapped reads after paired-end mapping were re-aligned to the same reference in single-end mode. We validated that this mapping mode only reported uniquely mapped reads. The methylation level of each CpG site was calculated as follows: (number of methylated reads/number of total reads). Only CpG sites with at least five reads were used for all analyses. Only nearly full-length retroelements, whose length is 90% or more of the length of the consensus sequence of each subtype, were used for DNA methylation analysis of retroelements. We also included solo-LTR transposons in the DNA methylation analysis if they also possessed more than 90% of the consensus LTR sequence. Retroelement information was obtained from the UCSC Genome Browser (http://genome.ucsc.edu/). For the DNA methylation analysis of retroelements, we used retroelements containing at least 10 CpG sites with a read depth of at least five reads. The methylation level of each retroelement copy was calculated by averaging the methylation levels of CpG sites within the copy.

Classification of retroelement copy according to DNA methylation levels.

Retroelement copies were classified according to their DNA methylation levels as follows: low < 20%, 20% ≦ medium < 60%, high 60%.

  • Association of KRAB-ZFP peaks, binding motifs, and retroelements.

We obtained the peak regions of 250 KRAB-ZFP, which were previously reported (Helleboid et al., 2019; Imbeault et al., 2017), from the gene expression omnibus GSE78099 and GSE120539. Overlap of the KRAB-ZFP peak and retroelement copy was investigated using bedtools v2.15.0 (Quinlan and Hall, 2010). The binding motif of each KRAB-ZFP was predicted by the findMotifsGenome.pl program in Homer v4.8.3 (Heinz et al., 2010). The KRAB-ZFP-binding motifs along retroelement copies were searched using FIMO (Grant et al., 2011). We used predicted motif sites with a q-value of 0.00005 or less for ZNF257/ZNF28/ZNF850 and with a q-value of 0.05 or less for ZNF93 and ZNF649 in this study.

  • DMR identification

DMR candidates were identified using the ‘Commet’ command in BisulFighter (Saito et al., 2014). To enhance the confidence of DMR call, we calculated the average methylation levels of the candidates using CpG sites with ≥5 reads in both sperm donors, and among the candidates, those containing ≥10 successive analyzable CpG sites and showing a ≥40% methylation difference were determined as DMRs.

Phylogenetic analysis of retroelement copies

Request a detailed protocol

The evolutionary history was inferred using the maximum likelihood method based on the Tamura-Nei model (Tamura and Nei, 1993). The initial tree(s) for the heuristic search were obtained by applying the neighbor-joining method to a matrix of pairwise distances estimated using the maximum composite likelihood approach. The tree was drawn to scale, with branch lengths measured as the number of substitutions per site. There were 10,153 positions in the final dataset. Evolutionary analyses were conducted using MEGA6 (Tamura et al., 2013).

Calculation of read mappability of each retroelement copy

Request a detailed protocol

We generated 100 bp reads from each position along the retroelement copy and aligned the simulated reads to the human genome using Bowtie with –m 1 or Bismark. Then, the mappability of each copy was calculated by dividing the number of properly mapped reads by the total number of reads derived from each copy.

RNA-seq analysis

Request a detailed protocol

We reanalyzed previously reported single-cell RNA-seq data from the testes (Sohni et al., 2019), hPGCs, and somatic cells next to hPGCs (Guo et al., 2015). Read count data of genes and cell type annotation of each cell were obtained from the Gene Expression Omnibus under accession numbers GSE124263 and GSE63818. Reads per million mapped reads (RPM) for the genes were calculated for each cell. We used the average RPM of spermatogonial stem cells 2 (Figure 5C). We also reanalyzed previously reported RNA-seq data from undifferentiated spermatogonia (Tan et al., 2020), which was deposited in the Gene Expression Omnibus under accession number GSE144085. Low-quality bases and adaptor sequences were trimmed using Trim Galore version 0.3.7. Then, trimmed reads were aligned to the hg19 genome using Bowtie v0.12.7 with -m 1 to remove multiple mapped reads. Enrichment of RNA-seq reads around SVAs was visualized using ngsplot (Shen et al., 2014).

ChIP-seq analysis

Request a detailed protocol

We reanalyzed previously reported KRAB-ZFP ChIP-exo data (Helleboid et al., 2019; Imbeault et al., 2017), which were deposited in the Sequence Read Archive SRP070561 and SRP162756. Low-quality bases and adaptor sequences were trimmed using the Trim Galore version 0.3.7. Then, trimmed reads were aligned to the hg19 genome using Bowtie v0.12.7 with -m 1 to avoid multiple mapped reads. Enrichment of ChIP-exo reads around retroelements was visualized using ngsplot (Shen et al., 2014).

Visualization of NGS data

Request a detailed protocol

The Integrative Genomics Viewer (Robinson et al., 2011) was used to visualize the NGS data. Enrichment of RNA-seq reads and KRAB-ZFPs was visualized using ngsplot (Shen et al., 2014). Scatter plot and violin plot analyses were performed using the ggplot2 package in R.

Data access

Request a detailed protocol

All reads from amplicon-seq in this study have been submitted to the Gene Expression Omnibus under accession number GSE174562.

Data availability

All reads from amplicon-seq in this study have been submitted to the Gene Expression Omnibus under accession number GSE174562.

The following data sets were generated
    1. Fukuda K
    2. Shinkai Y
    (2021) NCBI Gene Expression Omnibus
    ID GSE174562. Amplicon-seq of SVA methylation in human sperm.
The following previously published data sets were used
    1. Imbeault M
    2. Helleboid PY
    3. Trono D
    (2017) NCBI Gene Expression Omnibus
    ID GSE78099. ChIP-exo of human KRAB-ZNFs transduced in HEK 293T cells and KAP1 in hES H1 cells.
    1. Guo F
    2. Guo H
    3. Li L
    4. Tang F
    (2015) NCBI Gene Expression Omnibus
    ID GSE63818. The Transcriptome and DNA Methylome Landscapes of Human Primordial Germ Cells.
    1. Low DH
    2. Hammoud SS
    (2014) NCBI Gene Expression Omnibus
    ID GSE49624. Chromatin and Transcription Transitions of Mammalian Adult Germline Stem Cells and Spermatogenesis.
    1. Okae H
    2. Chiba H
    3. Hiura H
    4. Hamada H
    5. Sato A
    6. Utsunomiya T
    7. Kikuchi H
    8. Yoshida H
    9. Tanaka A
    10. Suyama M
    11. Arima T
    (2014) the Japanese Genotype-phenotype Archive
    ID JGAS00000000006. Genome-Wide Analysis of DNA Methylation Dynamics during Early Human Development.
    1. Sohni A
    2. Tan K
    3. Song H
    4. Burow D
    5. Wilkinson MF
    (2018) NCBI Gene Expression Omnibus
    ID GSE124263. Neonatal and adult human testis defined at the single-cell level.

References

    1. Lander ES
    2. Linton LM
    3. Birren B
    4. Nusbaum C
    5. Zody MC
    6. Baldwin J
    7. Devon K
    8. Dewar K
    9. Doyle M
    10. FitzHugh W
    11. Funke R
    12. Gage D
    13. Harris K
    14. Heaford A
    15. Howland J
    16. Kann L
    17. Lehoczky J
    18. LeVine R
    19. McEwan P
    20. McKernan K
    21. Meldrim J
    22. Mesirov JP
    23. Miranda C
    24. Morris W
    25. Naylor J
    26. Raymond C
    27. Rosetti M
    28. Santos R
    29. Sheridan A
    30. Sougnez C
    31. Stange-Thomann Y
    32. Stojanovic N
    33. Subramanian A
    34. Wyman D
    35. Rogers J
    36. Sulston J
    37. Ainscough R
    38. Beck S
    39. Bentley D
    40. Burton J
    41. Clee C
    42. Carter N
    43. Coulson A
    44. Deadman R
    45. Deloukas P
    46. Dunham A
    47. Dunham I
    48. Durbin R
    49. French L
    50. Grafham D
    51. Gregory S
    52. Hubbard T
    53. Humphray S
    54. Hunt A
    55. Jones M
    56. Lloyd C
    57. McMurray A
    58. Matthews L
    59. Mercer S
    60. Milne S
    61. Mullikin JC
    62. Mungall A
    63. Plumb R
    64. Ross M
    65. Shownkeen R
    66. Sims S
    67. Waterston RH
    68. Wilson RK
    69. Hillier LW
    70. McPherson JD
    71. Marra MA
    72. Mardis ER
    73. Fulton LA
    74. Chinwalla AT
    75. Pepin KH
    76. Gish WR
    77. Chissoe SL
    78. Wendl MC
    79. Delehaunty KD
    80. Miner TL
    81. Delehaunty A
    82. Kramer JB
    83. Cook LL
    84. Fulton RS
    85. Johnson DL
    86. Minx PJ
    87. Clifton SW
    88. Hawkins T
    89. Branscomb E
    90. Predki P
    91. Richardson P
    92. Wenning S
    93. Slezak T
    94. Doggett N
    95. Cheng JF
    96. Olsen A
    97. Lucas S
    98. Elkin C
    99. Uberbacher E
    100. Frazier M
    101. Gibbs RA
    102. Muzny DM
    103. Scherer SE
    104. Bouck JB
    105. Sodergren EJ
    106. Worley KC
    107. Rives CM
    108. Gorrell JH
    109. Metzker ML
    110. Naylor SL
    111. Kucherlapati RS
    112. Nelson DL
    113. Weinstock GM
    114. Sakaki Y
    115. Fujiyama A
    116. Hattori M
    117. Yada T
    118. Toyoda A
    119. Itoh T
    120. Kawagoe C
    121. Watanabe H
    122. Totoki Y
    123. Taylor T
    124. Weissenbach J
    125. Heilig R
    126. Saurin W
    127. Artiguenave F
    128. Brottier P
    129. Bruls T
    130. Pelletier E
    131. Robert C
    132. Wincker P
    133. Smith DR
    134. Doucette-Stamm L
    135. Rubenfield M
    136. Weinstock K
    137. Lee HM
    138. Dubois J
    139. Rosenthal A
    140. Platzer M
    141. Nyakatura G
    142. Taudien S
    143. Rump A
    144. Yang H
    145. Yu J
    146. Wang J
    147. Huang G
    148. Gu J
    149. Hood L
    150. Rowen L
    151. Madan A
    152. Qin S
    153. Davis RW
    154. Federspiel NA
    155. Abola AP
    156. Proctor MJ
    157. Myers RM
    158. Schmutz J
    159. Dickson M
    160. Grimwood J
    161. Cox DR
    162. Olson MV
    163. Kaul R
    164. Raymond C
    165. Shimizu N
    166. Kawasaki K
    167. Minoshima S
    168. Evans GA
    169. Athanasiou M
    170. Schultz R
    171. Roe BA
    172. Chen F
    173. Pan H
    174. Ramser J
    175. Lehrach H
    176. Reinhardt R
    177. McCombie WR
    178. de la Bastide M
    179. Dedhia N
    180. Blöcker H
    181. Hornischer K
    182. Nordsiek G
    183. Agarwala R
    184. Aravind L
    185. Bailey JA
    186. Bateman A
    187. Batzoglou S
    188. Birney E
    189. Bork P
    190. Brown DG
    191. Burge CB
    192. Cerutti L
    193. Chen HC
    194. Church D
    195. Clamp M
    196. Copley RR
    197. Doerks T
    198. Eddy SR
    199. Eichler EE
    200. Furey TS
    201. Galagan J
    202. Gilbert JG
    203. Harmon C
    204. Hayashizaki Y
    205. Haussler D
    206. Hermjakob H
    207. Hokamp K
    208. Jang W
    209. Johnson LS
    210. Jones TA
    211. Kasif S
    212. Kaspryzk A
    213. Kennedy S
    214. Kent WJ
    215. Kitts P
    216. Koonin EV
    217. Korf I
    218. Kulp D
    219. Lancet D
    220. Lowe TM
    221. McLysaght A
    222. Mikkelsen T
    223. Moran JV
    224. Mulder N
    225. Pollara VJ
    226. Ponting CP
    227. Schuler G
    228. Schultz J
    229. Slater G
    230. Smit AF
    231. Stupka E
    232. Szustakowki J
    233. Thierry-Mieg D
    234. Thierry-Mieg J
    235. Wagner L
    236. Wallis J
    237. Wheeler R
    238. Williams A
    239. Wolf YI
    240. Wolfe KH
    241. Yang SP
    242. Yeh RF
    243. Collins F
    244. Guyer MS
    245. Peterson J
    246. Felsenfeld A
    247. Wetterstrand KA
    248. Patrinos A
    249. Morgan MJ
    250. de Jong P
    251. Catanese JJ
    252. Osoegawa K
    253. Shizuya H
    254. Choi S
    255. Chen YJ
    256. Szustakowki J
    257. International Human Genome Sequencing Consortium
    (2001) Initial sequencing and analysis of the human genome
    Nature 409:860–921.
    https://doi.org/10.1038/35057062

Decision letter

  1. Deborah Bourc'his
    Reviewing Editor; Institut Curie, France
  2. Molly Przeworski
    Senior Editor; Columbia University, United States
  3. Michael Imbeault
    Reviewer; École Polytechnique Fédérale de Lausanne, Switzerland
  4. Geoffrey Faulkner
    Reviewer

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Transcriptional states of retroelement-inserted regions and specific KRAB zinc finger protein association are correlated with DNA methylation of retroelements in human male germ cells" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Molly Przeworski as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Michael Imbeault (Reviewer #1); Geoffrey Faulkner (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The reviewers agreed the study was well conducted and although was initially done on a limited number of available samples, the important observation regarding the inter-individual variability of SVA methylation in sperm was replicated in an additional dataset with more donors and further confirmed by targeted amplicon sequencing. The conclusions regarding the DNA methylation dependency towards specific KZFPs is not novel, although it was extended here to the context of primordial germ cells. The most intriguing conclusion is certainly the description of inter-individual differences in SVA methylation in the sperm of different donors. A link with being positioned inside a highly transcribed gene and in reverse orientation towards the host gene was drawn regarding the propensity of being methylated for an SVA, but the explanation for this is very uncertain. Results appear as preliminary on this matter.

Please address the following points in a revised version of your manuscript, and answer the individual comments and questions raised by the reviewers in a detailed rebuttal letter, with possible corrections and additional figures to be included in the manuscript.

1 – Integrate all the KZFP datasets that are available in GSE120539 (but that were not published) in your analysis

2 – Is it know whether interindividual differences in SVA methylation are also found in somatic tissues, or is it a specific feature of sperm DNA? If this has not been described, please analyze available WGBS datasets focusing on one tissue of several donors, or use your amplicon sequencing strategy. This would be an important addition to the biological meaning of such differences, whether this variability also occurs during embryonic reprogramming, not only germ cell reprogramming.

3 – Provide more information as to the methodology for full length retroelement analysis and ages of sperm donors

4 – Please show screen shots of individual loci that follow the stated correlation of DNA methylation and KZFP binding

5 – Review all sentences with over statements, as outlined by reviewer #1.

Additionally, please correct the following points:

– Line 304: "… subsequently, MIWI2-interacting protein SPOCD1 recruits the chromatin remodeling complex DNMT3A and DNMT3L to…". This sentence is wrong twice: (1) DNMT3A-DNMT3L does not have chromatin remodeling ability per se, it is a de novo DNA methylation complex, and (2) DNMT3C methylates piRNA-targeted retroelements in mice, not DNMT3A (see Barau et al., 2016). Please correct as you are referring here to the process that has been described in mice.

– The link between intragenic position of SVAs and their ability to undergo methylation by transcriptional readthrough/H3K36 methylation is interesting, as it means it would happen in a piRNA-independent manner. However, it does not explain why only the SVA elements with a reverse orientation to the host gene would be affected. Please report this.

Reviewer #1 (Recommendations for the authors):

I am a researcher active in the field of KRAB zinc finger proteins for about 12 years. In my opinion, the science in this manuscript is high quality and the findings as novel as they are interesting.

I have a few small recommendations to improve some sections of the manuscript and would appreciate if they can be implemented before publication.

Early in the manuscript it is stated that "we focused on full-length copies of retroelements to analyze 85 DNA methylation for at least 30 copies.". What did you consider as 'full-length', especially for LTR-containing retroviruses? Did you analyze only copies containing the internal part, discarding solo-LTRs? I could not find details describing this in the methods section. Also, this decision is probably biasing the analysis toward younger elements – not a problem in itself, but it should be stated clearly.

With very young elements (L1Hs) the mappability is probably not very good with 100bp non-paired end reads. Could you provide a supplemental table with average mappability per family of transposons, or as a supplemental figure as a violinplot of mappability per family.

There's more data of KZFP that is available that was not included – notably at GEO accession GSE120539 from the Trono lab – these are from the same experimental series as the ones published initially but didn't make it through analysis before publication – it would be great if you can include them.

For figure 2, I would like to see statistics of enrichment (p-value) for overlaps with specific KZFPs / families of transposons and DNA methylation categories.

Figure 2A – a heatmap is not the best visualization here – it could be a simple sorted barchart of hits zoomed in on the first few members with the highest scores. Same comment for figure 3H.

Review all sentences that are related to conclusions to avoid overly strong wording, considering that most findings of this manuscript are purely correlative and causation has not been demonstrated. As an example (out of many): "Therefore, SVAs are methylated during spermatogenesis if these are inserted into actively transcribed genes." – this is too strong, you might want to add 'suggest, might, potentially'.

Please discuss somewhere in the manuscript the potential for multiple KZFPs binding the same elements having a concerted effect on the elements.

Finally, I would like to see the age distribution of the sperm donors, and some analysis to see if variability is correlating with age in any way.

Reviewer #2 (Recommendations for the authors):

My two major reservations are the claims around inter-individual variability being difficult to distinguish from technical variablity, which I don't have a reasonable suggestion for how to address, and the first specific point above, namely that uniquely mapped WGBS reads are unlikely to measure methylation in the core VNTR region of an SVA (or the L1HS 5'UTR CpG island). The authors could address this point by showing a composite profile of WGBS coverage and methylation levels compared to L1/SVA consensus sequences, showing the inner parts of L1 and SVA. They could also show examples of individual loci that follow the stated patterns of DNA methylation and ZFP binding that supports a correlation between the two. Another option would be to do nanopore long-read sequencing, which obviously would take time and substantial resources, but would provide a comprehensive picture of the situation. Note that this issue affects some of the more high profile mouse retrotransposons, such as IAP, to a much lesser degree because their LTRs are more accessible to WGBS.

I think the figure legends for Figure 2E and Figure 2F are swapped.

https://doi.org/10.7554/eLife.76822.sa1

Author response

1 – Integrate all the KZFP datasets that are available in GSE120539 (but that were not published) in your analysis

Thanks for your suggestion. We reanalyzed KZFP peak data from GSE120539 and found that ZNF850 was enriched in highly methylated LTR12C copies in hPGCs. We included the analysis of ZNF850 in Supplementary Figure S4.

2 – Is it know whether interindividual differences in SVA methylation are also found in somatic tissues, or is it a specific feature of sperm DNA? If this has not been described, please analyze available WGBS datasets focusing on one tissue of several donors, or use your amplicon sequencing strategy. This would be an important addition to the biological meaning of such differences, whether this variability also occurs during embryonic reprogramming, not only germ cell reprogramming.

Thanks for the suggestion. This is also important point.

To address this issue, we reanalyzed IHEC DNA methylation data from adult skeletal muscle (15 samples) and CD4-positive helper T cells (18 samples). As all SVA groups showed high DNA methylation in all individuals tested, inter-individual variation of SVA DNA methylation in sperm would be cancelled during development (Supplementary Figure S5A, B). We stated this result in the text.

“In contrast to the inter-individual epigenetic variation of SVAs in sperm, a reanalysis of WGBS data of adult skeletal muscle from 15 individuals and of helper CD4-positive T cells from 18 individuals, which was deposited in the International Human Epigenome Consortium (IHEC) portal, showed high DNA methylation of SVAs in all individuals (Bujold et al., 2016) (Figure 6—figure supplement 1A, B). Thus, inter-individual variations in DNA methylation of SVAs in sperm are essentially canceled during development.” P15. line 276-281.

3 – Provide more information as to the methodology for full length retroelement analysis and ages of sperm donors

We added following sentences in the method section for the full length retroelement analysis.

“Only nearly full-length retroelements, whose length is 90% or more of the length of the consensus sequence of each subtype, were used for DNA methylation analysis of retroelements. We also included solo-LTR transposons in the DNA methylation analysis if they also possessed more than 90% of the consensus LTR sequence.” P34, line 606-609.

The ages of sperm donors used in Hammoud et al., data were described in the manuscript. On the other hand, we could not access clinical information of our samples for amplicon-seq, so age data for Okae et al., was not available.

4 – Please show screen shots of individual loci that follow the stated correlation of DNA methylation and KZFP binding

I added the representative loci in Figure S1J, S4G, 3H.

5 – Review all sentences with over statements, as outlined by reviewer #1.

We revised following sentences to be not over described as possible as we can.

Therefore, SVAs are methylated during spermatogenesis if these are inserted into actively transcribed genes.

“Therefore, SVAs located in actively transcribed regions in the antisense orientation are efficiently subjected to de novo DNA methylation during spermatogenesis.” p14, line 245-246

The data indicate that the effectiveness of the de novo DNA methylation varies among individuals.

“These results implicate the possibility that SVAs acquire DNA methylation during DNA methylation via transcription-directed machinery, and that the effectiveness of de novo DNA methylation varies among individuals.” p15, line 255-257

SVAs constitute a major source of inter-individual epigenetic variation in sperm

“SVAs are a potential source of inter-individual epigenetic variation in sperm”

p15, line 259

In this study, we showed that KRAB-ZFPs are associated with the DNA demethylation resistance of retroelements in hPGCs.

“Thus, the involvement of other KRAB-ZFPs in DNA demethylation resistance of retroelements in hPGCs is possible.” p19, line 330-331

SVAs are subjected to de novo DNA methylation by the transcription-directed DNA methylation machinery.

“the SVAs located in transcription-active regions in the antisense orientation are prone to methylation during spermatogenesis, which implies that the transcription-directed DNA methylation machinery might contribute to de novo DNA methylation of SVAs in male germ cells.” p18, line 302-304

Our data showed that the acquisition of DNA methylation of SVAs during spermatogenesis was regulated by the inserted regions and not by the nucleotide sequence. SVAs inserted in transcriptionally active regions in the antisense direction were targeted by de novo DNA methylation during spermatogenesis.

“Our data showed that the acquisition of DNA methylation of SVAs during spermatogenesis correlated with the inserted regions and not with the nucleotide sequence. SVAs inserted in transcriptionally active regions in the antisense direction are efficiently targeted for de novo DNA methylation during spermatogenesis.” P20, line 354- p21, 357

Additionally, please correct the following points:

– Line 304: "… subsequently, MIWI2-interacting protein SPOCD1 recruits the chromatin remodeling complex DNMT3A and DNMT3L to…". This sentence is wrong twice: (1) DNMT3A-DNMT3L does not have chromatin remodeling ability per se, it is a de novo DNA methylation complex, and (2) DNMT3C methylates piRNA-targeted retroelements in mice, not DNMT3A (see Barau et al., 2016). Please correct as you are referring here to the process that has been described in mice.

Thanks for the suggestion. We changed the sentences as bellow.

“In mouse spermatogenesis, MIWI2 binds piRNAs and is recruited to the nascent transcribed regions that are complementary to piRNAs (Watanabe et al., 2018). Subsequently, MIWI2 interacting protein SPOCD1, which forms a complex with DNMT3A and DNMT3L, and potentially with a rodent-specific DNA methyltransferase DNMT3C (Barau et al., 2016), induces DNA methylation on transposons (Zoch et al., 2020)..” P21, line 357- 361.

– The link between intragenic position of SVAs and their ability to undergo methylation by transcriptional readthrough/H3K36 methylation is interesting, as it means it would happen in a piRNA-independent manner. However, it does not explain why only the SVA elements with a reverse orientation to the host gene would be affected. Please report this.

Thank you for this comment. We described this issue as bellow.

“Furthermore, antisense RNAs embedded within protein-coding genes are selectively silenced by H3K36 methyltransferase SET2 in Saccharomyces cerevisiae (Venkatesh, Li, Gogol, and Workman, 2016). Thus, the machinery for gene body DNA methylation regulated by SETD2 is also a candidate for the de novo DNA methylation of SVAs during spermatogenesis.” p22, line 379-383.

Reviewer #1 (Recommendations for the authors):

I am a researcher active in the field of KRAB zinc finger proteins for about 12 years. In my opinion, the science in this manuscript is high quality and the findings as novel as they are interesting.

I have a few small recommendations to improve some sections of the manuscript and would appreciate if they can be implemented before publication.

Early in the manuscript it is stated that "we focused on full-length copies of retroelements to analyze 85 DNA methylation for at least 30 copies.". What did you consider as 'full-length', especially for LTR-containing retroviruses? Did you analyze only copies containing the internal part, discarding solo-LTRs? I could not find details describing this in the methods section. Also, this decision is probably biasing the analysis toward younger elements – not a problem in itself, but it should be stated clearly.

Thank you for your valuable and most KRAB-ZFP expert comments on our manuscript. As already pointed out, we also included solo-LTR copies in our DNA methylation analysis, if they possess more than 90% of the consensus LTR sequence. Length of consensus sequence was obtained from RepeatMasker data in UCSC genome browser. We described this in the Method section.

In addition, we stated our analysis was biased toward younger elements in Discussion section as bellow.

“Because we targeted full-length transposons, our analysis was biased towards relatively young transposons. Thus, it is possible that some fragmented older transposons may also be resistant to DNA demethylation in hPGCs.” P17, line 309-311.

With very young elements (L1Hs) the mappability is probably not very good with 100bp non-paired end reads. Could you provide a supplemental table with average mappability per family of transposons, or as a supplemental figure as a violinplot of mappability per family.

Thanks for the suggestion. To investigate mappability of WGBS read, we generated 100-bp sequences from each full-length copy in silico and mapped the reads to human genome by bismark. As the reviewer pointed out, the younger transposon, the lower mapping efficiency. This result was included in Supplementary Figure S1A.

There's more data of KZFP that is available that was not included – notably at GEO accession GSE120539 from the Trono lab – these are from the same experimental series as the ones published initially but didn't make it through analysis before publication – it would be great if you can include them.

Thanks for your good suggestion. We reanalyzed KZFP peak data from GSE120539 and found that ZNF850 was enriched in highly methylated LTR12C copies in hPGCs. We included the analysis of ZNF850 in Supplementary Figure S4. As we found KRAB-ZFP associated with LTR12 methylation, we removed old Figure S4 and Figure 3H.

For figure 2, I would like to see statistics of enrichment (p-value) for overlaps with specific KZFPs / families of transposons and DNA methylation categories.

Expected and observed numbers of SVA_A and LTR12C that overlap with KRAB-ZNF binding peaks were calculated and P-value for enrichment of retroelement in KRAB-ZNF binding peaks was calculated by prop.test by R. We included these data in Supplementary Figure S1C.

Figure 2A – a heatmap is not the best visualization here – it could be a simple sorted barchart of hits zoomed in on the first few members with the highest scores. Same comment for figure 3H.

Thanks for the comment.

We changed the heatmap in Figure 2A and Figure 3H to scatter plot.

Review all sentences that are related to conclusions to avoid overly strong wording, considering that most findings of this manuscript are purely correlative and causation has not been demonstrated. As an example (out of many): "Therefore, SVAs are methylated during spermatogenesis if these are inserted into actively transcribed genes." – this is too strong, you might want to add 'suggest, might, potentially'.

Thanks for the suggestion. We revised following sentences to be not over described as possible as we can.

Therefore, SVAs are methylated during spermatogenesis if these are inserted into actively transcribed genes.

“Therefore, SVAs located in actively transcribed regions in the antisense orientation are efficiently subjected to de novo DNA methylation during spermatogenesis.” p14, line 245-246

The data indicate that the effectiveness of the de novo DNA methylation varies among individuals.

“These results implicate the possibility that SVAs acquire DNA methylation during DNA methylation via transcription-directed machinery, and that the effectiveness of de novo DNA methylation varies among individuals.” p15, line 255-257

SVAs constitute a major source of inter-individual epigenetic variation in sperm

“SVAs are a potential source of inter-individual epigenetic variation in sperm”

p15, line 259

In this study, we showed that KRAB-ZFPs are associated with the DNA demethylation resistance of retroelements in hPGCs.

“Thus, the involvement of other KRAB-ZFPs in DNA demethylation resistance of retroelements in hPGCs is possible.” p19, line 330-331

SVAs are subjected to de novo DNA methylation by the transcription-directed DNA methylation machinery.

“the SVAs located in transcription-active regions in the antisense orientation are prone to methylation during spermatogenesis, which implies that the transcription-directed DNA methylation machinery might contribute to de novo DNA methylation of SVAs in male germ cells.” p18, line 302-304

Our data showed that the acquisition of DNA methylation of SVAs during spermatogenesis was regulated by the inserted regions and not by the nucleotide sequence. SVAs inserted in transcriptionally active regions in the antisense direction were targeted by de novo DNA methylation during spermatogenesis.

“Our data showed that the acquisition of DNA methylation of SVAs during spermatogenesis correlated with the inserted regions and not with the nucleotide sequence. SVAs inserted in transcriptionally active regions in the antisense direction are efficiently targeted for de novo DNA methylation during spermatogenesis.” P20, line 354- p21, 357

Please discuss somewhere in the manuscript the potential for multiple KZFPs binding the same elements having a concerted effect on the elements.

We added the following sentence in the Discussion section.

“In hPGCs, multiple KRAB-ZNPs were correlated with DNA demethylation resistance in the same retroelements, which may contribute to more robust or cooperative retroelement suppression.” P18, line 319- p19, line 321

Finally, I would like to see the age distribution of the sperm donors, and some analysis to see if variability is correlating with age in any way.

The ages of sperm donors used in Hammoud et al., data were described in the manuscript (Donor 1 -32 yo, Donor2 -37 yo). On the other hand, we could not access clinical information of our samples for amplicon-seq, so age data for Okae et al., was not available. To investigate the correlation of SVA methylation with age, we further analyzed publicly available sperm WGBS data from five healthy individuals and six oligozoospermic patients, which are deposited in the European Nucleotide Archive (ENA) under the accession number PRJEB34432 (Leitao et al., Clinical Epigenetics, 2020) and are linked to clinical information including age. We found that both healthy and patients showed variable SVA methylation and age was not associated with SVA methylation state. We also found that testosterone level in blood was negatively correlated with SVA methylation in sperm, implicating potential link between SVA methylation in sperm and physiological conditions. These data were included in Supplementary Figure S5C, D and described in the text as follow.

“Finally, to investigate whether inter-individual DNA methylation variations are associated with physiological or disease conditions, we reanalyzed publicly available WGBS data from five healthy donors and six oligozoospermic patients (European Nucleotide Archive (ENA) under the accession number PRJEB34432) (Leitao et al., 2020). The disease condition was not associated with inter-individual DNA methylation variations in SVAs, because both healthy donors and oligozoospermic patients showed inter-individual variations of DNA methylation in “High and Low” SVA_B–F copies (Figure 6—figure supplement 1C). On the other hand, a comparison of various physiological conditions between highly methylated individuals and lowly methylated ones (median methylation levels of “High and Low” > 50% vs < 50%) revealed that blood testosterone levels were significantly higher in lowly methylated individuals than in highly methylated ones (Figure 6—figure supplement 1D). However, prolactin, follicle stimulating hormone (FSH), luteinizing hormone (LH), sex hormone-binding globulin (SHBG) blood levels and age were not significantly different between the two groups (Figure 6—figure supplement 1D). Although further validation of this correlation is required, DNA methylation of SVAs in sperm may be associated with physiological conditions.” P16, line 282- P17, line 296.

Reviewer #2 (Recommendations for the authors):

My two major reservations are the claims around inter-individual variability being difficult to distinguish from technical variablity, which I don't have a reasonable suggestion for how to address, and the first specific point above, namely that uniquely mapped WGBS reads are unlikely to measure methylation in the core VNTR region of an SVA (or the L1HS 5'UTR CpG island). The authors could address this point by showing a composite profile of WGBS coverage and methylation levels compared to L1/SVA consensus sequences, showing the inner parts of L1 and SVA. They could also show examples of individual loci that follow the stated patterns of DNA methylation and ZFP binding that supports a correlation between the two. Another option would be to do nanopore long-read sequencing, which obviously would take time and substantial resources, but would provide a comprehensive picture of the situation. Note that this issue affects some of the more high profile mouse retrotransposons, such as IAP, to a much lesser degree because their LTRs are more accessible to WGBS.

Thanks for your valuable comment/suggestion. As the reviewer pointing out, it is difficult to measure DNA methylation in core VNTR regions of a younger SVA. But we were able to do it in the oldest SVA type, SVA_A. We showed the DNA methylation levels in VNTR regions of SVA_A and overlap with KRAB-ZFP peaks in new Supplementary Figure S1H. We also showed examples of individual loci showing correlation of DNA methylation in hPGCs and KRAB ZFP association in SVA_A and L1PA4 in new Supplementary Figure S1I,J and 3H. And we described about this issue as bellow.

“We also confirmed that DNA methylation levels within the VNTR were correlated with ZNF257 or ZNF28 association (Figure 2—figure supplement 1G, H).” p9, 149-151

“We also confirmed high DNA methylation in the ZNF649 binding motifs at individual loci (Figure 3H).” p11, 185-186

I think the figure legends for Figure 2E and Figure 2F are swapped.

Thank you. We fixed it.

https://doi.org/10.7554/eLife.76822.sa2

Article and author information

Author details

  1. Kei Fukuda

    Cellular Memory Laboratory, RIKEN Cluster for Pioneering Research, Wako, Japan
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Validation, Writing - original draft, Writing - review and editing
    For correspondence
    kei.fukuda@riken.jp
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4111-9409
  2. Yoshinori Makino

    Laboratory of Pathology and Development, Institute for Quantitative Biosciences, The University of Tokyo, Tokyo, Japan
    Contribution
    Resources
    Competing interests
    No competing interests declared
  3. Satoru Kaneko

    Department of Obstetrics and Gynecology, Ichikawa General Hospital, Tokyo Dental College, Ichikawa, Japan
    Contribution
    Resources
    Competing interests
    No competing interests declared
  4. Chikako Shimura

    Cellular Memory Laboratory, RIKEN Cluster for Pioneering Research, Wako, Japan
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  5. Yuki Okada

    Laboratory of Pathology and Development, Institute for Quantitative Biosciences, The University of Tokyo, Tokyo, Japan
    Contribution
    Resources
    Competing interests
    No competing interests declared
  6. Kenji Ichiyanagi

    Laboratory of Genome and Epigenome Dynamics, Department of Animal Sciences, Graduate School of Bioagricultural Sciences, Nagoya University, Nagoya, Japan
    Contribution
    Conceptualization
    Competing interests
    No competing interests declared
  7. Yoichi Shinkai

    Cellular Memory Laboratory, RIKEN Cluster for Pioneering Research, Wako, Japan
    Contribution
    Conceptualization, Funding acquisition, Resources, Supervision, Writing - original draft, Writing - review and editing
    For correspondence
    yshinkai@riken.jp
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6051-2484

Funding

Japan Society for the Promotion of Science (18H05530)

  • Yoichi Shinkai

Japan Society for the Promotion of Science (18H03991)

  • Yoichi Shinkai

RIKEN (SPDR)

  • Kei Fukuda

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank the staff of the Support Unit for Bio-Material Analysis (BMA) at the RIKEN Center for Brain Science (CBS) Research Resources Division (RRD) for NGS library construction. This research was supported by the Special Postdoctoral Researcher (SPDR) Program of RIKEN to KF, the Japan Ministry of Education, Culture, Sports, Science, and Technology Grant-in-Aid for Scientific Research (18H05530, 18H03991) to YS, and RIKEN internal research funds to YS.

Ethics

Human subjects: This study was approved by the ethics committees of RIKEN, Tokyo University, and Ichikawa General Hospital.All study participants were briefed about the aims of the study and the parameters to be measured, and consent was obtained.

Senior Editor

  1. Molly Przeworski, Columbia University, United States

Reviewing Editor

  1. Deborah Bourc'his, Institut Curie, France

Reviewers

  1. Michael Imbeault, École Polytechnique Fédérale de Lausanne, Switzerland
  2. Geoffrey Faulkner

Publication history

  1. Preprint posted: May 19, 2021 (view preprint)
  2. Received: January 6, 2022
  3. Accepted: March 21, 2022
  4. Accepted Manuscript published: March 22, 2022 (version 1)
  5. Version of Record published: March 30, 2022 (version 2)

Copyright

© 2022, Fukuda et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 860
    Page views
  • 155
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Kei Fukuda
  2. Yoshinori Makino
  3. Satoru Kaneko
  4. Chikako Shimura
  5. Yuki Okada
  6. Kenji Ichiyanagi
  7. Yoichi Shinkai
(2022)
Potential role of KRAB-ZFP binding and transcriptional states on DNA methylation of retroelements in human male germ cells
eLife 11:e76822.
https://doi.org/10.7554/eLife.76822
  1. Further reading

Further reading

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics
    Meng Huang, Minjie Hong ... Xuezhu Feng
    Research Article Updated

    Histone methylation plays crucial roles in the development, gene regulation, and maintenance of stem cell pluripotency in mammals. Recent work shows that histone methylation is associated with aging, yet the underlying mechanism remains unclear. In this work, we identified a class of putative histone 3 lysine 9 mono/dimethyltransferase genes (met-2, set-6, set-19, set-20, set-21, set-32, and set-33), mutations in which induce synergistic lifespan extension in the long-lived DAF-2 (insulin growth factor 1 [IGF-1] receptor) mutant in Caenorhabditis elegans. These putative histone methyltransferase plus daf-2 double mutants not only exhibited an average lifespan nearly three times that of wild-type animals and a maximal lifespan of approximately 100 days, but also significantly increased resistance to oxidative and heat stress. Synergistic lifespan extension depends on the transcription factor DAF-16 (FOXO). mRNA-seq experiments revealed that the mRNA levels of DAF-16 Class I genes, which are activated by DAF-16, were further elevated in the daf-2;set double mutants. Among these genes, tts-1, F35E8.7, ins-35, nhr-62, sod-3, asm-2, and Y39G8B.7 are required for the lifespan extension of the daf-2;set-21 double mutant. In addition, treating daf-2 animals with the H3K9me1/2 methyltransferase G9a inhibitor also extends lifespan and increases stress resistance. Therefore, investigation of DAF-2 and H3K9me1/2 deficiency-mediated synergistic longevity will contribute to a better understanding of the molecular mechanisms of aging and therapeutic applications.

    1. Chromosomes and Gene Expression
    2. Structural Biology and Molecular Biophysics
    Lena Maria Muckenfuss, Anabel Carmen Migenda Herranz ... Martin Jinek
    Research Article Updated

    3′ end formation of most eukaryotic mRNAs is dependent on the assembly of a ~1.5 MDa multiprotein complex, that catalyzes the coupled reaction of pre-mRNA cleavage and polyadenylation. In mammals, the cleavage and polyadenylation specificity factor (CPSF) constitutes the core of the 3′ end processing machinery onto which the remaining factors, including cleavage stimulation factor (CstF) and poly(A) polymerase (PAP), assemble. These interactions are mediated by Fip1, a CPSF subunit characterized by high degree of intrinsic disorder. Here, we report two crystal structures revealing the interactions of human Fip1 (hFip1) with CPSF30 and CstF77. We demonstrate that CPSF contains two copies of hFip1, each binding to the zinc finger (ZF) domains 4 and 5 of CPSF30. Using polyadenylation assays we show that the two hFip1 copies are functionally redundant in recruiting one copy of PAP, thereby increasing the processivity of RNA polyadenylation. We further show that the interaction between hFip1 and CstF77 is mediated via a short motif in the N-terminal ‘acidic’ region of hFip1. In turn, CstF77 competitively inhibits CPSF-dependent PAP recruitment and 3′ polyadenylation. Taken together, these results provide a structural basis for the multivalent scaffolding and regulatory functions of hFip1 in 3′ end processing.