Research ArticleGENETICS

Allelic H3K27me3 to allelic DNA methylation switch maintains noncanonical imprinting in extraembryonic cells

See allHide authors and affiliations

Science Advances  20 Dec 2019:
Vol. 5, no. 12, eaay7246
DOI: 10.1126/sciadv.aay7246

Abstract

Faithful maintenance of genomic imprinting is essential for mammalian development. While germline DNA methylation–dependent (canonical) imprinting is relatively stable during development, the recently found oocyte-derived H3K27me3-mediated noncanonical imprinting is mostly transient in early embryos, with some genes important for placental development maintaining imprinted expression in the extraembryonic lineage. How these noncanonical imprinted genes maintain their extraembryonic-specific imprinting is unknown. Here, we report that maintenance of noncanonical imprinting requires maternal allele–specific de novo DNA methylation [i.e., somatic differentially methylated regions (DMRs)] at implantation. The somatic DMRs are located at the gene promoters, with paternal allele–specific H3K4me3 established during preimplantation development. Genetic manipulation revealed that both maternal EED and zygotic DNMT3A/3B are required for establishing somatic DMRs and maintaining noncanonical imprinting. Thus, our study not only reveals the mechanism underlying noncanonical imprinting maintenance but also sheds light on how histone modifications in oocytes may shape somatic DMRs in postimplantation embryos.

INTRODUCTION

Genomic imprinting in mammals refers to a series of precisely regulated epigenetic processes that lead to parental allele–specific gene expression (1). Genomic imprinting plays critical roles in fetal development, placental function, homeostasis, and nervous system development (2). Since the discovery of genomic imprinting, allele-specific DNA methylation at the imprinting control regions (ICRs) has been the primary epigenetic modification that governs monoallelic gene expression within an imprinted cluster (1). The differential DNA methylation at the ICRs is established during gametogenesis [i.e., germline or primary differentially methylated regions (DMRs)], maintained in preimplantation embryos when global epigenetic reprogramming takes place, and eventually erased in primordial germ cells to reset the sperm- or oocyte-specific imprints (1). In addition to the primary DMRs, additional epigenetic mechanisms such as somatic DMRs (also known as secondary DMRs) and allele-specific histone modifications established in later embryonic development are also involved in imprinted gene expression (3, 4).

Although most imprinted loci harbor germline DMRs and show loss of imprinting in the absence of germline DNA methylation (5, 6), a few paternally expressed genes (PEGs) can maintain imprinted expression in the absence of oocyte DNA methylation (7, 8). This maternal DNA methylation–independent imprinted expression indicated the existence of a noncanonical imprinting mechanism. We recently demonstrated that oocyte-specific H3K27me3 serves as the primary imprinting mark to repress the maternal transcription of several dozens of PEGs in early embryos, including all imprinted genes that were previously known to be independent of oocyte DNA methylation (9, 10). Maternal H3K27me3 also regulates imprinted X inactivation by repressing the maternal transcription of Xist in preimplantation embryos (1012). Notably, maternal allele–biased H3K27me3 correlates with paternal-specific gene expression in human morulae (13). However, contrary to mice in which oocyte-inherited H3K27me3 is largely maintained in preimplantation development (14), H3K27me3 undergoes global depletion at four- to eight-cell stage in human early embryos (15). Therefore, the previously observed maternal allele–biased H3K27me3 might be reestablished in morulae rather than directly inherited from oocytes.

Canonical and noncanonical imprintings are two distinct epigenetic mechanisms. These two types of imprinting are established independently in oocytes, as imprinted expression of noncanonical imprinting, but not canonical imprinting, is lost in EED maternal knockouts (matKOs) (10). Furthermore, while canonical imprinting is relatively stable in development, noncanonical imprinting is largely transient in early embryos, with some genes maintaining imprinted expression in extraembryonic lineage (9). In this study, by analyzing several mouse models [i.e., DNMT3L matKO, EED matKO, and zygotic DNMT3A/3B double knockout (DKO)], we dissected the mechanism of noncanonical imprinting maintenance. Through integrative analyses of DNA methylome, transcriptome, and H3K27me3/H3K4me3 profiles, we demonstrate that oocyte H3K27me3 is required for establishing allele-specific de novo DNA methylation at implantation (somatic DMRs) and the somatic DMRs are essential for maintaining noncanonical imprinting in extraembryonic cells. Therefore, we conclude that although noncanonical imprinting initially relies on maternal H3K27me3 to regulate imprinted gene expression in preimplantation embryos, it eventually becomes DNA methylation dependent in postimplantation development.

RESULTS

Noncanonical imprinting is independent of oocyte DNA methylation

We recently reported that oocyte-derived H3K27me3 can silence the maternal allele for imprinted expression of dozens of genes during mouse preimplantation development (9, 11). However, only few of them have been demonstrated to be oocyte DNA methylation independent previously (8). To further confirm that noncanonical imprinting is independent of oocyte DNA methylation, we performed allele-specific gene expression analyses in DNMT3L matKO embryos by RNA sequencing (RNA-seq) (Fig. 1A). DNMT3L is a noncatalytic regulatory factor of de novo DNA methyltransferases and is essential for establishing oocyte DNA methylation (6). To differentiate parental origins of the transcripts, oocytes from DNMT3L KO (i.e., DNMT3LG/G) female mice (B6/129 mixed background) (6) were fertilized with wild-type (WT) PWK strain sperm and morulae [~78 hours postfertilization (hpf)] were collected for RNA-seq analyses. After confirming the RNA-seq data reproducibility (fig. S1A), we performed comparative gene expression analyses, which revealed that the absence of oocyte DNA methylation only caused 77 (0.74%) and 128 (1.23%) genes, respectively, up- and down-regulated in DNMT3L matKO morulae [RPKM (reads per kilobase million) > 1, fold change (FC) > 2, and false discovery rate (FDR) < 0.05] (fig. S1B and table S1). This limited gene expression alteration is in line with the fact that these embryos can survive to midgestation, when lethality takes place because of loss of oocyte DNA methylation–dependent imprinting (6).

Fig. 1 Noncanonical imprinting is independent of oocyte DNA methylation.

(A) Schematic diagram of the experimental strategy for assessing allele-specific gene expression in morulae. (+) Dnmt3l WT allele; (−) Dnmt3l KO allele. Dashed boxes indicate the samples collected for RNA-seq. (B) Heat maps showing the allelic expression bias of oocyte DNA methylation–dependent and maternal H3K27me3–dependent PEGs in morulae. Two independent samples were used for RNA-seq for each genotype. Genes with >20 SNP reads in both replicates are shown. (C and D) Representative genome browser views of an oocyte DNA methylation–dependent PEG Snrpn (C) and a maternal H3K27me3–dependent PEG Gab1 (D). For DNA methylation tracks, gray bars indicate the CpG sites with >5 reads coverage. Published H3K27me3 ChIP-seq data were from (14, 23), and DNA methylation data were from (44, 45).

Consistent with the known role of DNMT3L in regulating maternal DNA methylation imprints (6), the paternal allele expression bias of all oocyte DNA methylation–dependent PEGs [>20 single-nucleotide polymorphism (SNP)–containing reads] were either lost or became milder in DNMT3L matKO embryos (Fig. 1, B and C; fig. S1C; and table S1). In contrast, all the maternal H3K27me3–dependent PEGs (>20 SNP-containing reads), which were demonstrated to be biallelically expressed in EED maternal KO embryos (10), largely maintained their paternal allele–biased expression in DNMT3L matKO morulae (Fig. 1, B and D; fig. S1, D and E; and table S1). Although the reciprocal cross (i.e., PWK × B6/129) is not feasible under this experimental setting as the DNMT3L KO mice were maintained in a B6/129 background, we have previously shown that these maternal H3K27me3–dependent PEGs were paternally expressed at morula stage in both B6 × CAST and CAST × B6 crosses (9). Notably, although the Slc38a4 promoter overlaps with a germline differential DNA methylated region (DMR), its paternal allele–specific expression is unaffected in DNMT3L matKO embryos (Fig. 1B and fig. S1E), indicating that differential DNA methylation is not a key determinant for its allelic expression at this stage. This observation is in agreement with the previous reports that imprinted expression of Slc38a4 is mainly regulated by maternal H3K27me3 (8, 10). Collectively, these data demonstrate that the lack of DNMT3L only causes loss of canonical imprinting, but not noncanonical imprinting.

Maternal H3K27me3 is required for establishing allele-specific DNA methylation at implantation

The fact that Gab1, one of the maternal H3K27me3–dependent PEGs, harbors a somatic DMR at implantation and is imprinted in placenta (7, 8) prompted us to postulate that acquisition of somatic DMR might be the underlying mechanism for maintaining noncanonical imprinting. To test this possibility, we performed whole-genome bisulfite sequencing (WGBS) analyses in extraembryonic ectoderm (ExE) of embryonic day 6.5 (E6.5) embryos, cells known to maintain noncanonical imprinting (9). To assess whether the potential somatic DMRs are dependent on maternal H3K27me3, we also analyzed ExE of EED matKO embryos obtained by crossing WT PWK strain male mice and Gdf9Cre, Eedfl/fl female mice (B6/129 mixed background) whose oocyte H3K27me3 is depleted (10). The EED control (Eedfl/fl × PWK, CTR) WGBS dataset is highly correlated with a public dataset (fig. S2, A and B) (16). Furthermore, allelic DNA methylation of known ICRs was largely maintained in EED matKO ExE with the exception of the Zrsr1/Commond1 locus (fig. S2C and table S2). These data validate the correct separation of parental-specific WGBS reads and are consistent with the previous observation that canonical imprinting is generally not regulated by maternal H3K27me3 (10).

We noticed that three (Gab1, Phf17, and Smoc1) of the five oocyte H3K27me3–dependent PEGs that maintain imprinting in ExE (9) acquire maternal allele–specific DNA methylation in CTR at implantation (Fig. 2A). Of the three somatic DMRs, Gab1 DMR and Phf17 DMR are covered by allele-specific reduced representation bisulfite sequencing (RRBS) reads, which also revealed maternal allele–biased DNA methylation (see later). These somatic DMRs are lost in EED matKO ExE, indicating that establishment of these allelic DMRs depends on oocyte H3K27me3 (Fig. 2A). Notably, the germline DMR of Slc38a4 is hypomethylated in EED matKO ExE (Fig. 2B), suggesting that maintenance of Slc38a4 DMR also requires oocyte-derived H3K27me3. No somatic DMR was identified for Sfmbt2 in E6.5 ExE, suggesting that imprinted expression of this gene is unlikely to be regulated by DNA methylation at this developmental stage.

Fig. 2 Oocyte H3K27me3 controls somatic DMR establishment.

(A) Genome browser views of DNA methylation levels at oocyte H3K27me3–dependent somatic DMRs. Areas shaded by purple indicate the identified somatic DMRs. DNA methylation data of E3.5 trophectoderm (TE) were obtained from a public dataset (16). (B) Genome browser view of DNA methylation levels at the Slc38a4 locus. No SNPs are available to differentiate the allele-specific reads at Slc38a4 promoter. (C and D) Heat maps showing the allelic expression bias of identified noncanonical imprinting loci in B6 and PWK reciprocal crosses in E6.5 ExE (C) and EED CTR and matKO E6.5 ExE (D). The E6.5 ExE RNA-seq data were obtained from previous studies (9, 10). B, B6; P, PWK. The names of the newly identified noncanonical imprinting loci are bolded. (E) Genome browser views of the two newly identified noncanonical imprinting loci. Shaded area (black bars) indicates the associated somatic DMRs. Oocyte H3K27me3 and DNA methylation data were obtained from previous studies (23, 44).

We next sought to identify additional somatic DMRs whose establishment depends on oocyte H3K27me3. We reasoned that the candidate DMRs should be hypomethylated in EED matKO when compared to CTR and are associated with maternal H3K27me3–dependent PEGs. In addition, the candidate DMRs should also show maternal allele–specific DNA methylation in ExE if parental-specific SNPs are available. Assessing genome-wide DNA methylation in 100–base pair (bp) tiles (>2 CpGs and >10 reads coverage in both samples) revealed that global DNA methylation is mildly decreased in EED matKO ExE (mean: 0.58 versus 0.51) (fig. S2D). Consistently, of the 2412 DMRs identified between the two groups (>10 CpGs and >25% methylation differences), 2333 showed lower DNA methylation (hypo-DMRs) in EED matKO ExE compared with CTR, while only 79 were identified as hyper-DMRs (fig. S2E and table S2). Following the procedure described (fig. S2F), we identified six genes/regions that include four (Gab1, Phf17, Smoc1, and Slc38a4) known maternal H3K27me3–dependent PEGs (Fig. 2, C and D), which further validates our strategy. The rest of hypo-DMRs that were not associated with any imprinted genes are likely contributed by maternal H3K27me3 depletion through imprinting-independent mechanisms. Of note, two genes, Platr20 and Gm32885, were identified as previously unknown maternal H3K27me3–dependent PEGs. These two genes overlap with H3K27me3 domains, but not DNA hypermethylated regions, in oocytes and exhibit loss (or become weaker) of biased allelic expression in EED matKO ExE samples (Fig. 2, C to E). Together, these analyses allowed us to identify two new noncanonical imprinted genes and demonstrate that maternal H3K27me3 is required for establishing somatic DMRs at implantation (fig. S3A).

Somatic DMRs are essential for maintaining noncanonical imprinting

We next sought to determine whether the somatic DMRs are required for maintaining noncanonical imprinting in ExE. Analyses of public datasets showed that removal of DNMT3A, DNMT3B, or DNMT3L individually does not completely abolish the establishment of the somatic DMRs at the noncanonical imprinting loci (fig. S3B) (16). Thus, preventing establishment of somatic DMRs at noncanonical imprinted genes may require simultaneous disruption of both de novo DNA methyltransferases DNMT3A and DNMT3B (17).

To this end, we attempted to deplete both DNMT3A and DNMT3B by zygotic CRISPR-Cas9 injection, followed by isolating ExE (Fig. 3A). Disruption of Dnmt3a and Dnmt3b was first confirmed by reverse transcription quantitative polymerase chain reaction (RT-qPCR) and then by RNA-seq analyses (Fig. 3B and fig. S4A). RNA-seq data quality was validated by the high reproducibility between replicates, and the purity of the ExE cells was verified by cell lineage–specific marker expression analysis (fig. S4, B and C). Low-input RRBS analysis revealed that the global DNA methylation level of DNMT3A/3B DKO ExEs is lower than that of DNMT3A or DNMT3B single-KO ExE (16) and is as low as that of E3.5 trophectoderm (TE) (Fig. 3C). Previous studies have established that global de novo methylation has not yet taken place at E3.5 (16, 18). We also confirmed that regions that normally acquire DNA methylation at implantation (100-bp tiles: <25% in E3.5 TE and >75% in E6.5 ExE) remain hypomethylated in the DNMT3A/3B DKO ExE (fig. S4D). Collectively, these results demonstrate that de novo DNA methylation in ExE is largely abolished by depletion of zygotic DNMT3A/3B.

Fig. 3 Somatic DMRs are essential for maintaining noncanonical imprinting.

(A) Schematic diagram of the experimental strategy. Cas9 mRNA and six single-guide RNAs (sgRNAs) (three each targeting Dnmt3a or Dnmt3b) were injected into BDF1/PWK zygotes. Zygotes injected with water were used as the WT control. (B) RT-qPCR results showing relative expression levels of Dnmt3a and Dnmt3b in WT and DKO ExE. The expression levels in WT embryo #1 were set as 1.0. Complementary DNA (cDNA) was synthesized using a small piece of ExE as illustrated in (A). After verification of Dnmt3a/3b depletion by RT-qPCR, cDNA was then used to prepare RNA-seq libraries for next-generation sequencing. (C) Global view of DNA methylation status in TE and ExE of the indicated genotypes. Percentages of 100-bp tiles at four methylation levels are shown. Common CpGs between E3.5 TE WGBS data and E6.5 ExE RRBS data were retrieved, and tiles with >2 CpGs and >10 reads coverage in all samples were used to generate this plot (n = 317,106). For both WT and DKO, three low-input RRBS replicates were pooled and analyzed. Embryos #4, #6, and #7 shown in (B) were used for the RRBS experiments. WGBS data of E3.5 TE and RRBS data of DNMT3A and DNMT3B single-KO ExE were obtained from a previous study (16). (D) Genome browser views of DNA methylation levels covering the somatic DMRs of three maternal H3K27me3–dependent imprinted genes. (E) Genome browser views of DNA methylation levels at representative ICRs of oocyte DNA methylation–dependent imprinted loci KvDMR1 and Peg3. (F) Scatter plot comparing the gene expression levels between WT and DNMT3A/3B DKO E6.5 ExE. Two (WT) and seven (DKO) RNA-seq replicates were used for differential gene expression analyses. (G) Heat map showing the allelic expression bias of canonical and noncanonical imprinting loci. Genes with >20 SNP-containing RNA-seq reads are shown.

Next, we asked whether the somatic DMRs of the noncanonical-imprinted genes are disrupted by zygotic DNMT3A/3B depletion. We found that, of the five genes that acquire somatic DMRs at implantation (Fig. 2), Gab1, Phf17, and Platr20 promoters exhibited maternal allele–specific DNA methylation in WT, but not in the DNMT3A/3B DKO ExE (Fig. 3D). Although allelic DNA methylation analyses at the promoters of Gm32885 and Smoc1 are not feasible because of the lack of parental-specific SNPs and/or allelic RRBS reads coverage, the overall methylation levels at these loci were reduced in the DKO samples (fig. S4E). The Slc38a4 germline DMR also exhibited reduced DNA methylation in DNMT3A/3B DKO ExE (fig. S4E), which is consistent with the modest decrease in DNA methylation at this locus in DNMT3B single-KO ExE (fig. S3B) (16). In contrast to noncanonical imprinting, canonical imprints were largely unaffected by zygotic depletion of DNMT3A/3B (Fig. 3E, fig. S4F, and table S3). This observation supports the notion that germline DNA methylation imprints generally do not undergo demethylation or remethylation after fertilization. Collectively, these data demonstrate that DNMT3A/3B removal abolishes the establishment of somatic DMRs in ExE at the noncanonical imprinting genes.

We next asked whether disruption of the somatic DMR establishment affects the allelic expression of the noncanonical imprinted genes in the ExE cells. Comparative gene expression analyses revealed that DNMT3A/3B DKO only modestly affected E6.5 ExE [209 (2.04%) up-regulated and 90 (0.88%) down-regulated genes, RPKM > 1, FC > 2, and FDR < 0.05] (Fig. 3F and table S3), which is in line with the fact that DNMT3A/3B DKO embryos can undergo normal gastrulation (17). With the exception of Sfmbt2, the overall transcript abundance of maternal H3K27me3–dependent PEGs increased by about twofold in DKO when compared to WT (Fig. 3F). Furthermore, the paternal allele–biased expression of these genes, except Sfmbt2, was either lost or became milder in DNMT3A/3B DKO ExE (Fig. 3G and table S3), indicating that zygotic DNMT3A/3B is required for maintaining noncanonical imprinting. In contrast, ablation of zygotic DNMT3A/3B only mildly affected the allelic biased expression of germline DNA methylation–dependent imprinted genes (Fig. 3G), which is consistent with the largely unaffected allelic DNA methylation at their ICRs (Fig. 3E, fig. S4F, and table S3). The variable loss of imprinting of the maternal H3K27me3–dependent PEGs could be due to genetic mosaicisms caused by CRISPR-Cas9–mediated acute depletion (19). In addition, as these genes are expected to be paternally expressed in DNMT3A/3B DKO preimplantation embryos, the incomplete degradation of the transcripts might lead to an overestimation of the paternal allele–biased biased expression. Nonetheless, these results demonstrate that DNMT3A/3B-dependent somatic DMRs are essential for maintaining noncanonical imprinting in E6.5 ExE.

Somatic DMRs take over maternal H3K27me3 to maintain imprinting in E6.5 ExE

Given that the somatic DMRs are essential for maintaining noncanonical imprinting in postimplantation embryos, we wondered whether maternal H3K27me3 still has a role in regulating allelic expression of noncanonical imprinted genes at this developmental stage. Because oocyte H3K27me3 domains start to diminish after fertilization and only show weak enrichment at blastocyst stage (9), we hypothesize that oocyte-derived H3K27me3 has lost by E6.5 and the somatic DMRs take over to maintain noncanonical imprinting. To test this hypothesis, we profiled H3K27me3 in single E6.5 ExE using Cleavage Under Targets and Release Using Nuclease (CUT&RUN) assay (20). To probe the potential relationship between de novo DNA methylation and H3K27me3 in postimplantation embryos, H3K27me3 CUT&RUN was also performed in DNMT3A/3B DKO ExE samples.

The H3K27me3 CUT&RUN replicates showed high reproducibility, and H3K27me3 enrichment at canonical Polycomb targets, such as Hox clusters, was not affected in the absence of de novo DNA methylation (fig. S5, A and B). Consistent with the modest gene expression changes in DNMT3A/3B DKO ExE (Fig. 3F), global H3K27me3 level and distribution are comparable between WT and DKO samples (fig. S5, C and D). In addition, the germline DMR-dependent allele-specific H3K27me3 enrichment at canonical imprinting loci was either unaffected or mildly affected by DNMT3A/3B DKO (fig. S5E), which is consistent with the largely maintained imprinted expression at these loci in the DKO samples (Fig. 3G).

We next analyzed the maternal H3K27me3 domains in ExE and found that these domains, which are evident in morulae, are completely lost at E6.5 (Fig. 4, A and B). Nevertheless, Phf17 and Sfmbt2 retained maternal allele–biased H3K27me3 peaks nearby their promoters (Fig. 4, C and D). The peaks at the Phf17 locus are lost in DNMT3A/3B DKO (Fig. 4C), indicating that this allelic enrichment of H3K27me3 depends on de novo DNA methylation. In contrast, maternal H3K27me3 domain at the Sfmbt2 locus was unaffected by zygotic DNMT3A/3B DKO (Fig. 4D), which was in line with our observation that the somatic DMR is not acquired until E6.5 and its imprinted expression is not altered by DNMT3A/3B DKO (Fig. 3, F and G). Thus, imprinted expression of Sfmbt2 may be still regulated by the maternal H3K27me3 domain at E6.5, which may then be taken over by a somatic DMR that was recently identified in E14.5 placenta (21). Together with the indispensable role of somatic DMRs in noncanonical imprinting regulation, these data suggest that noncanonical imprinting is regulated by maternal H3K27me3 in preimplantation embryos and then by somatic DMRs after implantation.

Fig. 4 Oocyte-inherited H3K27me3 domains are absent in E6.5 ExE.

(A) Average enrichment of H3K27me3 [log2(maternal/paternal)] at the maternal H3K27me3 domains in morulae, WT ExE, and DNMT3A/3B DKO ExE. Maternal H3K27me3 domains previously identified in BDF1 × PWK IVF morula embryos (n = 4134) (31) were used to generate the plot. Note that the genetic background of the morula sample is the same as the ExE samples. For both WT and DNMT3A/3B DKO, two H3K27me3 CUT&RUN replicates were pooled for the analyses. For the DKO group, embryos #1 and #5 (Fig. 3B) were used for H3K27me3 CUT&RUN experiments separately and pooled for analysis. (B to D) Genome browser views of allelic H3K27me3 enrichment at Gab1, Platr20, Slc38a4, Smoc1, Phf17, and Sfmbt2 loci. The arrows in (C) points to the maternal allele–biased H3K27me3 peaks that are observed in WT but not in the DNMT3A/3B DKO ExE.

Oocyte H3K27me3 regulates allelic H3K4me3 in preimplantation embryos

After demonstrating the essential role of somatic DMRs in noncanonical imprinting maintenance, we set out to understand how oocyte H3K27me3 instructs allele-specific de novo DNA methylation at implantation. Because the somatic DMRs overlap with promoters of the maternal H3K27me3–dependent PEGs and promoter-associated H3K4me3 can antagonize DNA methylation (22), we asked whether the genes that acquire somatic DMRs exhibit allelic H3K4me3 in preimplantation embryos. To this end, we profiled H3K4me3 using CUT&RUN assay (20) in EED CTR and matKO morula and TE. The CUT&RUN replicates of morulae showed high reproducibility and are well correlated with the publicly available dataset (fig. S6, A and B) (23). Furthermore, comparative analyses revealed that only a small percentage of H3K4me3 peaks were gained or lost in EED matKO morula (fig. S6, C and D), which is consistent with the modest gene expression changes in EED matKO morulae (10). Similar results were also obtained in TE (fig. S6, E to H).

Of the five maternal H3K27me3–dependent PEGs that establish somatic DMRs, Gab1, Phf17, Smoc1, and Platr20 have parental-specific SNPs for allelic analysis and exhibited paternal allele–biased H3K4me3 in both morulae and TE at regions that acquire somatic DMRs in E6.5 ExE (Fig. 5A). These regions gained H3K4me3 on the maternal allele to a level comparable to that of the paternal allele in EED matKO morulae and TE (Fig. 5A), indicating that the allelic-biased H3K4me3 distribution is dependent on oocyte H3K27me3. In addition, analyses of publicly available H3K4me3 chromatin immunoprecipitation sequencing (ChIP-seq) datasets of one-, two-, four-, and eight-cell embryos (24) revealed that the paternal allele–biased H3K4me3 enrichment was established around the four- to eight-cell stages (Fig. 5B). Given the known role of H3K4me3 in repelling DNA methyltransferases (22), these data suggest that oocyte-inherited H3K27me3 instructs somatic DMR establishment potentially by establishing allelic H3K4me3 in preimplantation embryos.

Fig. 5 Oocyte H3K27me3 regulates allelic H3K4me3 in preimplantation embryos.

(A) Genome browser views of DNA methylation and allelic H3K4me3 enrichment in EED CTR and matKO morulae and TE of blastocyst embryos. The shaded areas indicate the somatic DMRs. Morulae were collected at ~78 hpf, and TE samples were dissected from ~120 hpf blastocysts. (B) Genome browser views of allelic H3K4me3 levels at Gab1, Phf17, Platr20, and Smoc1 loci in each stage of preimplantation embryos. The H3K4me3 ChIP-seq data were obtained from a previous publication (24). (C) Heat map showing the allelic enrichment of promoter H3K4me3 peaks. Only noncanonical imprinting loci with sufficient (>50) SNP-containing reads of the associated promoter H3K4me3 peaks are shown. “+”: locus establishes somatic DMR in ExE or overlap with CGI; “−”: locus has no somatic DMR in ExE or does not overlap with CGI. (D) Genome browser views of DNA methylation and allelic H3K4me3 levels at G730013B05Rik, Tle3, and Grik3.

To investigate whether allelic H3K4me3 plays a role in selecting genes to acquire somatic DMRs and maintain imprinted expression in extraembryonic lineage, we asked to what extent the maternal H3K27me3–dependent PEGs correlate with allelic H3K4me3 in morulae and TE. Of the 78 putative noncanonical imprinted genes identified in preimplantation embryos (Fig. 2, C and D) (9), 52 (66.6%) were associated with promoter H3K4me3 peaks (see Materials and Methods). Of those, 17 (32.6%) had sufficient parental-specific SNPs (>50 SNP-containing reads) in both morulae and TE for allelic H3K4me3 enrichment analyses (Fig. 5C and table S4). Of the eight genes exhibiting oocyte H3K27me3–dependent paternal H3K4me3 bias in both morula and TE (Fig. 5C), four genes (i.e., Sfmbt2, G730013B05Rik, Tle3, and Grik3) do not acquire somatic DMR in E6.5 ExE (Fig. 5, C and D), indicating that allelic H3K4me3 in preimplantation embryos alone is not sufficient to instruct allelic DNA methylation at implantation. Notably, G730013B05Rik, Tle3, and Grik3 promoters overlap with CpG islands (CGIs) and remain hypomethylated in E6.5 ExE (Fig. 5, C and D). It is possible that the intrinsic property of CGIs to antagonize DNA methyltransferases (18) may prevent both alleles from gaining DNA methylation at implantation. In summary, although paternal allele–biased enrichment of H3K4me3 in TE coincides with some somatic DMRs at E6.5, additional mechanisms other than allelic H3K4me3 are likely involved in regulating somatic DMR establishment at implantation.

DISCUSSION

Germline DMRs have been the only primary epigenetic marks that can direct somatic DMRs (~20) establishment in postimplantation development (4, 21, 25). Here, we identified oocyte H3K27me3 as another primary epigenetic signal that contributes to the establishment of at least five somatic DMRs at implantation (i.e., Gab1, Phf17, Smoc1, Platr20, and Gm32885) (Fig. 2). These somatic DMRs are essential for maintaining noncanonical imprinting in extraembryonic lineage (Fig. 3). Together with the observation that the oocyte-inherited H3K27me3 domains are largely lost by E6.5 (Fig. 4), our data support the idea that the somatic DMRs take over the maternal H3K27me3 domains to maintain noncanonical imprinting in postimplantation development (Fig. 6).

Fig. 6 Model illustrating that maintenance of maternal H3K27me3–dependent imprinting requires somatic DMR establishment in postimplantation embryos.

Oocyte H3K27me3 regulates somatic DMR establishment

The mechanism underlying somatic DMR establishment is largely unknown. It was previously reported that the transient transcription of an RNA Liz in early embryos locally promotes a somatic DMR establishment (26). However, this mechanism so far only applies to the Gpr1/Zdbf2 imprinted locus (26). Here, we demonstrate that oocyte-specific H3K27me3 can regulate allelic H3K4me3 deposition in preimplantation embryos (Fig. 5), which, in turn, might regulate somatic DMR acquisition at implantation (Fig. 6). Similar interplays between these epigenetic modifications have recently been observed in mouse embryonic fibroblast cells (27), in which H3K27me3-enriched regions rapidly gain DNA methylation upon Dnmt3b overexpression, while preexisting H3K4me3 shields CGI from aberrant DNA methylation (27).

It should be noted, however, that the presence of paternal-specific H3K4me3 is not sufficient for establishing somatic DMR, as Sfmbt2 promoter has a paternal-specific H3K4me3 in TE, but it does not acquire DMR in extraembryonic lineage until after E6.5 (21). This variation in timing when to acquire somatic DMRs has been previously noticed. For example, establishment of germline DNA methylation–dependent somatic DMRs varies from E6.5 to E13.5 for different loci (4). Notably, the other three genes (i.e., G730013B05Rik, Tle3, and Grik3) with paternal allele–biased H3K4me3 in TE (Fig. 5) showed hypo-DNA methylation in E6.5 ExE and did not maintain imprinted expression or acquire somatic DMRs at later stages (9, 21). Therefore, preexisting allelic H3K4me3 is not the only determinate for future somatic DMRs. Establishment of an in vitro system that fully recapitulates this unique epigenetic regulation can help to understand the mechanism of allelic H3K27me3 to allelic DNA methylation switch at implantation.

Placenta-specific maintenance of noncanonical imprinting

Genomic imprinting, in general, is more prevalent in placental tissues (2). In particular, oocyte H3K27me3–dependent noncanonical imprinting is only maintained in the extraembryonic lineage (9). Analyses of the publicly available datasets (16) revealed that all five oocyte H3K27me3–dependent somatic DMRs are hypermethylated in E6.5 epiblast and the imprinted genes are expressed at a very low level in this lineage despite both alleles being transcribed (9). It is possible that the lack of a transcription activator(s) in epiblast causes loss of active chromatin marks, including H3K4me3, and ultimately leads to gain of DNA methylation on both alleles. In the extraembryonic lineage, however, this factor(s) maintains an active chromatin environment on the paternal allele to antagonize DNA methyltransferases in the wave of global de novo methylation at implantation. Identification of a relevant trophoblast-specific transcription factors is needed to test this hypothesis.

Precise control of the transcript dosage of noncanonical imprinted genes is important for placenta development. On the one hand, loss-of-function studies have revealed that depletion of some of the noncanonical imprinted genes severely impairs placental development (2830). For example, placentas of GAB1 KO embryos are generally smaller than WT, and GAB1 has been shown to regulate trophoblast cell number in the labyrinth region (30). On the other hand, biallelic expression and increased transcript abundance of these genes may lead to placental overgrowth phenotypes in cloned mice, which show loss of noncanonical imprinting (31).

Complementary roles of DNA methylation and H3K27me3 in regulating imprinting

DNA methylation and H3K27me3 represent two independent epigenetic silencing markers, and both can serve as the primary imprints to initiate genomic imprinting. Accumulating evidence indicates that these two imprinting mechanisms cooperate to maintain genomic imprinting. In addition to the allelic H3K27me3 to allelic DNA methylation switch described here, a reverse switch from primary DNA methylation–mediated imprinting to secondary allele–specific H3K27me3-mediated imprinting has also been described in both embryonic and extraembryonic lineage (3, 3234). It has been shown that secondary allele–specific H3K27me3 can take over the primary DMR to maintain imprinted expression for a few imprinted genes in extraembryonic lineage (33). It is fascinating that DNA methylation and H3K27me3 are used complementally to serve as primary and/or secondary imprinting marks to safeguard the normal development of placenta. Last, unlike the autosomal noncanonical imprinted loci, maintenance of imprinted X inactivation does not require DNA methylation in placenta (35). Therefore, although these loci and Xist use maternal H3K27me3 to initiate imprinted expression in preimplantation embryos, they adopt distinct mechanisms for imprinting maintenance in postimplantation development.

MATERIALS AND METHODS

All animal studies were performed in accordance with the guidelines of the Institutional Animal Care and Use Committee at Harvard Medical School.

Collection of DNMT3L WT and matKO morulae

The Dnmt3lG/+ (B6.129S6-Dnmt3ltm1Bes/J) mouse line (6) was purchased from the Jackson Laboratory (004267). The Dnmt3l+/+ (i.e., WT) and Dnmt3lG/G (i.e., KO) female mice were obtained by intercrossing of the Dnmt3lG/+ mice. For superovulation, DNMT3L WT and KO female mice (6 to 8 weeks old) were injected intraperitoneally with 7.5 IU of pregnant mare serum gonadotropin (PMSG; Millipore) on day 1 and 7.5 IU of human chorionic gonadotropin (hCG; Millipore) on day 3 (44 to 48 hours after PMSG injection). For in vitro fertilization (IVF), both DNMT3L WT and KO oocytes were fertilized with PWK/PhJ (PWK, Jackson stock no. 003715) strain sperm. Specifically, spermatozoa collected from the caudal epididymis of PWK male mice (9 to 12 weeks old) were first capacitated by incubation for 1 hour in human tubal fluid (HTF) medium supplemented with bovine serum albumin (BSA; 10 mg/ml; Sigma). The activated spermatozoa were then used to inseminate oocytes collected 12 to 16 hours after hCG injection. Approximately 6 hours post-IVF (hpf), zygotes were transferred from HTF medium to KSOM (Millipore) and cultured at 37°C under 5% CO2 with air. The time when sperm was added to the cumulus oocyte complexes was considered as 0 hpf. By 2 hpf, most oocytes were fertilized, judged by the presence of second polar body or fertilization cone. At 8 to 11 hpf, fertilized zygotes with two pronuclei were selected for further development. At ~78 hpf, ~7 WT (Dnmt3l+/+) and matKO (Dnmt3lG/+) morulae were collected for preparation of RNA-seq libraries. For sample collection, embryos were briefly incubated in Acid Tyrode (Millipore) to remove zona pellucida and then washed three times in 0.2% BSA/phosphate-buffered saline (PBS) before being transferred to a 0.2-ml tube with ~0.5-μl carryover. Two independent IVF experiments were performed for the RNA-seq replicates.

Collection of EED matKO morulae, TE, and E6.5 ExE

Generation of oocyte-specific EED KO mice using the Eedflox and Gdf9Cre mice was as described (10). For collection of EED CTR (Eedflox/flox × PWK) and EED matKO (Gdf9Cre, Eedflox/flox × PWK) embryos, superovulation and IVF were performed as described in the Dnmt3l section. Morulae and TE were collected at ~78 and ~120 hpf, respectively. Mural TE was physically dissected from inner cell mass. Around 30 morulae/TE were collected for the preparation of each CUT&RUN library. Independent IVF experiments were performed for each CUT&RUN replicate. To collect E6.5 ExE, EED female mice were co-caged with PWK male mice, and the day of vaginal plug was counted as E0.5. Four ExEs were pooled for each WGBS library. All samples were washed three times in 0.2% BSA/PBS before being transferred into a 0.2-ml tube with ~0.5-μl carryover.

DNMT3A and DNMT3B depletion by zygotic CRISPR-Cas9 injection

Zygotic depletion of DNMT3A and DNMT3B by CRISPR-Cas9 was performed following a previous study (16) with modifications. Superovulation and IVF were similar as described in the Dnmt3l section except that B6D2F1 (BDF1, B6 × DBA) female mice (6 to 8 weeks old) (the Jackson Laboratory, 100006) were used. The single-guide RNAs (sgRNAs) used to target Dnmt3a and Dnmt3b were the same as previously reported (3× sgRNA for each gene) (16). For one-cell microinjection experiments, Cas9 mRNA (100 ng/μl) and six sgRNAs (50 ng/μl each) were injected into the cytoplasm of zygotes (~2 to 4 hpf) using a Piezo impact-driven micromanipulator (Primer Tech, Ibaraki, Japan). For WT group, water was injected to serve as a control. Following injection, embryos were cultured in HTF medium for another 4 hours and then transferred to KSOM (Millipore) and cultured at 37°C under 5% CO2 with air. At ~24 hpf, embryos were transferred into oviducts of surrogate ICR strain mothers (Charles River). At E6.5, each ExE was dissected and split into two pieces: The smaller one (about one-third of ExE) was used for complementary DNA (cDNA) synthesis, followed by RT-qPCR and RNA-seq, while the larger piece (about two-thirds of ExE) was collected for low-input RRBS or H3K27me3 CUT&RUN. Two independent microinjection experiments were performed.

cDNA synthesis, RT-qPCR, and RNA-seq

The SMARTer Ultra Low Input RNA cDNA Preparation Kit (Clontech) was used for reverse transcription and cDNA amplification following the manufacturer’s instructions. The numbers of cycles used for cDNA amplification were 10 and 13 for Dnmt3l morulae and E6.5 ExE, respectively, which yielded cDNA (0.5 to 2 ng/μl) in 20 μl. For E6.5 ExE samples, ~10 pg of cDNA was used as templates for SYBR Green gene expression assay (Invitrogen) in the ViiA 7 Real-Time PCR System (Thermo Fisher Scientific). The threshold cycles of Dnmt3a and Dnmt3b were normalized to Gapdh using the comparative CT method. Only samples that showed minimal level of both Dnmt3a and Dnmt3b were considered as DKO. For both Dnmt3l morula and E6.5 ExE, 400 pg of cDNA obtained from the SMARTer step was used to generate RNA-seq libraries. The Nextera XT DNA Library Preparation Kit (Illumina) was used for cDNA fragmentation, adaptor ligation, and amplification (12 cycles) according to the manufacturer’s instructions. Single-ended 100-bp sequencing was performed on a HiSeq 2500 sequencer (Illumina).

Low-input WGBS

For low-input WGBS library, E6.5 ExE samples (four for each library) lysed in Buffer RLT Plus (Qiagen) were subject to sodium bisulfite conversion using the Zymo EZ DNA Methylation-Direct Kit (Zymo Research) following the manufacturer’s instructions. The eluted bisulfite-converted DNA (BS-DNA) (19.5 μl) was then subject to random tagging by incubation at 65°C for 3 min after the addition of 2.5 μl of 10× Blue buffer (supplied with Klenow exo−, Enzymatics), 1 μl of 10 μM PreAmp oligo (5′-CTACACGACGCTCTTCCGATCTNNNNNNNNN-3′), and 1 μl of 10 mM deoxynucleotide triphosphates (dNTPs) (Thermo Fisher Scientific). After cooled down on ice, extension step was followed by adding 1 μl of Klenow exo− (50 U/μl; Enzymatics) and incubated at 4°C for 5 min and 37°C for 30 min (1°C/15 s from 4° to 37°C). After denaturation at 95°C for 1 min and cooled down on ice, the same extension program was repeated after the addition of 0.5 μl of Klenow exo− (50 U/μl), 0.1 μl of 10 mM dNTPs, 1 μl of 10 μM PreAmp primer, 0.25 μl of 10× Blue buffer, and 0.65 μl of H2O. The denaturation and extension steps were repeated for another three rounds, which brought the reaction volume to 35 μl. Any remaining oligos were removed by incubating at 37°C for 1 hour with the addition of 2 μl of Escherichia coli Exonuclease I (20 U/μl; New England BioLabs) and 63 μl of H2O. After preamplification, adaptor 2 tagging was followed. Specifically, purified DNA from the preamplification step (20 μl) was denatured at 95°C for 1 min after the addition of 20 μl of H2O, 5 μl of 10× Blue buffer, 2 μl of 10 μM Adaptor 2 oligo (5′-CAGACGTGTGCTCTTCCGATCTNNNNNNNNN-3′), and 2 μl of 10 mM dNTPs. After cooled down on ice, extension step was followed with the addition of 1 μl of Klenow exo− (50 U/μl) and incubation at 4°C for 5 min and 37°C for 90 min (+1°C/15 s from 4° to 37°C). Last, the double-tagged DNA was purified by SPRI beads (Beckman Coulter) and amplified by six cycles using KAPA HiFi Uracil+ ReadyMix (Kapa Biosystems) and NEBNext Multiplex Oligos for Illumina (New England BioLabs). The final libraries were subject to single-ended 150-bp sequencing on a HiSeq 2500 sequencer (Illumina).

Low-input RRBS

For low-input RRBS, around two-thirds of single ExE were lysed at 50°C for 3 hours and 75°C for 30 min in a 5-μl lysis buffer [20 mM tris-EDTA (Sigma), 20 mM KCl (Ambion), 0.3% Triton X-100 (Sigma), and protease (1 mg/ml; Qiagen)]. The lysed cells were then subject to Msp I digestion at 37°C for 3 hours and 80°C for 20 min after the addition of 10.2 μl of H2O, 1.8 μl of 10× TangoBuf (Thermo Fisher Scientific), and 1 μl of Msp I (10 U/μl; Thermo Fisher Scientific). Following enzymatic digestion, end-repair mix [1 μl Klenow frag exo− (5 U/μl), 0.2 μl of 10× TangoBuf, and 0.8 μl of dNTP mix (1 mM deoxyadenosine triphosphate, 0.1 mM deoxycytidine triphosphate, and 0.1 mM deoxyguanosine triphosphate)] was added and incubated at 37°C for 40 min and 75°C for 15 min. For adaptor ligation, the reaction was incubated at 16°C for 30 min, 4°C for at least 8 hours, and 65°C for 20 min after adding the ligation mix [2.25 μl of H2O, 1 μl of T4 ligase (30 U/μl; Thermo Fisher Scientific), 0.5 μl of 100 mM adenosine triphosphate (Thermo Fisher Scientific), and 0.75 μM methylated adaptor (forward: 5′-ACACTCTTTCCCTACACGACGCTCTTCCGATC*T-30; reverse: 5′-/5Phos/GATCGGAAGAGCACACGTCTGAACTCCAGTC-3′, where * indicates phosphorothioate bond)]. Bisulfite conversion reaction was then performed by the EpiTect Bisulfite Kit (Qiagen) following the manufacturer’s instructions. Note that “restriction digested DNA” protocol should be used for recovery of the BS-DNA. Last, adapter-tagged BS-DNA was amplified by 16 cycles using Kapa HiFi Uracil+ ReadyMix (Kapa Biosystems) and NEBNext Multiplex Oligos for Illumina (New England BioLabs) and purified using SPRI beads (Beckman Coulter). The RRBS libraries were subject to 100-bp paired-ended sequencing on an Illumina HiSeq 2500 system (Illumina).

H3K27me3 and H3K4me3 CUT&RUN

CUT&RUN libraries were prepared following the same protocol as previously described (10, 13). For H3K27me3 CUT&RUN, a rabbit anti-H3K27me3 antibody (1:100; Diagenode, C15410069) was used. For H3K4me3 CUT&RUN, a rabbit anti-H3K4me3 antibody (1:100; Cell Signaling, 9727S) was used. All CUT&RUN libraries were sequenced on a HiSeq2500 with paired-ended 100-bp reads (Illumina).

RNA-seq analyses

Low-quality bases and sequencing adaptors of raw RNA-seq reads were trimmed using Trim Galore (version 0.4.5) with parameters “-q 20 –length 20 bp.” Reads after trimming were then aligned to the mm9 reference genome using HISAT2 (version 2.1.0) (36) with parameters “--sp 1000,1000 --dta-cufflinks --no-unal.” Known exon-exon splicing site information was provided for HISAT2 to improve alignment of reads that span splicing junctions. To minimize the alignment bias toward the reference genome (B6), SNPs between B6/129 and PWK strains were masked as “N.” As both Dnmt3l and Eed transgenic mouse lines were from founders of a mixed background of B6 and 129, only SNPs that are common between B6 and 129 but different from PWK (n = 18,361,596) were used to differentiate parental alleles. For Dnmt3a/3b depletion experiments, the genetic background of mouse embryos is B6/DBA × PWK. Therefore, SNPs that are common between B6 and DBA but different from PWK (n = 18,317,863) were used for these datasets. The known SNPs were retrieved from www.sanger.ac.uk/science/data/mouse-genomes-project. Analyses of EED CTR and matKO E6.5 ExE RNA-seq data (10) were the same as Dnmt3l RNA-seq. Analyses of B6 × PWK and PWK × B6 E6.5 ExE RNA-seq data (9) were similar as described above except that the only SNPs between B6 and PWK were masked as “N” in the reference assembly.

Following reads mapping, uniquely aligned reads were obtained by selecting the alignment records with “NH:i:1” tag and mapping quality greater than 30. SNPsplit (version 0.3.2) (37) was used to determine the parental origins of uniquely aligned reads, and all SNP-containing reads for each gene were summarized for allele-specific expression analyses. The bigwig tracks were generated using bamCoverage from deepTools (version 3.0.2) (https://deeptools.readthedocs.io/en/develop/) with parameters “--skipNonCoveredRegions --binSize 10 --scaleFactor 1.” For differential gene expression analyses, raw read counts for each gene were generated using TEtranscripts (version 1.5.1) (38) and statistical significance (i.e., FDR) was computed using the DESeq R package (version 1.34.1) (39). Only genes with RPKM > 1 in either group, FDR < 0.05, and FC > 2 were considered as differentially expressed.

WGBS data analyses

WGBS raw reads were trimmed using Trim Galore (version 0.4.5) with parameters “-q 20 –length 35 bp -clip-R1 9 –three_prime_clip_R1 9.” Following reads trimming, Bismark (version 0.17.0) (40) was used to align the reads to the mm9 reference genome with default parameters. B6/129 and PWK SNPs were masked as “N” for the reference assembly. After excluding nonuniquely aligned reads and removing PCR duplicates using “MarkDuplicates” from Picard Tools (version 2.8.0) (https://broadinstitute.github.io/picard/), methylation levels for each CpGs were computed as “reads containing methylated C/(reads containing methylated C + reads containing unmethylated C).” For comparison with public E6.5 ExE WGBS data (16), methylation levels of 100-bp tiles with >2 CpGs and >10 reads coverage in both samples were used to compute the Pearson correlation coefficient. DMRs were identified using MethPipe (version 3.4.3) (41) with the cutoff of >10 CpGs and >25% methylation differences between two groups. The methylation level of each hypo-DMR and hyper-DMR is the average of all CpGs (>5 reads) within each region.

Parental origins of the WGBS reads were determined using SNPsplit (version 0.3.2) (37), and “bedGraphToBigWig” from ENCODE (www.encodeproject.org/software/bedgraphtobigwig/) was used to convert bedGraph format to bigwig tracks. The minimum reads coverage for nonallelic and allelic DNA methylation tracks is five and three, respectively. For nonallelic and allelic DNA methylation at known ICRs, the chromosome coordinates were retrieved from (25).

RRBS data analyses

RRBS paired-end reads were trimmed using Trim Galore (version 0.4.5) with parameters “--rrbs --paired -q 20 –length 20.” Bismark (version 0.17.0) (40) was used to align the trimmed reads to the mm9 reference genome with parameters “-X 1000.” B6/DBA and PWK SNPs were masked as “N” for the reference assembly. Uniquely aligned reads were used to compute methylation levels of each CpGs similar to the analyses of WGBS datasets. The RRBS reads from the three replicates of WT and DKO were pooled for further analyses. Generation of nonallelic and allelic DNA methylation tracks were the same as the WGBS datasets except that only CpGs with >5 reads coverage were used for both types of tracks. For comparison with WGBS data of TE and RRBS data of DNMT3A and DNMT3B single-KO ExE (16), only CpGs covered by >5 reads in all samples were retrieved and 100-bp tiles with >2 CpGs and >10 reads coverage in all samples were used for global DNA methylation analyses.

H3K27me3 and H3K4me3 CUT&RUN analyses

For both CUT&RUN datasets, paired-end reads were trimmed using Trim Galore (version 0.4.5) with parameters “--paired -q 20 -length 20.” Following reads trimming, Bowtie2 (version 2.2.9) (42) was used to align the reads to the mm9 reference genome with parameters “--no-unal -I 150 -X 800 --no-mixed --no-discordant.” For H3K27me3 datasets, B6/DBA and PWK SNPs were masked as “N” for the reference assembly. For H3K4me3 datasets, B6/129 and PWK SNPs were masked as “N” for the reference assembly. PCR duplicates were removed using “MarkDuplicates” from Picard Tools (version 2.8.0) (https://broadinstitute.github.io/picard/). Only non-PCR duplicates and uniquely aligned reads (alignment records without “XS” tag) were used for downstream analyses. For comparing reproducibility between replicates and/or publicly available datasets, read counts over each bin (2 kb for H3K4me3 and 5 kb for H3K27me3) were used. The RPKM values for each bin were calculated following the formula “read counts/((bin_length/1000) × (total_reads/106)).” Pearson correlation coefficient was used to assess the reproducibility between H3K27me3/H3K4me3 profiling experiments. After confirming the reproducibility, CUT&RUN replicates were pooled for downstream analyses. SNPsplit (version 0.3.2) (37) was used to determine the parental origins of the read pairs, and the bigwig tracks were generated using bamCoverage from deepTools (version 3.0.2) (https://deeptools.readthedocs.io/en/develop/) with parameters “--binSize 20 -e 250 --scaleFactor 1.”

To identify H3K27me3 domains in E6.5 ExE, macs2 (version 2.1.1) (43) was used for initial peak calling with parameters “-f BAMPE -g mm -q 0.01 --broad --broad-cutoff 0.1 --nolambda --SPMR --nomodel –B.” H3K27me3 peaks within 5-kb distance were then merged into H3K27me3 domains. To filter out weak domains, raw read counts for each H3K27me3 domain were generated and RPKM values were calculated as described above. Only H3K27me3 domains with RPKM values greater than the median of all domains were kept for further analyses (WT median RPKM = 1.2 and KO median RPKM = 1.0). H3K27me3 domains in WT and DKO ExE were further merged for the purpose of comparing average H3K27me3 level at the identified domains between the two groups. To investigate the inheritance of oocyte H3K27me3 domains in morulae and E6.5 ExE, previously identified 4134 maternal H3K27me3 domains in IVF morulae (B6/DBA × PWK) (31) were used.

To compare H3K4me3 enrichment in morulae and TE between EED CTR and matKO, macs2 (version 2.1.1) (43) was used for peak calling with parameters “-f BAMPE -g mm -q 0.01 --nolambda --SPMR –nomodel.” H3K4me3 peaks in CTR and matKO were then merged, and the RPKM value for each peak was computed as described above. Only peaks with RPKM >2 in either CTR or matKO were used for the comparative analyses. For allelic H3K4me3 peak analyses at the putative noncanonical imprinting loci (n = 78) (9), parental-specific read counts of promoter H3K4me3 peaks (raw macs2 peaks without any filtering) were generated. If multiple H3K4me3 peaks exist at the transcription start site, the strongest peak was selected for further analyses. For each peak with >50 SNP-containing reads, “log2(paternal/maternal)” was computed, as this internal comparison of allelic H3K4me3 enrichment minimizes variabilities between CUT&RUN experiments.

Statistical analyses and data visualization

Statistical analyses were performed in R (www.r-project.org/). Pearson correlation coefficient was computed using the “cor” function in R (figs. S1A, S2A, S4B, S5A, and S6, A, B, E, and F). Figures 1B, 2 (C and D), 3G, and 5C and fig. S4D were generated using the R function “pheatmap.” All sequencing tracks were visualized using Integrated Genome Browser (https://software.broadinstitute.org/software/igv/).

Data availability

All RNA-seq, low-input WGBS and RRBS, and H3K4me3 and H3K27me3 CUT&RUN datasets generated in this study are summarized in table S4 and have been deposited in the Gene Expression Omnibus under accession number GSE130115. Oocyte H3K27me3 ChIP-seq data were obtained from GSE73952 (23), and sperm H3K27me3 ChIP-seq data were from GSE76687 (14). DNMT3L WT and KO WGBS data were from DRA000570 (44), and sperm WGBS data were from DRA000484 (45). B6 × PWK and PWK × B6 E6.5 ExE RNA-seq data were from GSE92605 (9), and EED CTR and matKO ExE RNA-seq data were from GSE116713 (10). WGBS data of TE/ExE and RRBS data of DNMT3A/3B/3L single-KO ExE were from GSE84236 (16). Morula (B6/DBA × PWK) H3K27me3 ChIP-seq data were from GSE112546 (31), and H3K4me3 ChIP-seq data of one-, two-, four-, and eight-cell embryos were from GSE71434 (24). H3K4me3 ULI-NChIP-seq data of morula and TE were obtained from GSE73952 (23).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/12/eaay7246/DC1

Fig. S1. Modest gene expression alterations in DNMT3L matKO morulae.

Fig. S2. Characterization of DNA methylome of EED CTR and matKO E6.5 ExE.

Fig. S3. Both DNMT3A and DNMT3B contribute to somatic DMR establishment at noncanonical imprinting loci.

Fig. S4. Validation of RNA-seq and RRBS datasets of DNMT3A/3B DKO ExE.

Fig. S5. H3K27me3 profile in E6.5 ExE.

Fig. S6. H3K4me3 profiles in morula and TE.

Table S1. RNA-seq analyses of DNMT3L CTR and matKO morulae.

Table S2. DNA methylation analyses of EED CTR and matKO E6.5 ExE.

Table S3. DNA methylation and RNA-seq analyses of DNMT3A/3B DKO E6.5 ExE.

Table S4. Allelic read counts of promoter H3K4me3 in morula and TE and summary of datasets generated in this study.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We would like to acknowledge X. Wu for help in WGBS, L. Shen for help in low-input RRBS, and W. Zhang for helpful discussion. The protein A–MNase for CUT&RUN was a gift from S. Henikoff. Funding: This project was supported by the NIH (R01HD092465). Y.Z. is an Investigator of the Howard Hughes Medical Institute. Author contributions: Y.Z. conceived the project. Z.C. and Y.Z. designed the experiments. Z.C. performed most of the experiments and analyzed sequencing datasets. Q.Y., A.I., and C.Z. helped in some experiments. Z.C. and Y.Z. interpreted the data and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article