N6-methyladenosine regulates maternal RNA maintenance in oocytes and timely RNA decay during mouse maternal-to-zygotic transition

N6-methyladenosine (m6A) and its regulatory components play critical roles in various developmental processes in mammals. However, the landscape and function of m6A in early embryos remain unclear owing to limited materials. Here we developed a method of ultralow-input m6A RNA immunoprecipitation followed by sequencing to reveal the transcriptome-wide m6A landscape in mouse oocytes and early embryos and found unique enrichment and dynamics of m6A RNA modifications on maternal and zygotic RNAs, including the transcripts of transposable elements MTA and MERVL. Notably, we found that the maternal protein KIAA1429, a component of the m6A methyltransferase complex, was essential for m6A deposition on maternal mRNAs that undergo decay after zygotic genome activation and MTA transcripts to maintain their stability in oocytes. Interestingly, m6A methyltransferases, especially METTL3, deposited m6A on mRNAs transcribed during zygotic genome activation and ensured their decay after the two-cell stage, including Zscan4 and MERVL. Together, our findings uncover the essential functions of m6A in specific contexts during the maternal-to-zygotic transition, namely ensuring the stability of mRNAs in oocytes and the decay of two-cell-specific transcripts after fertilization. Wu et al. optimized the m6A mapping method for ultralow-input materials and characterized the transcriptome-wide landscape of m6A methylation in mouse oocytes and early embryos.

F ertilization triggers a remarkable and complex cell fate transition in oocytes from terminal differentiation to a totipotent state, termed the maternal-to-zygotic (MZT) transition, during which post-transcriptional regulation is considered to be important for the storage, timely decay and sequential translational activation of maternal and zygotic transcripts [1][2][3] . Although some details of the reprogramming of the epigenetic and chromatin landscape have been described in recent years, the various mechanisms of RNA control remain elusive. N 6 -methyladenosine (m 6 A), which is found on messenger RNAs, repeat RNAs and long noncoding RNAs [4][5][6][7][8] , is considered to be important for various developmental events by regulating RNA-related processes [9][10][11][12][13][14] . A previous study in zebrafish found that m 6 A can promote maternal RNA degradation during the MZT process 15 . In mice, different studies have revealed the essential roles of m 6 A writer and reader proteins in oogenesis and early embryogenesis. For instance, depletion of METTL3, METTL14 or METTL16 in embryos from crossed heterozygous mice results in embryonic lethality at the early stage in post-implantation embryos [16][17][18] . Maternal knockout (KO) of Kiaa1429 (also known as Virma), a component of the RNA m 6 A methyltransferase complex, in oocytes results in failure of germinal vesicle (GV) breakdown and meiotic resumption 19 . In addition, maternal deletion of m 6 A reader YTHDF2 or IGF2BP2 in mouse embryos results in two-cell (2C) arrest 20,21 , and YTHDC1, a nuclear m 6 A reader, is indispensable for gamete maturation and embryo development 11,22 . All these findings reveal that regulating RNA m 6 A and determining the m 6 A-marked (m 6 A + ) RNA fate is quite important for oocyte maturation, ZGA and post-implantation development. However, the molecular dynamics of m 6 A on RNA in oocytes and early embryos remain unknown owing to limited materials, and the detailed regulatory model is largely ambiguous.
During the MZT, dynamic stage-specific expression of endogenous retroviruses (ERVs) or repeats is observed and is considered to regulate oocyte maturation and ZGA 1,23 . In mice, the mouse transcript (MT) retrotransposon, a member of the MaLR family of non-autonomous Class III ERV retrotransposons, accounts for 12% of all transcripts in fully grown oocytes as measured by expressed sequence tags in a complementary DNA library 24 . Knockdown (KD) of MT retrotransposon elements in GV oocytes results in aberrant meiotic resumption 25 . After fertilization, murine ERV-like (MERVL), a LTR Class III molecule, is transiently activated in the late one-cell (L1C) stage and considered to contribute to the activation of downstream genes in ZGA 1 . Recently, the transcripts of transposable elements (TEs) including LINE1 and IAP are found to N 6 -methyladenosine regulates maternal RNA maintenance in oocytes and timely RNA decay during mouse maternal-to-zygotic transition be marked by m 6 A and decayed through YTH proteins-mediated manner 10,26 and involved in the regulation of transcriptional activity and chromatin status 11,14,22,27 . However, the existence and role of m 6 A on TEs in early embryos are not well understood.
In this Article, we first developed a method of ultralow-input (ULI) m 6 A RNA immunoprecipitation (MeRIP) followed by sequencing (ULI-MeRIP-seq) that enabled the profiling of m 6 A RNA methylation with 50 ng of total RNA. With this method, we revealed the dynamics of the m 6 A landscape of the transcriptome during MZT progression from GV oocytes to four-cell (4C)-stage embryos. Here we defined 1,089 maternal genes and 756 zygotic genes covered by m 6 A in oocytes and after fertilization, respectively, and observed the preference of m 6 A on zygotic decay (Z-decay) genes and ZGA genes with high abundance. Among repetitive elements, we observed highly enriched m 6 A on the MTA and MERVL. With the advantage of maternal KO mice, we found that m 6 A established through KIAA1429 may contribute to the stabilization of Z-decay mRNAs and MTA transcripts in oocytes. Furthermore, we found that METTL3 played an important role in the establishment of m 6 A on ZGA transcripts and MERVL, which further promoted the timely turnover of 2C-specific markers and ensured the developmental progression of pre-implantation embryos.

ULI-MeRIP-seq for m 6 A.
To investigate the role of RNA m 6 A methylation during the MZT in mammals, we developed ULI-MeRIPseq based on a published MeRIP-seq method that required 500 ng to 2 µg of total RNA 28 . Here we downscaled the total RNA amount from 2 µg to 50 ng with several improvements. Briefly, MaXtract high-density tubes and glycogen were used to increase the recovery of RNA and the removal of genomic DNA during TRIzol-based RNA extraction and immunoprecipitation (IP) product purification. To avoid frequent purification and increase IP efficiency, we used sonication-RNA fragmentation instead of the zinc-mediated assay (Extended Data Fig. 1a, b). We further optimized the usage of the m 6 A antibody for ULI RNAs to increase the IP efficiency, which was measured by the signal-to-noise (SN) ratio of endogenous transcripts (Extended Data Fig. 1c). With these improvements, we were able to examine the m 6 A RNA methylome of mRNAs and noncoding RNAs using 50 ng of total RNA per IP reaction, which can be extracted from approximately 150 mouse oocytes (Fig. 1a), and 1/20 of the sonicated RNA were used to construct the input libraries.
To evaluate our data quality, we used 50 ng and 2 µg of total RNA from mouse embryonic stem cells (mESCs) for ULI-MeRIP-seq, respectively. Quality control analysis by quantitative PCR (qPCR) showed that both 2 µgand 50 ng-input RNA libraries exhibited relatively high SN ratios of spike-in RNAs Glaussia luciferase (GLuc)/Cypridina luciferase (CLuc) as well as endogenous Klf4/ Gapdh transcripts compared with the published SN values 28 (Fig. 1b). The modification levels on pluripotency marker genes were consistent with previous reports 16,29 (Extended Data Fig. 1d). We calculated the ratio of antisense reads on coding genes, and confirmed that the genomic DNA contamination in our data was comparable to that in a previous study with DNase I treatment 26 (Extended Data Fig. 1e). The data generated using 50 ng of total RNA recapitulated the results generated using 2 µg of total RNA, with high similarity according to peak, motif and genome browser tracks (Fig. 1c,d and Extended Data Fig. 1f-i). These results supported the use of ULI-MeRIP-seq for oocyte and early embryo samples.

Dynamic characteristics of the m 6 A methylome during the MZT.
Recent studies have revealed that m 6 A on mRNAs may be involved in mouse oocyte and embryo development 18,20,21 . To detect the quantity changes of RNA m 6 A during the mouse MZT, we performed liquid chromatography with tandem mass spectrometry (LC-MS/MS) of total RNA samples extracted from oocytes at the GV and MII (metaphase of second meiosis) stages, embryos at the L1C, L2C and 4C stages, as well as ESCs (Fig. 1e). We found that the overall abundance of m 6 A continuously decreased during the MZT (Fig. 1f). To describe the whole-transcriptome RNA m 6 A methylome profile, we employed ULI-MeRIP-seq using 300-750 oocytes or embryos from each stage, which included two or three replicates for IP and input libraries (Supplementary Table 1). The SN ratio of GLuc/CLuc spike-in RNA was checked to ensure high enrichment for m 6 A for each IP reaction (Extended Data Fig. 2a and Supplementary Table  2). Library complexity analysis indicated that PCR duplicates were less than the 500 ng input MeRIP-seq previously reported and none of the input or IP samples showed severe PCR bottlenecking (PCR bottleneck coefficient 2 >1) according to ENCODE standards 28 (Extended Data Fig. 2b-e), and all replicates were highly comparable (Extended Data Fig. 2f, g). Saturation analysis showed that our data were able to detect both high-level and low-level transcripts of genes (Extended Data Fig. 2h). Principal component analysis (PCA) of genes showed that both the transcriptome and m 6 A methylome diverged after the 2C stage when major ZGA occurred (Fig. 1g). To identify m 6 A-enriched sites, we compared peak calling using model-based analysis of ChIP-seq (MACS) 30 and MeTPeak 31 model-based analyses. The results showed that MACS captured most (more than 80%) of the m 6 A-modified transcripts that were detected by MeTPeak in all samples (Extended Data Fig.  2i). Considering that MACS can also be used to identify peaks in intergenic regions that may include highly active retrotransposon elements in early embryos 5,26 , we chose MACS for peak calling in this study.
As expected, the peak length distribution and the enrichment of the canonical RRACH motif (R = G/A; H = A/C/U) on m 6 A peaks at each stage were consistent with those in previous studies (Extended Data Fig. 2j,k) 4,5 ; however, the density of RRACH at exon peaks decreased throughout the MZT process (Extended Data Fig. 2l). In addition to the well-characterized enrichment of m 6 A methylation in coding sequences (CDSs), in untranslated regions (UTRs), and near stop codons (Fig. 1h), we found that a considerable number of m 6 A peaks mapped to distal intergenic regions, which was particularly high at the L1C stage (Extended Data Fig. 2m).
To evaluate the details of m 6 A dynamics, we analysed the m 6 A + coding transcripts and peaks at each stage. Overall, the number of m 6 A + transcripts increased gradually after fertilization, which doubled at L1C (Extended Data Fig. 3a-d). By comparing the dynamics of m 6 A + transcripts and peaks, we found that the major gain and loss of m 6 A + peaks at the MII-to-L1C and L1C-to-L2C transition was not consistent with that on m 6 A transcripts (Fig. 1i, j). Further analysis showed that the location of gained peaks at L1C stage and lost peaks at L2C and 4C stage preferred transposon regions (Fig. 1k), which suggests that m 6 A may regulate the fate of transposon RNAs during the MZT process.
We noticed that, while a continuously decreased m 6 A/A ratio was detected by LC-MS/MS, m 6 A + transcripts and peaks increased after fertilization. To identify the potential reason, we examined the expression patterns of the transcripts gained m 6 A at L1C and L2C, and found most of them were 2C-activated, suggesting that de novo m 6 A modification was established along with ZGA (Extended Data Fig. 3e). These results prompt us to further profile the m 6 A dynamics of maternal decay (MD) and ZGA transcripts. m 6 A marks Z-decay and ZGA transcripts during the MZT. During the MZT process, the majority of maternal RNAs decay during oocyte maturation and fertilization, and the first wave of zygotic genome activation occurred from L1C to L2C stage. Using transcriptome data, we identified 2,533 MD genes that were highly expressed in oocytes but dramatically downregulated in embryos, and 1,115 ZGA genes that showed no or low expression in oocytes but were upregulated after fertilization ( Fig. 2a Table 3; see Methods for details). We then found that 43% of MD transcripts were marked by m 6 A in GV or MII oocytes (Fig. 2a,b). Approximately 68% of ZGA transcripts gained m 6 A in the L1C or L2C stage, almost along with the transcription activation ( Fig. 2a,b). Gene Ontology (GO) analyses showed that m 6 A-unmarked (m 6 A − ) MD transcripts were significantly enriched in energy metabolism pathways, and m 6 A + MD transcripts were mainly associated with TGF-β, epithelial-to-mesenchymal transition and GTPase signalling pathways, which play essential roles in oogenesis such as Gdf9, Tgfb2, Jun, Mos and Bmp15 (Fig. 2c,d and Extended Data Fig. 3f).
The m 6 A − ZGA transcripts predominantly enriched for acid chemical and amino sugar metabolic processes, and m 6 A + ZGA transcripts were strongly enriched for processes involved in blastocyst development and transcription factor activity ( Fig. 2c and Extended Data Fig. 3g). Notably, the transcripts of 2C marker genes, including Dux and Zscan4, were found to be marked by m 6 A ( Fig. 2d and Extended Data Fig. 3g). Analysis of RRACH density around the stop codon revealed that the m 6 A + MD and ZGA transcripts possessed higher RRACH densities (Extended Data Fig. 3h), which partly explained the preference for m 6 A.   Previous studies have revealed that maternally encoded transcript decay occurs before or after ZGA, which is termed maternally encoded mRNA decay before ZGA (M-decay) or Z-decay (decay after ZGA) 32,33 , respectively. We thus wondered if the enrichment of m 6 A was different between M-decay and Z-decay transcripts. We applied previously described criteria to define M-decay and Z-decay genes 32 (Methods) and found that the majority of previously identified MD genes fell into the M-decay and Z-decay classes (Extended Data Fig. 3i). The percentage of m 6 A + transcripts in Z-decay classes was significantly higher than that of M-decay classes (Fig. 2e). Further analysis showed that the RNA level of m 6 A + transcripts was higher than the unmarked counterparts in Z-decay and ZGA classes, but the difference was not significant in M-decay classes (Fig. 2f,g). These results indicated that m 6 A methylation was preferentially performed on Z-decay and ZGA transcripts with elevated expression levels. m 6 A is enriched on stage-specifically expressed TEs. Previous studies have reported that certain retrotransposons such as MaLR, MERVL and LINE exhibited stage-specific expression and regulatory function during oocyte and embryo development 1,23,24 . In addition, we observed the existence of considerable m 6 A turnover at intergenic TEs (Fig.1k). We therefore analysed the dynamics of m 6 A peaks on transposons at certain stages, which required more detailed investigation on TEs.
To start this analysis, we first isolated the reads located on repetitive RNAs and then calculated the enrichment of m 6 A on different classes of repeats (Extended Data Fig. 4a and Methods). Generally, m 6 A methylation was mainly enriched in LTR RNAs in all MZT samples (Fig. 3a). However, m 6 A in ESCs was mainly enriched in LINE RNAs (Fig. 3a). Subsequently, we analysed the m 6 A enrichment and expression of the main families of LTRs and LINEs. We found that MaLR and ERVL, which are both LTR class III retrotransposons, possessed high m 6 A enrichment in oocytes and early cleavage embryos, respectively (Fig. 3b), consistent with their stage-specific expression features (Extended Data Fig. 4b) 24,34 .
To be more precise, we defined all the expressed TE subfamilies during MZT progression (Methods) and classified them into three clusters, termed maternal TEs, ZGA TEs and mid-pre-implantation gene activation (MGA) TEs, based on their expression patterns from GV oocytes to blastocyst embryos (Fig. 3c). Most of the expressed maternal TEs, mainly members of the MaLR family, possessed high RNA abundance and m 6 A enrichment from the oocyte to L1C embryo stages (Fig. 3c). Compared with MGA TEs, almost all ZGA TEs, including MERVL, MT-int and ORR1As, exhibited high RNA and m 6 A levels from L1C to 4C embryos (Fig. 3c). This observation revealed that m 6 A marked the transcripts of active TEs along with their expression-stage specificity during the MZT process. To investigate the detailed m 6 A dynamics, we divided all individual full-length MTA and MERVL copies into five clusters based on their expression levels in GV and 2C stage, respectively. We found that the occupancy of m 6 A positively correlated with the expression level, and almost all transcripts of high-expressed MTA and MERVL were marked by m 6 A (Fig. 3f,g and Extended Data Fig. 4e,f). Motif analysis showed that UGAC, which was reported as a common m 6 A motif variant on mRNA 35 , was highly enriched within m 6 A peaks on MTA, and the classic GGAC motif was enriched within m 6 A peaks on MERVL (Extended Data Fig. 4g,h).
The profiling of m 6 A on MTA during the MZT showed not only continuous demethylation after the zygote stage but also a significant m 6 A density shift from transcription end sites (TESs) to transcription start sites (TSSs) in MII oocytes and L1C embryos ( Fig. 3h and Extended Data Fig. 4i). Considering that the MTA RNAs stored in GV oocytes are readily decayed in matured or fertilized oocytes, the shift in m 6 A peaks suggests that region-specific m 6 A deposition might be important for m 6 A-mediated regulation. For MERVL, the deposition of m 6 A appeared at the L1C stage and increased in the 2C and 4C stages ( Fig. 3i and Extended Data Fig. 4j), which seemed to be established in a co-transcriptional manner. We further confirmed these results using the unique mapping reads (Extended Data Fig. 4k-n). These results reveal that highly abundant m 6 A methylation is established on MTA and MERVL, which may be further involved in the post-transcriptional regulation of these maternal and ZGA TEs.

KIAA1429-mediated m 6 A stabilizes Z-decay mRNAs in oocytes.
To explore the regulatory function of m 6 A on maternal RNA, we took advantage of Kiaa1429 conditional KO (cKO) mouse model, which exhibited overall m 6 A reduction, abnormal RNA metabolism and GV breakdown failure in oocytes 19 . To underline the defects caused by the depletion of KIAA1429 in oocytes, we performed RNA sequencing (RNA-seq) and ULI-MeRIP-seq of GV oocytes from Kiaa1429 Zp3 control and cKO mice (Fig. 4a). We found the number of m 6 A + transcripts in GV oocytes decreased dramatically with widespread transcriptome disorder occurred, when KIAA1429 was deleted ( Fig. 4b and Extended Data Fig. 5a-d). Further analysis revealed that the abundance of m 6 A + Z-decay transcripts was reduced significantly in KIAA1429-deficient oocytes, but the abundance of M-decay and unmodified Z-decay transcripts was almost unchanged (Fig. 4c). We then identified the differentially expressed genes (DEGs; see cut-off in Methods), 1,869 downregulated and 2,228 upregulated genes, upon KIAA1429 depletion (Extended Data Fig. 5e and Supplementary Table 4). We found that a very limited number of DEGs overlapped with m 6 A + M-decay transcripts. In contrast, 94.9% of differentially expressed m 6 A + Z-decay transcripts were downregulated (282/297), suggesting the potential role of m 6 A in maintaining high expression levels of these genes before fertilization (Fig. 4d).
Previous studies have indicated that m 6 A deficiency in chromosome-associated regulatory RNA caused by Mettl3 KO promotes the negative regulation of chromatin openness and the transcription activity in ESC lines 10 . We therefore investigated whether the repression of m 6 A + Z-decay transcripts in Kiaa1429 Zp3 -cKO oocytes were combined with the change on chromatin openness. We therefore performed the assay for transposase-accessible chromatin using sequencing (ATAC-seq) with the nuclei of control and Kiaa1429 Zp3 -cKO GV oocytes, in which the non-surrounded nucleolus (NSN) and partly surrounded nucleolus (PSN) GV oocytes were collected separately (30-50 nuclei per reaction). In general, the ATAC-seq data show high repeatability and genome coverage to those in published studies 36 (Extended Data Fig. 5f,g). Detailed calculation revealed that the ATAC-seq signals at the promoter regions (+1.5 kb and −0.5 kb around TSS) of the M-decay and Z-decay genes did not change significantly upon depletion of Kiaa1429 (Fig. 4e). Therefore, the reduction of m 6 A + Z-decay transcripts are probably caused by reduced RNA stability rather than transcription activity.
KIAA1429-mediated m 6 A stabilizes MTA RNAs in oocytes. In addition to the Z-decay transcripts, we also observed high enrichment of m 6 A on MTA RNAs in oocytes. We therefore assessed the behaviour of MTA transcripts in Kiaa1429 Zp3 -cKO oocytes. As expected, the m 6 A modification on MTA RNA, especially the peak near the TES, decreased greatly when KIAA1429 was depleted (Fig. 4f). In addition to the loss of m 6 A, the dramatic decrease in RNA level for almost all MTA copies was remarkable and much more drastic than that observed for m 6 A + transcripts (Fig. 4g,h). We also compared the ATAC-seq signals of the clustered MTA with different expression levels and found that the transcriptional activity of MTA remained high in both NSN and PSN oocytes when KIAA1429 was depleted (Extended Data Fig.  5h). These results indicate that KIAA1429-mediated m 6 A modi-fication of MTA elements is extremely important for maintaining a high abundance of MTA RNA in GV oocytes. As KD of MTA in GV oocytes results in GV arrest 25 , the dramatic loss of MTA RNA in Kiaa1429 Zp3 -cKO oocytes may lead to abolished oocyte competence. Potential m 6 A readers for maternal RNA stabilization. Next, we tried to identify the m 6 A readers involved in stabilizing maternal RNA. We first focused on YTHDF2, which is involved in the decay of maternal RNAs in mouse oocytes and in the zebrafish MZT process 15,20 . We used previously established Ythdf2 Zp3 -cKO mice and observed that most of their oocytes were able to develop to the MII Peak enrichment log2(obs/exp)   stage and be fertilized but were ultimately arrested at the 1C or 2C stage, which was similar to what was observed in a published study 20 (Extended Data Fig. 6a,b). Then, single-cell RNA-seq (scRNA-seq) was performed using MII oocytes and 2C embryos from Ythdf2 Zp3 control and cKO mice. Further analysis showed that the transcriptome differences mainly occurred in 2C stage but not in MII oocytes (Extended Data Fig. 6c,d). Compared with the unchanged M-decay transcripts after fertilization, the degradation of Z-decay transcripts was repressed upon YTHDF2 depletion in the 2C stage (Fig. 4i). However, to our surprise, the degradation dysfunction caused by YTHDF2 depletion was also observed on the m 6 A − Z-decay transcripts (Fig. 4i). These results indicate that YTHDF2 has limited effects on maternal mRNA stabilization but affects the degradation of Z-decay transcripts in an m 6 A-independent manner. We then assessed the RNA levels of MTA in MII and 2C samples and found that the decay of MTA was blocked in Ythdf2 Zp3 -cKO embryos, suggesting that YTHDF2 regulates the degradation of MTA in fertilized embryos (Fig. 4j). Besides YTHDF2, we wondered whether another m 6 A reader protein family, IGF2BP, was involved in the stabilization of m 6 A + transcripts in oocytes, similar to the observation in cancer cells 37 .
We assessed the expression patterns of Igf2bp mRNAs and found that Igf2bp2 and Igf2bp3 were highly expressed in oocyte samples (Extended Data Fig. 6e). We further found that both IGF2BP2 and IGF2BP3 could target MTA in oocytes through RIP-qPCR; additionally, the abundance of MTA RNA decreased modestly in IGF2BP2-depleted oocytes from Igf2bp2 KO mice (Extended Data Fig. 6f,g). Together, our data indicate that YTHDF2 may be involved in the decay of MTA after fertilization and that IGF2BP proteins may recognize MTA RNAs and help to stabilize their high abundance in oocytes.
KD of METTLs impacts the decay of m 6 A + ZGA transcripts.
Our data revealed that more than half of the ZGA transcripts and MERVL RNAs were marked with m 6 A when they were expressed after the L1C stage, which has also been observed in zebrafish but has rarely been discussed 15 Fig. 7a and Methods). Although Mettl3 and Mettl14 were highly expressed in MII oocytes, all three genes possessed considerable RNA levels in the L1C stage (Extended Data Fig. 7b). Upon injection, the RNA level of Mettl3/14/16 in morulae was reduced to less than half, and the developmental rate decreased modestly (Extended Data Fig. 7c-e). We then performed RNA-seq with the m 6 A-inhibited and control embryos at the 2C and morula stages (Extended Data Fig. 7a). We found that at the middle 2C stage the level of m 6 A + and m 6 A − ZGA transcripts was comparable, but at the morula stage, the expression of m 6 A + ZGA transcripts increased in KD embryos (Extended Data Fig. 7f). Transcriptome analysis showed that the difference was almost undetectable at the 2C stage, but a group of genes were significantly upregulated in KD embryos at the morula stage (Extended Data Fig. 7g,h and Supplementary Table 5). Interestingly, most of these upregulated genes were 2C-specific ZGA genes with m 6 A modification on their transcripts, including many 2C marker genes, such as Zscan4 family members (Extended Data Fig. 7h-j). Moreover, the RNA level of MERVL also increased significantly in KD embryos at the morula stage (Extended Data Fig. 7k,l).
After ZGA, some 2C-specific mRNAs and TEs, such as MERVL, decay through an unclear mechanism, and histone modification H3K9me3 is established at these regions at the eight-cell stage to shut down their transcriptional activity 39 . Therefore we examined the half-lives of Zsacn4 and MERVL transcripts as well as H3K9me3 modification in control and KD embryos. After the treatment of actinomycin D, the levels of Zscan4 and MERVL transcripts remained stable in KD embryos but decreased dramatically in the control group (Extended Data Fig. 7m). We then performed chromatin immunoprecipitation (ChIP)-qPCR of H3K9me3 on MERVL and Zscan4 at the morula stage, and found the installing of H3K9me3 was not obviously impacted (Extended Data Fig. 7n). All results support the idea that the upregulation of the 2C transcripts and MERVL upon KD of writers results from failure of RNA decay rather than transcriptional activation.

METTL3 ensures the timely decay of ZGA transcripts.
In the METTLs KD experiments, although we confirmed the downregulated mRNAs of Mettls, we found the protein level of METTL3 was not reduced significantly at 2C stage, and the efficiency of the antibody was difficult to evaluate (Extended Data Fig. 7o). Considering that METTL3 is an abundant maternal protein with high methyltransferase activity 38 , we used the reported inhibitor STM2457 (ref. 40 ) to fully block the function of both maternal and zygotic METTL3 proteins (Fig. 5a). When 10 μM STM2457 (ref. 40 ) was present in the culture medium, the development of embryos appeared normal at the beginning of the 2C stage (in vitro fertilization (IVF) 19 h) (Extended Data Fig. 8a) but delayed at the beginning of the 4C stage (IVF 44 h) (Fig. 5b). In addition to the developmental delay, the blastocyst formation ratio was reduced dramatically at embryonic day (E) 3.5 with STM2457 treatment, which is more serious than the phenotype of METTLs KD (Fig. 5c and Extended Data Fig. 8b,c).
To explore the function of METTL3, we tested the m 6 A level after STM2457 treatment by LC-MS/MS and ULI-MeRIP-seq (Fig. 5a). We confirmed that the m 6 A/A ratio of total RNA at the 4C stage was decreased after STM2457 treatment according to LC-MS/MS (Fig.  5d). The ULI-MeRIP-seq results showed that the m 6 A + transcript number and average m 6 A level were decreased by approximately 1/4 (Extended Data Fig. 8d,e). Detailed analysis showed that the ZGA m 6 A + transcripts, but not the maternal transcripts, lost m 6 A upon METTL3 inhibition, which indicates that METTL3 is particularly response for the de novo establishment of m 6 A on ZGA transcripts ( Fig. 5e and Extended Data Fig. 8f).
To investigate the transcriptional effects of inhibition of m 6 A, we further performed total RNA-seq at the 2C, 4C and morula stages (Fig. 5a). A large number of genes were differentially expressed between the DMSO and STM2457 treatment groups at the 4C and morula stages, but the gene expression of 2C embryos was almost unaffected (Fig. 5f, Extended Data Fig. 8g and Supplementary Table 6). In addition, only the m 6 A + ZGA genes were upregulated at the 4C and morula stages (Fig. 5g). As we found that writer KD affects the decay of 2C transcripts, we classified the ZGA genes into genes with sustained expression (sustained or upregulated after the 2C stage) and genes with transient expression (downregulated after the 2C stage) (Fig. 5h and Methods). The m 6 A + / m 6 A − gene ratio for transiently expressed genes was significantly higher than that for sustained genes (Extended Data Fig. 8h). ULI-MeRIP-seq for the 2C stage showed that the inhibition of m 6 A installing occurred on both sustained and transiently expressed m 6 A + ZGA transcripts after treatment with STM2457 ( Fig. 5i and Extended Data Fig. 8i). Overall gene expression of sustained ZGA genes changed modestly, with m 6 A + genes slightly increased after STM2457 treatment (Fig. 5j). However, for transiently expressed ZGA genes, a substantial increase after STM2457 treatment occurred in the 4C and morula stages, consistent with delayed development, and m 6 A + genes were upregulated more than m 6 A − genes, especially at the morula stage (Fig. 5j). We further confirmed the decreased m 6 A levels and defective decay of Zscan4 by qPCR (Extended Data Fig. 8j,k). Furthermore, we assessed the m 6 A and RNA dynamics on TEs via STM2457 treatment. Similarly, the ZGA transposon MERVL, but not the maternal transposon MTA, showed reduced m 6 A modification, similar to the finding for the transiently m 6 A + genes (Fig.  6a-c). In addition, the degradation of MERVL under STM2457 treatment was impaired, and the RNA level at the morula stage was increased ( Fig. 6d and Extended Data Fig. 8l,m).
Overall, the results from our experiment with the METTL3 inhibitor STM2457 strongly suggest that m 6 A modification promotes the RNA degradation of both transient ZGA genes and MERVL.

Discussion
In mouse oocytes and embryos, m 6 A is believed to play pivotal roles in these processes, as proven by the many developmental defects in different writer-or reader-KO mouse models, yet we lack detailed knowledge of the molecular mechanisms [19][20][21]41 . In this study, we developed ULI-MeRIP-seq for m 6 A mapping and revealed the highly dynamic landscape of the m 6 A methylomes on both coding and repetitive RNAs. Importantly, our study elucidates the essential and distinct roles of m 6 A during the MZT: (1) maintaining the stability of maternal Z-decay mRNAs and (2) promoting the decay of 2C-specific mRNAs (Fig. 6e). Importantly, our data on MTA and MERVL highlight the timely storage and clearance of TEs that occurs through m 6 A-mediated post-transcriptional regulation, which may also play important roles in multiple developmental events.
We found that the upregulation of MD transcripts including MTA is dramatic in YTHDF2-deficient 2C embryo; however, this upregulation seemed to be independent of m 6 A. This phenomenon indicates that YTHDF2's function in regulating coding genes during the MZT is not limited to the direct recognition of m 6 A. Previous studies have shown that binding of YTHDF2 to m 6 A-modified RNA can lead to phase separation in vitro and promote stress granule formation 42,43 . Intriguingly, a recent study reported that FMR1 regulated maternal RNA decay by m 6 A-instructed granule phase switch in Drosophila 44 . Whether a similar mechanism exists in mouse MZT process merits further investigation.
We used the METTL3 inhibitor STM2457 to inhibit m 6 A establishment. Compared with the KD assay, we observed more significant developmental defects, including a 2C exit barrier, delayed blastocyst formation and differentiation defects. This observation indicates that the m 6 A writing mediated by maternal METTL3 protein after fertilization is quite important for further development.
Here we found that most ZGA transcripts including MERVL gained m 6 A methylation in a co-transcriptional manner. Additionally, the expression level of MERVL was downregulated at the 2C stage when m 6 A was inhibited, indicating that m 6 A modification of MERVL may have other roles besides promoting decay. In a previous study based on cell lines, m 6 A was found to be co-transcriptionally incorporated by METTL3, which further impact chromatin state and transcription activity 10,27 . It will be interesting to determine whether co-transcriptional incorporation of m 6 A in ZGA is also related with the efficiently reactivation of zygotic genome, and to determine whether the dual function of m 6 A exists in different developmental systems.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41556-022-00915-x.

Methods
Our research complies with all relevant ethical regulations, and the study protocol is approved by Tongji University.
Animal and cell lines. All wild-type (WT) mice were housed in a specific-pathogen-free animal facility at Tongji University, Shanghai, China. All animal maintenance and experimental procedures were performed according to the Tongji University guide for the use of laboratory animals. All animal experiments were approved by the Biological Research Ethics Committee of Tongji University. Kiaa1429 Zp3-cKO mice were obtained from Bin Shen's laboratory 19 . Mouse strain, sex, number and age of animals in every experiment is listed in Supplementary Table 7. R1 mES cell line (ATCC SCRC-1011) was used in this study.
Culturing and collection of oocytes and embryos. For the total RNA-seq and ULI-MeRIP-seq of WT embryos, we collected embryos from C57BL/6 female mice and DBA/2 male mice (Supplementary For ULI-MeRIP-seq of WT oocytes and embryos, 500 GV oocytes, 540 MII oocytes, 750 L1C embryos, 580 L2C embryos and 300 4C embryos were collected. For Kiaa1429 control and cKO GV oocytes, 630 oocytes each were collected. For STM2457-treated embryos, 300 L2C embryos in DMSO and 300 in STM2457 were collected.
For STM2457 treatment, we used 10 μM STM2457 in CZB medium and 1:2,000-diluted DMSO in CZB medium as a control. Embryos were transferred to the medium after 4 h of IVF and cultured in the medium.
To test the half-lives of specific RNAs, we transferred 4C embryos (IVF 48 h) into CZB medium containing 1 μg ml −1 actinomycin D (Sigma, A9415). Embryos were collected at 0 h, 3 h and 6 h with actinomycin D treatment.
Detailed information including the backgrounds of the embryos and the times of collection is listed in Supplementary Table 7. ULI-MeRIP. ULI-MeRIP was performed according to a previously described protocol 28 with some modifications. The step-by-step ULI-MeRIP-seq protocol can be found at Nature Protocol Exchange 45 .
RNA extraction and purification for ULI-MeRIP. First, embryos were lysed in RNAiso Plus reagent. A 1/5 volume of chloroform was added to the samples and mixed thoroughly, and then the mixture was transferred to a MaXtract high-density tube (Qiagen, 129046) for centrifugation at full speed for 10 min. The aqueous phase was recovered and transferred to a new tube. A 1/10 volume of 3 M NaAc, 2 μl of glycogen (5 mg ml −1 , Roche, 10901393001) and an equal volume of isopropanol were added and mixed for RNA precipitation. The mixture was incubated at −80 °C for 30 min. The RNA was precipitated after centrifugation for 15 min at full speed. After washing twice with 75% ethanol, the RNA pellet was resolved in RNase-free water. Upon assessment via Qubit RNA HS Assay (Thermo Fisher Scientific, Q32852), approximately 120 ng of total RNA was extracted from 500 embryos or oocytes.
MeRIP and purification. A 1/20 volume of fragmented RNA mixture was kept as input. The remaining sonicated sample was diluted in 600 μl of IP buffer and divided into three replicates of IP reactions. Each 200 μl dilution mixture was incubated with antibody-coated beads as one IP reaction and rotated at 4 °C for 4-5 h. After IP, the samples were washed with IP buffer twice, low-salt buffer (50 mM NaCl, 10 mM Tris-HCl (pH 7.5) and 0.1% NP-40 in nuclease-free H 2 O) twice and high-salt buffer (500 mM NaCl, 10 mM Tris-HCl (pH 7.5) and 0.1% NP-40 in nuclease-free H 2 O) twice for 5 min per wash. After washing, the RNA-antibody-coated beads were resuspended in 200 μl of RNAiso Plus reagent. The RNA was extracted with the same method as that used for RNA extraction from embryos. The eluted RNA was used for reverse transcription (RT)-qPCR or library preparation.
Library preparation for ULI-MeRIP-seq and total RNA-seq. Library preparation of ULI-MeRIP samples and total RNA samples was performed using a SMARTer Stranded Total RNA-Seq Kit version 2 (Takara, 634411) according to the manufacturer's protocol. Fragmented input RNA and IP-extracted RNA was reverse transcribed into cDNA without fragmentation. Total RNA of embryos was subjected to the fragmentation step according to the manufacturer's protocol. The libraries of IP and input samples (or total RNA samples) were amplified for 15-19 cycles and 14-15 cycles, respectively. The purified libraries were sequenced on an Illumina NovaSeq 6000 platform (Novogene and Nanjing Jiangbei New Area Biophamaceutical Public Service Platform).

RT-qPCR.
For ESC samples, RNA was reverse transcribed using an All-in-One Reverse Transcription Kit (ABM, G490) directly after ULI-MeRIP. cDNA was used for qPCR with primers produced by GENEWIZ (Supplementary Table 8). For embryos with STM2457 treatment and KD efficiency tests, the RNA levels were tested by qPCR using cDNA reverse transcribed with an All-in-One Reverse Transcription Kit. For other embryo samples, the RNA level and m 6 A enrichment level were tested using diluted sequencing libraries.
SN ratio calculation. The SN ratio was calculated according to the following equation: SN ratio = 2^(−Ct m6A-positive or target gene in IP + Ct m6A-positive or target gene in input )/ 2^(−Ct m6A-negative gene in IP + Ct m6A-negative gene in input ).

Quantitative analysis of m 6 A levels.
For quantification of m 6 A levels in embryos and ESCs, total RNA of samples (about 200 embryos each stage) was injected into an LC-MS/MS instrument capable of ultra-performance liquid chromatography with a C18 column and triple-quadrupole mass spectrometer (AB SCIEX QTRAP 5500). The m 6 A level was detected in positive ion multiple reaction-monitoring mode and quantified by nucleoside-to-base ion mass transitions (282.0 to 150.1 for m 6 A and 268.0 to 136.0 for A). The m 6 A concentration was calculated from the standard curve, which was generated from pure nucleoside standards.
Immunofluorescence. GV oocytes of control and cKO mice, 2C embryos with writer KD were fixed with 4% paraformaldehyde (Sigma, P6148) for 30 min at room temperature, followed by permeabilization with 0.5% Triton X-100 in PBS for 45 min. Then oocytes were blocked overnight in blocking buffer (2.5% BSA and 0.2% Triton X-100 in PBS) at 4 °C. Oocytes and embryos were then incubated with primary antibodies (anti-KIAA1429 antibody, Atlas Antibodies, HPA031530, 1:500 dilution; recombinant anti-YTHDF2 antibody, Abcam, ab220163, 1:500 dilution; anti-METTL3 antibody, Abcam, ab195352, 1:500) for 12 h at 4 °C. After washing with PBS three times, oocytes were incubated with secondary antibodies (Thermo Scientific, A-21206, 1:500) and Hoechst 33342 (KeyGen BioTECH, KGA212-10) for 2 h. Samples were mounted onto glass slides and scanned under a confocal laser scanning microscope (LSM 700, Carl Zeiss, Germany). miniATAC-seq. ATAC-seq libraries were prepared as previously described with minor modifications 46  Single-cell mRNA-seq. Single-cell mRNA-seq for Ythdf2 Zp3 -cKO embryos was performed as previously described 47,48 . Briefly, Ythdf2 fl/fl (control) and fl/ fl;Zp3-cre (cKO) oocytes and 2C embryos were collected and washed three times in 0.5% BSA/PBS. Then, the oocytes or embryos were transferred into 4 μl of lysis buffer using a mouth pipette. Reverse transcription and cDNA library construction was performed according to the standard protocol. The cDNA library was amplified by PCR (16-19 cycles). The amplified cDNA library was fragmented using a Covaris S220 sonicator to 200-500 nt. The sequencing libraries were prepared using a KAPA HyperPlus Library Preparation kit (Roche, kk8504). Single-end 50-bp sequencing was performed on an Illumina HiSeq 2500 platform at Berry Genomics Corporation.
KD of writers in early embryos. For m 6 A writer KD, the nontarget mixture contained an IgG antibody (67 ng μl −1 , Sigma, 12-371) and siNontarget (5 μM). The KD mixture contained an anti-METTL3 antibody (28 ng μl −1 , Abcam, ab195352), an anti-METTL14 antibody (13 ng μl −1 , Sigma, HPA-038002), an anti-METTL16 antibody (67 ng μl −1 , Abcam, ab186012), siMettl3 (5 μM), siMettl14 (5 μM) and siMettl16 (5 μM). The siRNA sequences are listed in Supplementary Table 8 and were synthesized by GenePharma Company. The mixture was prepared and centrifuged for 10 min. MII oocytes from BDF1 female mice (C57BL/6 female mice crossed DBA/2 male mice) were injected with approximately 10 pl of the mixture using a Piezo-driven micromanipulator. IVF proceeded as described in method of embryo collection. ULI-native ChIP-qPCR. Five hundred cells were used per reaction. All isolated cells from morula embryos were washed in 0.5% BSA in PBS solution three times. The ULI-native ChIP procedure was performed as previously described 49 . One microliter of H3K9me3 antibody (Active Motif, 39161) was used for each IP reaction. The sequencing libraries were generated using a KAPA Hyper Prep Kit for the Illumina platform (Roche, kk8504). The amplified library cDNA was used for qPCR. The SN ratios were calculated for MERVL/SINEb1 and Zscan4/SINEb1. RNA-seq data processing. The RNA-seq data were first subjected to Trim_galore (version 0.6.4) for adaptor trimming as well as quality control with the parameters -paired -j 7-basename. The trimmed paired-end reads were then aligned to the mm9 reference genome with random chromosomes cleaned by HISAT2 (ref. 50 ) (version 2.1.0) under the parameters -p 15-dta-cufflinks-no-mixedno-discordant. Gene expression was quantified as FPKM values by Cufflinks 51 (version 2.2.1). For the downstream data analyses, the FPKM values for each gene were averaged between replicates. The RefSeq gene annotation files were downloaded from the University of California, Santa Cruz (UCSC) browser. For genes with multiple isoforms, the longest transcripts were selected. The R package DESeq2 (ref. 52 ) (version 1.26.0) was used for gene differential expression analysis. A fold change >2 and a false discovery rate (FDR) <0.05 were used as cut-offs for downregulated and upregulated genes. Genome coverage bigWig files for the UCSC genome browser were generated by the deepTools 53 (version 3.5.0) bamCoverage tool with the parameters -normalizeUsing CPM -bs 1. Genome coverage bigWig files for aggregation plots were generated by bamCoverage with the parameters -normalizeUsing RPKM -bs 50. An aggregation plot for the input signal was plotted with the computeMatrix and plotProfile functions of the deepTools package.

RIP
ULI-MeRIP data processing. The strategy for sequencing read trimming and mapping of the ULI-MeRIP data was the same as that used for RNA-seq. The m 6 A peaks in each m 6 A IP sample were identified by the MACS2 (ref. 30 ) (version 2.2.4) callpeak tool with the corresponding input sample serving as a control. MACS2 was run with the options -g mm-nomodel-keep-dup all for each replicate. For downstream analysis, a peak was kept if it was shared among all replicates. m 6 A-modified genes were defined as genes in which any exon overlapped with an m 6 A peak. Genome coverage bigWig files for the UCSC genome browser were generated by the deepTools bamCoverage tool with the parameters -normalizeUsing CPM -bs 1. Genome coverage bigWig files for aggregation plots were generated by bamCoverage with the parameters -normalizeUsing RPKM -bs 50. The aggregation plot for the input signal was plotted with the computeMatrix and plotProfile functions of the deepTools package (version 3.5.0). The Pearson correlation coefficient was calculated using the top 2,000 genes ranked by Coefficient of Variation (CV) of fold enrichment levels of m 6 A at stop codon regions (±200 bp) 54 .

Consensus motif identification within m 6 A peaks.
High-confidence peaks with fold change values greater than 4 were used to search the enriched motifs for genes or repeats. The coordinates of peak-summit regions (50 nt upstream and downstream flanking the summit) of high-confidence peaks were retrieved, and those peak-summit regions were mapped to the annotated genes or repeats to fetch the strand information with bedtools (version 2.29.0). Then, the peak-summit regions with strand information were subjected to findMotifsGenome.pl of the HOMER 55 suite (version 4.11.1) to identify motifs under the parameters -rna -len 6.
Peak enrichment analysis. The enrichment of m 6 A peaks in gene-related genomic elements and repetitive elements was calculated as the log 2 ratio of the observed counts over expected counts of m 6 A peaks overlapping with selected genomic regions. The observed counts value was calculated as the number of m 6 A peaks overlapping with related genomic regions using intersectBed of bedtools by restricting the proportion of overlapped regions to greater than 50%. The expected counts value was calculated as the number of random peaks generated using shuffleBed overlapping with related genomic regions. The radar plot for peak enrichment was plotted with the R package ggradar (version 0.2).
Gene group categorization. MD transcripts were defined as transcripts with FPKM values in the GV or MII stage >5 and fold change values (GV/2C) >4. ZGA transcripts were defined as genes with FPKM values in the GV and MII stages <1, FPKM values in the zygote and 2C stages >1, and fold change (2C/GV) values >4. Transcripts meeting the criteria of both an input FPKM >1 and overlap with an exon m 6 A peak were defined as m 6 A-modified genes. An MD transcript was defined as m 6 A + if this transcript was modified in the GV or MII stage; otherwise, it was defined as m 6 A − . Similarly, a ZGA transcript was defined as m 6 A + if it was modified in the zygote or 2C stage; otherwise, it was defined as m 6 The MD transcripts were further classified according to published criteria 33 . Specifically, a value of 1 was added to the expression level of each gene, and the value was then log 2 -transformed in the following analysis. The MD transcripts were classified by the following formulas: M-decay: expression (GV) > expression (zygote) + 1; expression (zygote) < expression (2C) + 1. Z-decay 1: expression (GV) < expression (zygote) + 1; expression (GV) > expression (zygote) − 1; expression (zygote) > expression (2C) + 1. Z-decay 2: expression (GV) > expression (zygote) + 1; expression (zygote) > expression (2C) + 1. Stable: expression (GV) < expression (zygote) + 1; expression (GV) > expression (zygote) − 1; expression (zygote) < expression (2C) + 1; expression (zygote) > expression (2C) − 1. Z-decay 1 and Z-decay 2 were combined as the Z-decay transcripts.
GO analysis. GO analysis of genes was performed using the compareCluster function of the R package clusterProfiler 56 (version 3.14.3) under the parameters fun = 'enrichGO' , pAdjustMethod = 'BH' , ont = 'BP' , OrgDb= org.Mm.eg.db, keyType = 'SYMBOL' and qvalueCutoff = 0.05. Enrichment maps were constructed using the R package Complexheatmap (version 2.2.0). The GO terms for each functional cluster are summarized with a representative term.

Analysis of ULI-MeRIP-seq data for Kiaa1429-cKO and STM2457 treatment.
To determine the m 6 A level after Kiaa1429 KO and STM2457 treatment, the total m 6 A level was calculated as the number of reads mapped to the mouse genome divided by the number of reads mapped to the m 6 A-modified spike-in (GLuc). To normalize the input and IP signal, the genome coverage bigWig files of KO or treated IP samples and input were generated by bamCoverage with the parameters -normalizeUsing RPKM -bs 50 -scaleFactor f. The scale factor f was calculated as: where T represents the library sequence depth and N represents the read number of spike-ins. For input, N is the number of reads mapping to CLuc. For IP samples, N is the number of reads mapping to GLuc. The first step of peak calling was the same as previously described, and the peaks were further filtered by calculating the IP/input signal ratio by using scaled genome coverage bigWig files. Peaks with IP/input ratios greater than 1.5 were kept for further analysis.
For defining transient and sustained transcripts, ZGA transcripts with a fold change value (morula/2C) >3 were defined as transient transcripts; otherwise, they were defined as sustained transcripts.

TE analysis.
Repeat elements were identified with RepeatMasker according to mm9 annotations downloaded directly from the UCSC browser (http:// hgdownload.cse.ucsc.edu/goldenPath/mm9/database/rmsk.txt.gz). Alignment was performed using HISAT2 against mm9 under the parameters -p 15-dta-cufflinksno-mixed -no-discordant. For multiple mapping, up to five alignments for each read were reported. We excluded repeats that overlapped with genic regions to rule out the impact of gene expression. To quantify the expression level for each locus, the RPKM value for each genome-wide nonoverlapping 50 bp window was calculated as the signal density in the sequencing data by the bamCoverage program with the parameters -normalizeUsing RPKM -bs 50, and the RPKM values were averaged for windows covering the locus and further averaged between replicates. The expression level of each family was calculated as the average level of all copies. For downstream analysis, families with expression levels greater than 5 in at least one stage were kept.
To plot the aggregation profiles, full-length MTA and MERVL annotations were applied. Full-length MTA was defined by the presence of MTA_Mm, MTA_ Mm-int and MTA_Mm, as described in a previously published paper 57 , using mergeBed, allowing at most a 20 bp gap between adjacent annotations. Full-length MERVL was defined by the presence of MT2_Mm, MERVL-int, ORR1A3-int, MERVL-int and MT2_Mm, as described in a previously published paper 58 .
ATAC-seq data analysis. ATAC-seq data were aligned by Bowtie 2 (version 2.3.5.1) (ref. 59 ) with the parameters -no-mixed-no-discordant. PCR duplicates were removed by the Picard (version 2.21.1) function MarkDuplicates with the default parameters. MTA elements were divided into five groups according to their RNA expression levels (0-1, 1-5, 5-10, 10-50 and >50). The ATAC-seq signal of each group was measured by the average RPKM value. One MTA element (chr2:98509266-98510936) contained an unusually high and short peak of read density and was excluded from analysis. To calculate ATAC signals at promoters, genome coverage bigWig files for aggregation plots were generated by bamCoverage with the parameters -normalizeUsing RPKM -bs 50. The RPKM values were averaged for windows covering the locus and further averaged between replicates.

scRNA-seq data analysis for oocytes and embryos of Ythdf2 Zp3 -cKO mice.
Ythdf2-cKO scRNA-seq data with single-end 50 bp reads were mapped by using TopHat2 (version 2.1.1) without adaptor trimming. Gene expression was quantified as the FPKM value by Cufflinks (version 2.2.1). For the downstream data analyses, the FPKM values for each gene were averaged between replicates. Considering the potential batch effects caused by different library construction pipelines, we re-defined MD transcripts with the same criteria as those used for WT samples by using Ythdf2 control scRNA-seq data. The MD transcripts were further classified into M-decay and Z-decay transcripts by applying the criteria described above.
Statistics and reproducibility. R (version 3.6.1) was used for statistical analysis. Pearson's correlation coefficient was calculated using the cor.test function with the default parameters to evaluate the reproducibility of replicates. Two-sided Wilcoxon rank sum test was used for group comparisons using the Wilcoxon test function in R. False discovery rate was used for GO enrichment analysis in clusterProfiler. The numbers of embryos or biological replicates used in this study can be found in relative figure legends and Methods. No statistical methods were used to pre-determine the sample size. No data were excluded from the analyses. All WT embryos used in the data were randomly allocated to control and treatment groups. The investigators were not blinded to group allocation during data collection, as the treatment of different groups are identical during all the experiments. The ULI-MeRIP-seq and RNA-seq was performed using three biological replicates, unless otherwise specified in the legends. The m 6 A quantification by LC-MS/MS was performed once with three technical replicates because of the material limitation. ULI-MeRIP-qPCR was performed at least twice with similar results, and a representative result is shown. RIP-qPCR and native ChIP-qPCR were performed with two biological replicates. RT-qPCR of writer KD with actinomycin D treatment was performed once owing to material limitation. RT-qPCR of gene expression in Igf2bp2 KO GV and Ythdf2 cKO GV was performed with two biological replicates. ATAC-seq for Kiaa1429 control and cKO GV nucleus were performed with at least three biological replicates, except cKO PSN type GV, which was performed with two biological replicates.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All the MeRIP-seq, RNA-seq and ATAC-seq data generated in this study are summarized in Supplementary Table 1. The sequencing data reported in this paper have been deposited in the Genome Sequence Archive under project PRJCA004536, and the accession code is CRA003985. Previously published MeRIP-seq data of 500 ng starting RNA that were re-analysed here are available under accession code GSE116002. MeRIP-seq data of mESCs with DNaseI treatment re-analysed here are available under accession code GSE145315. mm9 reference genome (https://hgdownload.soe.ucsc.edu/goldenPath/mm9/ bigZips/mm9.fa.gz) was used in all the data. Repeat elements were identified with RepeatMasker from the UCSC browser (http://hgdownload.cse.ucsc.edu/ goldenPath/mm9/database/rmsk.txt.gz). Source data are provided with this paper. All other data supporting the findings of this study are available from the corresponding author on reasonable request.

Code availability
The code is provided at https://github.com/xiaocuixu/mouse_embryo_m6a. Fig. 2 | Validation of uLi-MeRiP-seq data quality in oocytes and embryos. a, High enrichment of m 6 A in embryo IP samples tested by qPCR of GLuc versus CLuc. b, Scatterplots illustrating the distribution of PDP values in input (x-axis) and IP (y-axis) samples. The axes represent the ratio of the number of PCR duplicate reads to the total number of mapped reads. c, Scatterplots illustrating the distribution of NRF. The axes represent the ratios of the number of distinct genomic locations covered by reads to the total number of genomic locations covered by reads. d, Scatterplots illustrating the distribution of PBC1. The axes represent the ratios of the number of genomic locations covered by only one read to the number of distinct genomic locations covered by reads. e, Scatterplots illustrating the distribution of PBC2. The axes represent the ratios of the number of genomic locations covered by only one read to the number of genomic locations covered by only two reads. f, PCC between three independent IP replicates (and two replicates for 4 C) calculated by reads count. Data are presented as mean values. g, Heatmap depicting the Pearson correlation of different samples of the top 2000 transcripts ranked by CVs of fold enrichment (IP/input) levels of m 6 A at ± 200 bp around the stop codons. h, Saturation plot showing that high (FPKM > 1), medium (0.1<FPKM < 1) and low (0.01<FPKM < 1) levels of gene expression can be covered in input. i, m 6 A + transcripts detected by MACS and MeTPeak. Unique represents m 6 A + transcripts detected by either MACS or MeTPeak. Overlap represents m 6 A + transcripts shared by MACS and MeTPeak. j, Density of m 6 A peak length in oocytes, early embryos and mESCs. k, Sequence logo and p-values of the consensus motif of m 6 A peak centers in each sample. l, [RRACH] density at exon peaks in each stage. Box plots are presented with horizontal line, median; box, interquartile range (IQR); whiskers, most extreme value within ±1.5 × IQR. m, Bar chart presenting the fraction of m 6 A peaks in different genomic regions. Data in a-m represent results from three independent IP experiments of GV to 2-cell and ESC, and two independent IP experiments of 4-cell. Fig. 3 | m 6 A dynamics during MZt. a, The number of transcripts harboring m 6 A (m 6 A + ) in each stage. b, Alluvial diagram showing the global dynamics of m 6 A + transcripts during MZT(n = 6940). Each line represents a transcript that is classified as a m 6 A + class in at least one stage. c, Bar chart showing m 6 A peak number in each sample. d, Alluvial diagram showing the global dynamics of m 6 A peaks during MZT (n = 43980). Each line represents a m 6 6 As. Box plots are presented with horizontal line, median; box, interquartile range (IQR); whiskers, most extreme value within ±1.5 × IQR. Two-tailed unpaired Wilcoxon test was used to calculate the p-values. i, Pie plot showing M-decay and Z-decay transcripts ratios of MD transcripts. Data in a-i represent results from three independent IP experiments of GV to 2-cell and ESC, and two independent IP experiments of 4-cell. Peaks of each stage are the intersect peaks of replicates or triplicates.