Genome-wide Bisulfite Sequencing in Zygotes Identifies Demethylation Targets and Maps the Contribution of TET3 Oxidation

Summary Fertilization triggers global erasure of paternal 5-methylcytosine as part of epigenetic reprogramming during the transition from gametic specialization to totipotency. This involves oxidation by TET3, but our understanding of its targets and the wider context of demethylation is limited to a small fraction of the genome. We employed an optimized bisulfite strategy to generate genome-wide methylation profiles of control and TET3-deficient zygotes, using SNPs to access paternal alleles. This revealed that in addition to pervasive removal from intergenic sequences and most retrotransposons, gene bodies constitute a major target of zygotic demethylation. Methylation loss is associated with zygotic genome activation and at gene bodies is also linked to increased transcriptional noise in early development. Our data map the primary contribution of oxidative demethylation to a subset of gene bodies and intergenic sequences and implicate redundant pathways at many loci. Unexpectedly, we demonstrate that TET3 activity also protects certain CpG islands against methylation buildup.


In Brief
Peat et al. employ whole-genome bisulfite sequencing to generate a global picture of the profound DNA methylation reprogramming that is triggered by fertilization. This identifies gene bodies as a major demethylation target where TET3 oxidation plays an important role and uncovers a protective TET3 function as well as redundancy in the demethylation machinery.

Accession Numbers GSE63417
INTRODUCTION Over a decade ago, pioneering studies demonstrated that fertilization triggers a global and active loss of DNA methylation from the paternal genome, but not its maternal counterpart Oswald et al., 2000). This methylation reprogramming occurs at a pivotal developmental time point, when a global transcriptional transition is required for the genesis of totipotency from specialized germ cell states (Hemberger et al., 2009). The role of methylation changes is still unclear, and a ''resetting'' of the epigenome may be necessary to generate a blank canvas on which to paint new regulatory marks for the development of the embryo.
The wave of zygotic erasure does not affect all regions of the paternal genome equally, as methylation at certain sequences, such as parental imprints and active retrotransposons (e.g., intracisternal A particles [IAPs]), must be maintained for embryonic viability (Seisenberger et al., 2012). Recent genome-scale profiling by reduced-representation bisulfite sequencing (RRBSseq) extended observations of variation in methylation dynamics across different classes of repetitive elements and even within their component families (Smith et al., 2012). However, as the great majority of 5-methylcytosine (5mC) in sperm lies outside the CpG islands (CGIs) enriched in RRBS-seq , which are largely hypomethylated (Smallwood et al., 2011), a comprehensive characterization of zygotic demethylation targets is still lacking. Such a profile would be particularly instructive in the context of a major transcriptional and fate transition, and would complement recent genome-wide studies in gametes and later preimplantation stages Wang et al., 2014).
We and others have demonstrated that oxidation of the paternal genome by TET3 constitutes an important mechanism for the active removal of its methylation in the zygote (Gu et al., 2011;Wossidlo et al., 2011), adding to the small repertoire of factors that are known to contribute to demethylation in the zygote, including the elongator complex (Okada et al., 2010) and the deaminase AID . However, the knowledge of these factors' involvement has been set against a limited understanding of normal demethylation dynamics across the genome, and an appreciation of the role of specific pathways that is restricted to a handful of loci (Gu et al., 2011;Hajkova et al., 2010;Okada et al., 2010). Therefore, the architecture of the demethylation machinery and the relationship between different mechanisms (synergy, redundancy, or specificity) remain almost entirely unknown.
Here, we set out to produce a comprehensive picture of methylation dynamics in the paternal pronucleus and thus provide the foundation on which to map the contribution of oxidation using genetic disruption of TET3 activity in the zygote.

RESULTS
In order to obtain genome-wide profiles from limited material, we optimized a strategy for generating whole-genome bisulfite sequencing (WGBS-seq) libraries wherein bisulfite treatment is used to both convert cytosines and fragment DNA at the start of the protocol . This enhances yield by eliminating the need for fragmentation by sonication and avoiding the degradation of adaptor-tagged sequences. We made substantial modifications (for the detailed protocol, see Supplemental Experimental Procedures) that allowed us to consistently generate high-quality libraries from fewer than 300 zygotes.
This enhanced WGBS-seq protocol was applied to zygotes collected from control females and females carrying a conditional deletion that abolishes TET3 activity in the germline (referred to as TET3 zygotes; see Supplemental Experimental Procedures and sequencing statistics in Table S1). 129S2/ SvHsd studs were used in matings in order to specifically trace the fate of paternal methylation in the zygote with SNPs. To permit examination of demethylation by oxidation as well as dilution at DNA replication, we collected late-stage zygotes (B) Frequency distribution of 20 kb tile methylation values. Tiles were allocated to bins centered every 5%. (C) Mean methylation levels in repetitive-element classes. ''Zygotes'' refers to overall levels, as natural sequence variability prohibits assignment of parental alleles using SNPs. LINE, long interspersed nuclear element; SINE, short interspersed nuclear element; MuLV, murine leukemia virus; IAP, intracisternal A particle; Etn, early transposon; MERV-L, murine endogenous retrovirus with leucine tRNA primer; ERV, endogenous retrovirus. See also Figure S1. that had completed S phase, after the global loss of 5mC and gain in 5-hydroxymethylcytosine (5hmC) had occurred Wossidlo et al., 2011). The concordance of imprint control region methylation levels with expected values indicated that the data sets generated by this strategy were of high quality ( Figure S1A).

The Zygotic Methylation Landscape
We first used 20 kb tiles to assess global 5mC values ( Figure 1A). The results demonstrated a large decrease in mean methylation from sperm (85.4%) to the paternal pronucleus (50.8%), in line with genome-wide demethylation. Methylation levels were slightly but significantly increased in TET3 zygotes (44.9% versus 46.8%) and this effect was more marked in paternal DNA (50.8% versus 58.5%), consistent with the role of TET3 in promoting paternal demethylation as shown by immunofluorescence and locus-specific bisulfite analysis (Gu et al., 2011;Wossidlo et al., 2011).
The distribution of methylation values in sperm and paternal pronuclear DNA further established a dramatic reprogramming of methylation patterns at fertilization: whereas the majority of sequences were highly methylated in sperm, few tiles (7.5%) remained over 80% methylated in the zygote ( Figure 1B), demonstrating that most of the paternal genome was affected by some methylation loss. The bulk of sequences in the paternal pronucleus possessed intermediate methylation in a broad distribution that shifted upward significantly in the absence of TET3 oxidation.
Examination of major repetitive-element classes confirmed previous observations of a transition to a more hypomethylated state in zygotes relative to sperm, with the established exception of IAPs ( Figure 1C; Lane et al., 2003;Smith et al., 2012;Wossidlo et al., 2010). The use of WGBS-seq allowed this analysis to be extended to largely uncharacterized elements that are poorly  covered by RRBS-seq due to the paucity of CCGG sites used for enrichment. This demonstrated that MERV-L elements (which contain only four CCGG sites in 6.5 kb) undergo substantial demethylation in the zygote. Removal of methylation may be important for the strong transcriptional activation of these repeats at fertilization (Kigami et al., 2003;Macfarlan et al., 2012). Furthermore, our data show that methylation at satellite repeats, which lack CCGG sites altogether, is relatively stable during zygotic reprogramming. This maintenance may be important for the heterochromatinization that is required for the first mitosis to proceed (Fadloun et al., 2013;Probst and Almouzni, 2011). These insights into methylation trajectories in the zygote underscore the value of the WGBS-seq approach. The absence of oxidation by TET3 had a limited impact on the methylation levels of repetitive elements (none were altered by more than 10%). This indicates either that TET3 is not targeted to these sequences or that alternative demethylation pathways operate redundantly or are recruited to noncanonical targets to compensate for the absence of TET3 oxidation.
These analyses paint a picture of a significant but moderate impact of TET3 at a global level, consistent with immunofluorescence data , and suggest that the major role of oxidative demethylation is at single-copy loci.

Methylation Trajectories in Different Genomic Contexts
We next examined methylation dynamics at specific genomic features more closely. Given the particular biology of methylation at CGIs and transcription start sites (Jones, 2012), we further subcategorized the CGIs and promoters based on whether they overlapped with a promoter region or CGI, respectively.
An examination of methylation patterns first established that loss of paternal methylation at fertilization was significant at all sequence features, but occurred to varying degrees ( Figure 2A). A quantitative comparison of methylation levels in sperm and paternal DNA revealed that 90% of intergenic and gene body sequences were significantly demethylated (decrease R 10%, p < 0.05), with 70% losing >25% methylation ( Figure 2B). As these regions contribute the great majority of paternally methylated cytosines ( Figure S1B), this demonstrates the pervasive removal of methylation from single-copy sequences and further identifies gene bodies as a major target of zygotic demethylation. In contrast, little demethylation was observed at CGI-associated sequences ( Figure 2B), reflecting the general hypomethylation of these regions in sperm (Figures 2A and S2A). The minority that were methylated underwent a similar degree of demethylation compared with other sequences ( Figure S2B), indicating that paternal methylation is generally targeted for removal regardless of the context. Although they were a less prevalent feature and thus not a large carrier of paternal methylation, non-CGI promoters followed a trajectory similar to that observed for intergenic and gene body sequences (Figures 2A and 2B).
In order to delineate the contribution of replication to the observed demethylation, we expressed absolute demethylation as a proportion of the initial methylation in sperm to determine the relative demethylation level (RDL) for each locus ( Figure 2C). An RDL > 0.5 cannot be explained simply by a 2-fold reduction in methylation at replication, and therefore implies the action of active removal pathways. Demethylation at loci with RDL < 0.5 may also be active but cannot be discriminated from replicative dilution; therefore, this analysis measures the minimum number of actively demethylated loci. We applied a conservative threshold of RDL > 0.6 and found that the degree of methylation loss at 25% and 30% of demethylated intergenic and gene body sequences, respectively, which constitute the majority of demethylation targets, implicates an active process. This demonstrates that while lack of maintenance at replication may provide an important mechanism for loss of methylation, active demethylation of the paternal genome is common. Strikingly, although few CGI-associated sequences lost 5mC, where demethylation did occur, it was almost exclusively active.
By annotating each demethylated locus with the nearest gene, we examined the relationship between zygotic demethylation and transcription in 2-cell embryos using published mRNA-seq data , and searched for functional enrichment in Gene Ontology (GO) databases. With the exception of nonpromoter CGIs, genes associated with demethylated loci have elevated expression at the 2-cell stage ( Figure S2C) and as a group are enriched for cytoskeletal, ion transport, signaling, protein modification, and RNA processing terms (Table S2). This expression profile and the enrichment for several known ZGA transcriptional modules, including GTPase signaling and RNA processing (Xue et al., 2013), suggest a role for demethylation in preparing for transcriptional activation in the early embryo. The prevalence of cytoskeletal and cell junction terms may be linked to the function of these components during embryo compaction (Ducibella and Anderson, 1975;Fierro-Gonzá lez et al., 2013). Genes associated with demethylated nonpromoter CGIs are not highly expressed at ZGA and are enriched for a largely nonoverlapping set of GO terms, including many neural functions.
We next questioned whether pervasive loss of intragenic methylation at fertilization was linked to variation in transcription p value < 0.05). n = 5,393 intergenic sequences, 7,405 gene body sequences, 4,213 CGI promoters, 1,366 non-CGI promoters, 526 promoter CGIs, and 257 nonpromoter CGIs. (D) Distribution of early embryo paternal transcriptional noise levels at genes allocated into quartiles according to average intragenic methylation levels in the paternal pronucleus. Normalized transcriptional noise refers to expression-corrected variation in transcription from the paternal allele between single cells as calculated in Supplemental Experimental Procedures. The p values shown are the result of an ANOVA multiple-comparison test between indicated data sets with Sidak correction. The differences between second and third quartiles, and third and fourth quartiles are not significant at any stage according to this test. The total number of genes that fell into each quartile was as follows: 0%-25%, 1,379; 25%-50%, 2,468; 50%-75%, 2,605; and 75%-100%, 270. Genes that were not expressed at a particular stage were excluded from analysis of that stage. The number of genes analyzed at each stage is given in Supplemental Experimental Procedures. (E) Relationship between the average paternal pronuclear gene body methylation and paternal transcriptional noise at the late 2-cell stage. A quadratic equation is a significantly better fit to the data than either a linear or no relationship (p < 0.0001 in both cases, extra sum-of-squares F test). See also Figure S2 and Table S2. between cells at subsequent stages. Using published allele-specific, single-cell RNA-seq data , we generated a metric of cell-to-cell transcriptional noise at paternal alleles for each expressed gene during early embryonic development (see Supplemental Experimental Procedures). We then examined the distribution of noise signals at genes whose average intragenic methylation levels in the paternal pronucleus fell into each of four quartiles ( Figure 2D). Strikingly, this revealed that intermediate levels of gene body methylation (25%-75%) are associated with significantly higher transcriptional noise than hypomethylated genes (<25%) throughout the preimplantation embryo. Although the low number of hypermethylated genes (>75%) in the paternal pronucleus precludes establishment of statistical significance, their transcriptional noise levels appear substantially lower than those of intermediate-methylated genes. Indeed, the relationship between intragenic methylation and transcriptional noise is significantly better modeled by a quadratic equation than a straight line at all developmental stages examined ( Figures 2E and S2D). The widespread demethylation of paternally hypermethylated gene bodies to intermediate methylation levels in the zygote ( Figure 2A) may therefore promote transcriptional noise in the early embryo. As our noise metric normalizes for expression levels, this association is independent of any influence of intragenic methylation on overall transcription. This analysis uncovers a connection between methylation and early embryonic transcription at a major demethylation target, further emphasizing the importance of genome-wide approaches to explore methylation reprogramming and its functional impact.

The Absence of TET3 Activity Disrupts Normal Methylation Dynamics
Having established a detailed picture of zygote methylation, we sought to understand the contribution of oxidation by TET3 to this landscape. The absence of TET3 oxidation significantly elevated the paternal pronuclear methylation levels of all features at identified demethylation targets ( Figure 3A) as well as these groups as a whole ( Figure S3A). By quantifying the impact of TET3 deletion on demethylated loci, we established that loss of methylation was significantly impaired (gain % 10%, p < 0.05) at 25%-30% of intergenic, gene body, and non-CGI promoter sequences, but found relatively little effect on the small number of demethylated CGI-associated sequences (Figure 3B). As intergenic and gene body sequences constitute the bulk of demethylated loci, this identifies these regions as the major target of oxidative demethylation by TET3. Hierarchical clustering analysis of all methylation values showed that for gene body and non-CGI promoter sequences, but not other features, paternal pronuclear DNA from zygotes lacking TET3 activity formed a strongly supported clade with sperm DNA rather than paternal pronuclear DNA from control zygotes ( Figures 3C and S3B). This underlines the importance of TET3 for reprogramming of sperm methylation patterns at these features.
Nevertheless, our data demonstrate that although TET3 makes a significant contribution, demethylation proceeds normally at the majority of loci in its absence. This indicates that either most demethylated loci are not targeted by TET3 or redundant or compensatory mechanisms minimize the impact of its absence. At CGI-associated sequences, where demethylation is almost exclusively active ( Figure 2C), this implicates the activity of alternative enzymatic pathways such as the elongator complex or oxidation-independent BER.
In general, the targets of TET3 demethylation reflected the characteristics of demethylated loci with regard to their association with elevated expression in 2-cell embryos ( Figure S3C) and, where there were a sufficient number to test, enrichment for similar GO terms (Table S2). This suggests that although the importance of oxidative 5mC removal varies with the sequence feature, demethylation by TET3 is not associated with specific biological functions.
We next focused on the loci that were demethylated by TET3 and examined the degree to which its absence affected normal demethylation in order to gauge the contribution of oxidation to their methylation trajectory ( Figure 3D). At a substantial number of loci, particularly the small number of targeted CGI promoters (e.g., the Col14a1 intragenic sequence and Plbd1 CGI promoter), TET3 was responsible for all or most of the observed demethylation ( Figure 3E). However, at most loci, TET3 deletion resulted in partial impairment of demethylation (e.g., the Zfp356 intragenic sequence and Coro1c non-CGI promoter; Figure 3E), indicating that multiple pathways can act at these sequences-either in synergy with TET3 oxidation or to compensate (incompletely) for its loss. The prevalence of partial impairment demonstrates that such redundant targeting or compensatory ability is a widespread feature of the demethylation machinery.
This analysis further revealed that at some demethylated loci, paternal methylation in TET3 zygotes actually increased to higher levels than in sperm (i.e., value > 1 on Figure 3D), implying that beyond its function in demethylation, TET3 is required to (C) Dendrogram derived from hierarchical clustering using the Pearson correlation distance for the indicated sequence features. The number above the branch denotes the approximately unbiased p values computed by multiscale bootstrap resampling (n = 10,000) using the pvclust package for R . This was identical for both features. (D) Proportion of normal demethylation affected by the loss of TET3 at TET3-demethylated loci, calculated by dividing the absolute paternal methylation change between control and TET3 zygotes (as in B) by the normal absolute demethylation in control zygotes (as in Figure 2B). TET3-demethylated loci are defined as demethylated loci that underwent an absolute paternal methylation gain of R10% in the absence of TET3 and were significant according to a chi-square test (corrected p value < 0.05). n = 1,310 intergenic sequences, 2,487 gene body sequences, 259 CGI promoters, and 398 non-CGI promoters. The number of promoter and nonpromoter CGIs that met these criteria was insufficient for meaningful analysis. *Proportional change > 1 indicates that paternal methylation in TET3 zygotes increased to higher levels than in sperm. (E) Examples of complete (Col14a1 intragenic sequence and Plbd1 CGI promoter) and partial (Zfp356 intragenic sequence and Coro1c non-CGI promoter) impairment of demethylation due to loss of TET3 oxidation. See also Figure S3 and Table S2. protect certain sequences from aberrant de novo methylation in the zygote. Indeed, an examination of all loci that gained paternal methylation in TET3 zygotes (''TET3 targets'') suggests that this phenomenon occurs frequently at CGI-associated sequences ( Figures 4A and S4A).
By determining whether paternal methylation in TET3 zygotes increased to significantly higher levels than in sperm (R10%, p < 0.05), we generated a list of loci protected by TET3 to assess their normal methylation behavior and establish the relative importance of this protective function ( Figures 4B and 4C). Figure 4. A Protective Role for TET3

A C B
(A) Distribution of methylation values for TET3 targets at CGI-associated sequences. TET3 targets are defined as loci where the absence of TET3 activity resulted in a paternal methylation gain of R10%, which is significant according to a chi-square test (corrected p value < 0.05).
(B) Function of TET3 at its targets and their methylation behavior. The area of each circle is proportional to the percentage of loci in that category. TET3 targets are defined as in (A). Protected, paternal methylation R10% higher in TET3 zygotes than in sperm and significant according to a chi-square test (corrected p value < 0.05); demethylated, criteria as in Figure 2C; normal gain methylation, paternal methylation R10% higher in control zygotes than in sperm and significant according to a chi-square test (corrected p value < 0.05). Percentages do not necessarily add to 100 due to the application of statistical criteria.
(C) Examples of loci protected by TET3 classified according to whether they are normally demethylated, stable, or gain methylation. Agrn and Gm106 are CGI promoters, Lmtk3 is an intragenic CGI, Kif6 is a non-CGI promoter, and Tbc1d2b is a promoter CGI. Chr8 Orphan CGI is located at chromosome 8 29088955-29089878 (NCBI37 assembly). See also Figure S4.
Although a subset of protected sequences are demethylation targets (e.g., Agrn CGI promoter and orphan CGI), the majority do not normally undergo methylation loss at fertilization and instead are stably maintained (e.g., Lmtk3 intragenic CGI and Kif6 non-CGI promoter; Figure 4C) or gain some methylation (e.g., Tbc1d2b promoter CGI and Gm106 CGI promoter; Figure 4C). Thus, TET3 targets a distinct set of loci specifically to prevent accumulation of methylation. Although they constitute a small fraction (2.4%) of total TET3 targets, almost twice as many promoter CGIs are protected than demethylated. A roughly equal number of nonpromoter CGIs (which account for 1.8% of TET3 targets) fall into each category ( Figure 4B). Protection is therefore an important function at these sequence features, in contrast to intergenic, gene body, and non-CGI promoter sequences, where demethylation is the sole function of TET3 at >95% of loci ( Figure 4B). Like its demethylation targets, the loci that are protected by TET3 are associated with elevated expression in the early embryo ( Figure S4B), suggesting that TET3 prevents accumulation of methylation to safeguard transcriptional activation.

WGBS-Seq Characterizes Methylation Trajectories across the Genome
Here, we present genome-wide methylation profiles of zygotes, allowing a characterization of the extensive methylation reprogramming that occurs at this profound developmental transition. This fills an important gap in WGBS-seq studies tracking methylation dynamics from gametes to the blastocyst stage Wang et al., 2014), and expands the perspective from enrichment-based profiling of zygotes, including two reports published while this study was under review (Guo et al., 2014;Shen et al., 2014;Smith et al., 2012). By extensively optimizing a noncanonical WGBS-seq strategy, we were able to overcome the technical barrier of low cell numbers, which enabled us to access previously poorly characterized regions of the genome.
Our data paint a picture of global remodeling of sperm methylation patterns in which the major carriers of 5mC (repetitive elements, intergenic regions, and gene bodies) are extensively demethylated. The observation of pervasive removal of 5mC from gene bodies adds a new dimension to our understanding of methylation reprogramming, which up to now has centered on promoters and CGIs. Furthermore, the WGBS-seq profiles allowed us to establish the behavior of functionally important repetitive sequences lacking the CCGG sites used in RRBSseq, including the demethylation of MERV-L elements and stable maintenance of satellite methylation.
The degree of methylation loss from many sequences necessarily implicates active removal processes and shows that such mechanisms operate at a large number of loci. The role of active and passive mechanisms at sequences undergoing lower proportional demethylation (RDL < 0.5) is unclear. Dissecting their relative contribution across the genome will necessitate inhibition of replication coupled with WGBS-seq to examine the generality of observations of a large role for passive demethylation from RRBS-seq (Guo et al., 2014;Shen et al., 2014).
We find that loss of 5mC from all functional sequence features except nonpromoter CGIs is associated with transcriptional activation in the early embryo, as assessed by expression levels and functional enrichment. Clearly, this relationship requires experimental testing. Given the existence of multiple demethylation pathways in the zygote that may act redundantly, including replicative dilution, it may be challenging to perturb demethylation sufficiently to evaluate any transcriptional effect. The potential involvement of regulatory mechanisms in addition to DNA methylation that operate to promote transcriptional activation could further complicate interpretation. Indeed, maternal deletion of TET3 had little impact on 2-cell and inner cell mass transcriptomes (Shen et al., 2014).
A connection between zygotic demethylation of gene bodies and ZGA seems at odds with the positive correlation between intragenic methylation and active transcription observed in many tissues, including gametes (Jones, 2012;. However, recent studies suggest that this relationship is more nuanced. A meta-analysis in human cell lines found that intragenic methylation was highest at moderately expressed genes and low at both weakly and strongly expressed genes (Jjingo et al., 2012). Functionally, methylation can modulate the use of alternative promoters that initiate within the gene body (Maunakea et al., 2010), and the methylation of exons regulates their incorporation during splicing in a mechanism that involves the opposing effects of CTCF and MeCP2 on RNA polymerase II (Maunakea et al., 2013;Shukla et al., 2011). This raises the possibility that demethylation of gene bodies following fertilization is involved in the transition from gamete-specific promoters and exons to those employed in the early embryo.
Beyond its impact on overall expression levels, we uncovered a connection between intragenic methylation in the paternal pronucleus and transcriptional variability that endures throughout the preimplantation embryo. In contrast to the negative correlation reported in some somatic cells (Huh et al., 2013), intermediate methylation in the paternal pronucleus was consistently associated with higher noise levels relative to hypomethylated and hypermethylated gene bodies. Paternal demethylation at fertilization results in a massive conversion of hypermethylation to intermediate methylation levels at gene bodies, and is therefore linked to elevated transcriptional noise. Whether this relationship is causal will need to be examined experimentally, which will involve challenges similar to those faced in attempting to understand the role of demethylation in ZGA. Intermediate methylation levels in postreplicative zygotes could be the result of averaging differentially methylated strands or chromatids generated by processes such as replicative dilution and uneven active demethylation. Such differences would create a variation in methylation between cells upon subsequent divisions that could drive transcriptional noise. The use of strand-specific and single-cell methylation profiling technologies will provide insight into this possibility.
It has been posited that methylation heterogeneity provides a means of breaking symmetry during fate allocation through its influence on transcriptional activity , and indeed the loss of methylation greatly increases the symmetry of embryonic stem cell divisions (Jasnos et al., 2013). The association of zygotic demethylation with transcriptional noise invites speculation that this serves to generate heterogeneity in the early embryo that contributes to the first lineage decisions. We look forward to further investigations into the relationship among epigenetic reprogramming, heterogeneity, and lineage specification.
Dissecting the Contribution of TET3 to the Zygotic Methylation Landscape Oxidation by TET3 has been shown to provide an important pathway for removal of 5mC from the paternal genome at fertilization (Gu et al., 2011;Wossidlo et al., 2011), adding to other identified active pathways such as BER (Hajkova et al., 2010; and the elongator complex (Okada et al., 2010), and the potential for passive demethylation at the first DNA replication. However, the targets and interplay of these pathways have remained almost entirely uncharacterized. By combining our whole-genome methylation profiling with maternal deletion of zygotic TET3 activity, we were able to examine the contribution of TET3 to paternal methylation trajectories across the genome. This broadens findings from recently published RRBS-seq studies using similar genetic deletion approaches (Guo et al., 2014;Shen et al., 2014) and represents an important step toward mapping the division of labor among components of the demethylation machinery.
Although loss of TET3 activity disrupts methylation patterns at all single-copy features, the major impact is seen at intergenic, gene body, and non-CGI promoter sequences, where demethylation is impaired at more than a quarter of loci. As these regions constitute the bulk of demethylated loci, this establishes oxidation as a prevalent pathway for removal of 5mC across the genome. However, these data imply that TET3 is not required for demethylation of most loci, particularly among CGI-associated sequences, and this is also true of repetitive elements, which are largely unaffected in deletion zygotes. TET3 activity may be restricted to this limited number of loci; alternatively, it is possible that TET3 in fact targets additional loci, but that other demethylation pathways operate redundantly or are recruited in compensation and thereby neutralize the effect of its absence. Discriminating among these possibilities will require the simultaneous disruption of multiple demethylation pathways or the ability to generate detailed maps of TET3 binding to paternal pronuclear DNA. However, the observation that loss of TET3 frequently results in partial impairment of demethylation already demonstrates redundancy or compensation at a large number of loci. This is an important insight into the architecture of the demethylation machinery. Such redundancy could serve to limit the impact of the loss of a single pathway on methylation reprogramming and may be especially prevalent at sequences where demethylation must be ensured for developmental integrity. This could include LINE1 and MERV-L retrotransposons, whose activation is critical for early cleavage divisions (Beraldi et al., 2006;Kigami et al., 2003).
In addition to its contribution to zygotic demethylation, we uncovered a role for TET3 in protecting against the accumulation of aberrant de novo methylation at CGI-associated sequences. Mechanistically, it seems plausible that the same catalytic activity achieves both functions, i.e., TET3 provides protection through the oxidation of any deposited methyl group-a mecha-nism for active maintenance of a hypomethylated state in the zygote. Alternatively, TET3 could prevent accumulation of methylation independently of its oxidase function, perhaps by inhibiting the de novo methylation machinery or recruiting other chromatin modifiers. The elevated expression in 2-cell embryos of both demethylation and protection targets of TET3 suggests that these roles serve a common function in facilitating transcriptional activation.
Because bisulfite treatment cannot discriminate between 5mC and 5hmC (Huang et al., 2010), it is important to consider that the extent of demethylation may be underestimated when assessed with this technique. Furthermore, conventional bisulfite analysis is not informative regarding events downstream of the initial 5mC oxidation, as 5fC, 5caC, and unmodified cytosine behave alike . Previous work has shown that oxidized bases may be diluted passively due to a lack of maintenance at DNA synthesis (Hashimoto et al., 2012), and 5fC and 5caC can also be actively processed to unmodified cytosine by excision and BER, or decarboxylation Ito et al., 2011;Schiesser et al., 2012Schiesser et al., , 2013. Although some studies have used immunofluorescence or locus-specific modified bisulfite techniques to examine the use of these demethylation pathways (Guo et al., 2014;, it will be important for future studies to employ quantitative and genome-wide techniques that are capable of discriminating each cytosine modification in order to assess the full extent of 5mC removal and track the fate of oxidation products across the genomes of zygotes and early embryos.

Mouse Lines and Zygote Collection
All experimental procedures were approved by the Animal Welfare, Experimentation and Ethics Committee at the Babraham Institute and were performed under licenses by the Home Office (UK) in accordance with the Animals (Scientific Procedures) Act 1986. Mice with a conditional deletion in the Tet3 gene were generated for a previous study details in Supplemental Experimental Procedures).
Female mice were superovulated by intraperitoneal injection of pregnant mare's serum followed by intraperitoneal injection of human chorionic gonadotrophin (hCG) 48 hr later. Superovulated females were naturally mated with a stud from a 129S2/SvHsd or C57BL/6J background and examined for a vaginal plug the following day. Zygotes were collected 12-13 hr after the presumed time of insemination (i.e., 24-25 hr after hCG injection), at which stage zygotes have completed S phase (Santos et al., 2005). Verification of the pronuclear stage in batches of zygotes was performed manually.

Preparation of WGBS-Seq Libraries
Three independent collections of zygotes were obtained for each genotype (control and TET3; see Supplemental Experimental Procedures for details), giving a total of 225 control and 237 TET3 zygotes, of which 120 and 129 were derived from 129S2/SvHsd studs for control and TET3 samples, respectively. Zygotes were pooled and whole-genome bisulfite libraries were prepared using a post-bisulfite adaptor tagging (PBAT) strategy optimized from . An outline of this approach and a detailed protocol are provided in Supplemental Experimental Procedures.

Analysis of WGBS-Seq Data Sets
Details regarding analysis of the WGBS-seq data sets are provided in Supplemental Experimental Procedures.
Bisulfite data for C57BL/6J MII oocytes and sperm were drawn from Smallwood et al. (2014) (bulk and all single cells merged; Gene Expression Omnibus accession number GSE56879) and  (DNA Data Bank of Japan accession number DRA000484), respectively.
For data sets generated in this study, raw reads were trimmed using Trim Galore, aligned, and deduplicated, and methylation calls were extracted with Bismark (Krueger and Andrews, 2011) using custom pipelines for nonallele-specific analysis (alignment to NCBIM37) and paternal-specific analysis (alignment to both NCBIM37 and a 129S1/SvImJ genome).
For defined genomic features, probes that did not meet the threshold coverage for reliability were excluded. The methylation level was expressed as the mean of all sufficiently covered individual CG sites within the defined region. Methylation dynamics were calculated by comparing data sets using the SeqMonk chi-square filter with a p value requirement of <0.05 after multiple-test correction, and a minimum of ten observations and 10% difference.
For repetitive elements, Bismark was used to map all reads from each data set against consensus sequences constructed from Repbase (Jurka et al., 2005). The methylation level was expressed as the mean of individual CG sites.
Functional enrichment was assessed by analyzing gene lists with the DAVID web tool (Huang et al., 2009) using a background of all mouse genes. This was performed with the Functional Annotation Table feature for level 5 terms in the GOTERM databases. Terms with p < 0.05 after Benjamini correction were considered to be significant.

ACCESSION NUMBERS
The Gene Expression Omnibus accession number for the WGBS-seq data reported here is GSE63417.  (A) Methylation levels of imprint control regions (ICRs) in datasets used in this study. Each dot represents a single ICR, the horizontal bar denotes the mean. Taking the polar body into consideration, the predicted methylation value in post-replicative zygotes is 60% for maternal ICRs (assuming 100% and 0% methylation of maternal and paternal alleles respectively), and 40% for paternal ICRs (assuming 0% and 100% methylation of maternal and paternal alleles respectively). In the paternal component of the zygote, ICR methylation should be similar to values in sperm.

Supplemental Information includes Supplemental Experimental
Insufficient data were available to analyse paternal-specific values for paternal ICRs.
(B) Distribution of cytosines in CG context methylated at least 60% in sperm, amongst single-copy genomic features.
Data was drawn from WGBS-seq by ; only cytosines covered by at least five reads were included.  bootstrap resampling (n=10,000) using the pvclust package for R ; colours indicate to which feature group these refer.
(C) Expression of TET3-demethylated loci in 2-cell embryos (data drawn from . The box plot shows the distribution of log (base 2) transformed counts of reads per million total reads (RPM), displaying the median (bar), interquartile range (box), maximum and minimum (whiskers). The p-values shown are the result of an unpaired ANOVA multiple comparison test with Sidak correction. 'All genes' refers to all annotated genes according to Ensembl. TET3-demethylated loci are defined as demethylated loci (according to criteria in Figure S1E) that underwent an absolute paternal methylation gain of greater than 10% in the absence of TET3 and were significant according to a Chi-Squared test (corrected p-value <0.05). N=1310 intergenic sequences, 2487 gene bodies sequences, 259 CGI promoters and 398 non-CGI promoters.   Mice with a conditional deletion in the Tet3 gene were generated for a previous study . Briefly, C57BL/6NTac ES cells (TaconicArtemis) were targeted with a vector introducing LoxP sites around exon five of Tet3 RefSeq NM_183138.2 (sequence: CCGGACCTGTGCTTGCCAAGGCAAAGACCCTAACACCTGCGGTGCCTCCTTCTCCTTCGGCTGTTCCTGG AGCATGTACTTCAACGGCTGCAAATATGCTCGGAGCAAGACGCCACGAAAGTTCCGCCTCACGGGAGAC AATCCGAAGGAG) which encodes residues required for chelation of Fe(II) and is upstream of exons containing other key catalytic residues (Tahiliani et al., 2009). Expression of Cre recombinase results in excision of this region and a frameshift from exon six affecting all downstream exons and introducing a premature stop codon in exon seven.
Animals heterozygous for the floxed allele were bred to a transgenic mouse line containing Zp3-Cre on a C57BL/6J background (Lewandoski et al., 1997) and homozygous floxed mice were generated by inter-crossing to give either control females (flox/flox, Zp3-) or conditional deletion (referred to as 'TET3') females (flox/flox, Zp3+). Complete efficacy of Cre excision was confirmed in the oocytes of conditional deletion females by RT-PCR between exons four and six followed by Sanger sequencing. Our previous study confirmed that TET3 activity is abolished in zygotes derived from flox/flox Zp3+ females .

Preparation of Whole-Genome Bisulfite Libraries
Whole-genome bisulfite libraries were prepared using a post-bisulfite adaptor tagging (PBAT) strategy optimised from . In this approach, bisulfite treatment is used to both convert cytosines and fragment DNA at the start of the protocol. This combination eliminates the need for fragmentation by sonication, and avoids degrading adaptor-tagged sequences. As adaptors cannot be ligated to bisulfite-treated DNA in the same manner as standard protocols, they are instead bound via a random tetramer sequence in two stages.

Bioinformatic Analysis of Whole-Genome Bisulfite Sequencing
For the datasets generated in this study, raw reads were trimmed to remove both poor quality calls and adapters using Trim Galore (v0.3.5, default parameters; www.bioinformatics.babraham.ac.uk/projects/trim_galore). In order to overcome hybridisation artefacts, trimmed reads were first aligned to the mouse genome (build NCBIM37) in pairedend mode (to be able to count overlapping parts of the reads only once) while writing out unmapped singleton reads; in a second step remaining singleton reads were aligned in single-end mode. Alignments were carried out with Bismark v0.10.0  with the following set of parameters: a) paired-end mode: --pbat; b) single-end mode for Read 1: --pbat; c) single-end mode for Read 2: defaults. Reads were then deduplicated with deduplicate_bismark selecting a random alignment for positions that were covered more than once. CG methylation calls were extracted from the deduplicated mapping output ignoring the first 4bp to reduce the methylation bias typically observed in PBAT libraries using the Bismark methylation extractor (v0.10.0) with the following parameters: a) paired-end mode: --no_overlap --report --ignore 4 --ignore_r2 4; b) single-end mode: --report --ignore 4.
To generate paternal-specific profiles for zygotes, alignments were carried out against both NCBIM37 as well as a 129S1/SvImJ genome which was constructed by injecting 4,334,345 high-confidence flagged SNPs from the Mouse Genomes Project (http://www.sanger.ac.uk/resources/mouse/genomes). SNPs close to indels, repeat element or those having two possible base calls for each genome were excluded. First, paired-end alignments were run and then sorted allele-specifically depending on the number of non-bisulfite mismatches to the two genomes. To allow the allelespecific sorting with a custom Perl script, Bismark was run with the option '--non_bisulfite_mismatches'. Sequences that aligned equally well to both genomes were classified as Unassigned alignments, others were partitioned into either C57BL/6J (NCBIM37) or 129S1/SvImJ genome specific alignments. Second, unaligned singleton reads were mapped against both genomes in single-end mode and any reads that were already aligned during the paired-end step were removed from the single-end alignment output. Next, reads 1 and 2 were also sorted allele-specifically. Once sorting had completed, all alignment outputs were deduplicated, methylation calls extracted and the resulting methylation call files were merged for further analysis. Analysis of mapped data was performed using SeqMonk (www.bioinformatics.babraham.ac.uk/projects/seqmonk) with the NCBI37 genome build.
CG methylation calls from uniquely mapped reads were analysed according to probes designed across features as follows. Specific genomic features were defined as follows. CGIs were defined by pull-down data (Illingworth et al., 2010) and were subdivided into promoter and non-promoter CGIs based on whether they overlapped with a promoter region (defined below). Promoters were defined as the region 1.5 kb upstream and 0.5 kb downstream from a TSS (annotation from Ensembl [Flicek et al., 2014]), and were subdivided into CGI and non-CGI promoters based on whether they overlapped with a CGI (defined above). For gene bodies, 5 kb probes were designed across the entire length of a gene and for intergenic sequences, 5 kb probes were designed over the entire genome excluding genes and regions 2 kb either side.
The coverage thresholds that were applied across all datasets for inclusion of features are as follows: CGIs, two reads at a minimum of four cytosines; promoters, two reads at a minimum of five cytosines; gene bodies, two reads at a minimum of 10 cytosines; intergenic sequences, two reads at a minimum of 10 cytosines.
Loci were defined as associated with a gene according to the following criteria: gene bodies, overlapping; promoters, overlapping; CGIs, within 20kb.

RNA-seq Dataset Analysis
For analysis of RNA-seq data generated by Park et al., 2013 (accession number DRA001066), raw data were mapped to the mouse genome assembly (build GRCm38) using TopHat (v1.4.1, options -g 1) in conjunction with gene models from Ensembl release 61. The number of reads mapping to each annotated mRNA transcript was quantified then corrected for total read count (reads per million, RPM) without length correction. The log 2 RPM value of genes associated with indicated loci (defined as for gene ontology analysis) was used to perform statistical tests and generate box-whisker plots using Prism 6 (GraphPad).

Analysis of Transcriptional Noise
We use the term 'transcriptional noise' to refer to the degree of variation in a gene's transcription level between individual cells. A metric of this noise was calculated for the paternal allele of each expressed gene at each stage of early development as follows. Allele-specific expression measurements were drawn from  and normalised for total read count (RPM). For all cells at each developmental stage, the mean expression and squared coefficient of variation (CV 2 ) for each gene were computed and both log10 transformed. Given the negative correlation between the CV 2 and a gene's expression level (Brennecke et al., 2013), we calculated a mean corrected value of noise by subtracting from each gene's empirical CV 2 the median CV 2 for genes expressed at that level computed using a running median (Newman et al., 2006). Genes not expressed at a particular stage were excluded from the analysis of that stage. Noise levels were compared to average intragenic methylation calculated from the control paternal pronucleus dataset using probes designed across the entire gene body; only genes meeting the threshold of two reads at a minimum of 10 cytosines were included in the analysis. The number of genes analysed in the 0-25%, 25-50%, 50-75% and 75-100% quartiles, respectively, at each stage is as follows; Late 2-cell: 445,795,686 and 70;808,684 and 65;877,783 and 83;865,759 and 78;Blastocyst: 509,972,892 and 87.