Intro
Life begins with fertilization, a transformative process in which terminally differentiated gametes are reprogrammed into totipotent embryos through the oocyte-to-embryo transition (OET) [ 1 ]. After entering meiosis, primary oocytes are arrested at prophase I and remain in this state until gonadotropin stimulation prompts full-grown oocytes (FGOs) to resume meiosis. Following this, oocytes are arrested again at metaphase II until fertilization occurs. During the period from FGOs to fertilized embryos, before zygotic genome activation (ZGA), the genome remains transcriptionally quiescent. As a result, post-transcriptional regulatory mechanisms predominantly control oocyte maturation and the OET [ 2 , 3 ]. These mechanisms include RNA modification networks that regulate pre-deposited maternal RNAs in oocytes, particularly internal N 6 -methyladenosine (m 6 A) and N 6 ,2′-O-dimethyladenosine (m 6 Am) modifications, the latter occurring at the first transcribed nucleotide of messenger RNA (mRNA) transcripts [ 4 , 5 ].
m 6 A, the most prevalent internal RNA modification on RNA polymerase II (Pol II) transcripts, is dynamically regulated by ‘writer’ proteins (methyltransferases) and ‘eraser’ proteins (demethylases), which control RNA fate through interactions with ‘reader’ proteins [ 6 ]. This modification is considered a master regulator of RNA metabolism, influencing mRNA transcription and processing, stability, translation, and localization [ 7 ]. Dysregulation of m 6 A has been mechanistically linked to various human diseases, from metabolic disorders like obesity to malignancies such as cancer [ 8 , 9 ]. In contrast, m 6 Am is a distinct modification localized to the adenosine nucleotide of the 5ʹ terminal transcription start site (TSS) adjacent to the mRNA cap structure. Its biosynthesis involves two sequential steps: (i) 2ʹ-O-methylation of the ribose by CMTR1, a modification conserved across most eukaryotes, followed by (ii) N 6 -methylation of the adenine base catalyzed by PCIF1 protein [ 10 ]. Methylation by PCIF1 seems efficient in human cell lines, as up to 48%–92% of m 7 G-capped RNAs with a TSS adenosine carry the m 6 Am mark, with the rest being Am [ 11 – 13 ]. Functional studies reveal that m 6 Am exhibits context-dependent roles, enhancing RNA stability in some systems while promoting or repressing cap-dependent translation depending on cellular conditions or experimental models [ 11 , 13 – 16 ].
Previous studies have explored the m 6 A methylome in mouse oocytes and early-stage embryos using various m 6 A detection strategies, including ULI-MeRIP-seq, pico-MeRIP, scm 6 A-seq, and SLIM-seq [ 17 – 21 ]. These studies have shown that m 6 A is widely deposited on maternally inherited transcripts, ZGA gene transcripts, key transcriptional regulators, and retrotransposon RNAs such as MTA and MERVL. Functionally, m 6 A plays a dual role in ensuring the stability of mRNAs in oocytes and facilitating the timely decay of two-cell-specific transcripts [ 17 ]. Although the dynamics and functional roles of m 6 A during human early development have been reported recently [ 22 ], the transcriptome-wide profiles of m 6 Am in both human and mouse oocytes and early embryos have yet to be characterized. In this study, we developed a method, MeLACE-seq [methylated RNA immunoprecipitation with linear amplification of complementary DNA (cDNA) ends and sequencing], which combines an optimized LACE-seq protocol [ 23 ] with ultraviolet (UV)-crosslinking of the m 6 A antibody to both m 6 A and m 6 Am sites. This approach enables simultaneous profiling of m 6 A and m 6 Am methylomes with high resolution and low cell input requirements. Using this approach, we profiled the dynamic m 6 A and m 6 Am methylomes across multiple developmental stages of human and mouse oocytes and early embryos. Our study comprehensively analyzes the conserved regulatory patterns of m 6 A and m 6 Am during the OET in both human and mouse species.
Results
To investigate the transcriptome-wide distribution and dynamics of m 6 A and m 6 Am in human and mouse oocytes and early embryos, we developed MeLACE-seq based on a similar strategy of LACE-seq that we recently reported [ 23 ], which features several key steps (Fig. 1A ). Briefly, (i) we crosslinked a N 6 -methyladenine-specific antibody to m 6 A and m 6 Am sites on RNA transcripts within cell lysates using UVC light; (ii) we ligated a 3′ adaptor directly to the 3′ end of reverse-transcribed RNA fragments by substituting the step of adding polyA tail to cDNA within one-tube reaction to convenient to paired-end sequencing; (iii) during linear amplification of cDNA fragments via in vitro transcription, we selectively depleted crRNA fragments to minimize ribosomal RNA interference and increase the percentages of other m 6 A and m 6 Am marked transcripts.
Genome-wide mapping of m 6 A and m 6 Am methylome in human and mouse oocytes and early embryos. (A) Schematic illustration of the MeLACE-seq method. After cell lysis, the m 6 A antibody-m 6 A RNA complexes were purified after UV-crosslinking and RNA fragmentation, followed by in vitro transcription, library preparation (see the ‘Materials and methods’ section), and sequencing. (B) Schematic of human and mouse oocytes and embryos used for m 6 A and m 6 Am methylome (MeLACE-seq data) comparison. (C) The distribution of the enriched cap m 6 Am and m 6 A signals derived from MeLACE-seq datasets in human and mouse oocytes and embryos. (D) Principal component analysis of cap m 6 Am and m 6 A methylomes derived from MeLACE-seq datasets for human and mouse oocytes and embryo stages. (E) Alluvial plot showing the global dynamics of cap m 6 Am + and m 6 A + genes during human and mouse OETs. Each line represents a transcript that is classified as a m 6 A + or m 6 Am + group in at least one stage. (F) Genes gaining and losing m 6 A/m 6 Am marks correspond to a change in gene expression between consecutive oocyte and embryo stages. The transcriptome datasets of human and mouse oocytes and early embryos are derived from previous publications [ 31 , 32 ]. Up, the log 2 of the expression fold change of one stage versus the next is > 1; Down, <−1; Unchanged, otherwise.
Given that Mettl3 serves as the primary m 6 A writer in mammals and cKO of Mettl3 results in significantly reduced ovary size and impaired oocyte maturation [ 24 ], we conducted MeLACE-seq on both control ( Mettl3 flox/flox ) and Mettl3 cKO ( Mettl3 flox/flox ; Zp3 -Cre) oocytes at the GV stage to assess the reproducibility and specificity of our method as a proof-of-concept. Our analysis revealed high reproducibility between biological replicates ( R = 0.93–0.95) and a dramatic reduction in m 6 A signals in Mettl3 cKO oocytes compared to controls ( Supplementary Fig. S1A and B ), demonstrating the specificity of MeLACE-seq in enriching m 6 A signals. Since the captured signals in MeLACE-seq rely on an antibody targeting N 6 -methyladenine, which binds to both m 6 A and m 6 Am sites, we further analyzed m 6 Am peaks by distinguishing them from m 6 A peaks based on their localization at TSS, as previously described [ 30 , 36 ] (see the ‘Materials and methods’ section). Then, we examined the internal m 6 A and cap m 6 Am signals in Mettl3 cKO oocytes separately. The results showed that internal m 6 A signals in gene bodies were dramatically reduced, whereas cap m 6 Am signals around TSS regions were not markedly altered in Mettl3 cKO oocytes ( Supplementary Fig. S1C ). These findings indicate that depletion of the m 6 A writer Mettl3 selectively impairs m 6 A deposition without affecting m 6 Am deposition. Consistent with established findings [ 37 ], the internal m 6 A peaks were predominantly enriched within the conserved DRACH motif (D = G/A/U, R = G/A, A = m 6 A, H = A/C/U) ( Supplementary Fig. S1D ), further validating the specificity of our MeLACE-seq approach. The first transcribed m 6 Am peaks were enriched within the canonical BCA motif (A = m 6 Am; B = C, G, or U) ( Supplementary Fig. S1D ), consistent with prior studies [ 37 ]. To further probe the performance of MeLACE-seq, we compared our data with existing datasets from MeRIP-seq [ 24 ], ULI-MeRIP-seq [ 17 ], and picoMeRIP-seq [ 18 ] assays performed on varying numbers of mouse GV oocytes ( Supplementary Table S3 ). Despite using fewer oocytes as input, MeLACE-seq detected a comparable number of m 6 A-marked genes and a greater number of m 6 Am-marked genes and achieved higher resolution than these MeRIP-based methods ( Supplementary Fig. S1E–J ). It is worth noting that the MeLACE-seq method preserves the 5′ end sequence information of RNA through template switching during reverse transcription, enabling the reliable identification of cap-specific m6Am sites. Additionally, truncated cDNA ends align precisely downstream of antibody- N 6 -methyladenine crosslinking sites, thereby enhancing the resolution of putative m 6 Am and internal m 6 A modifications.
To further validate the reliability of MeLACE-seq, we employed SELECT [ 37 ] to verify m 6 A sites identified from MeLACE-seq signals ( Supplementary Fig. S2A ). The results showed that putative m 6 A sites detected by MeLACE-seq could be efficiently validated using SELECT, as exemplified by MeLACE-seq signals on Btg4, Foxo1, Gdf9 , and Obox1 mRNAs ( Supplementary Fig. S2B–E ). To assess the reliability of m 6 Am signals identified by MeLACE-seq, we performed MeLACE-seq in HEK293 cells following knockdown of the m 6 Am writer PCIF1 via transfection with PCIF1 -targeted siRNA ( Supplementary Fig. S2F and G ). Upon efficient PCIF1 depletion, m 6 Am signals around TSS regions were dramatically reduced, whereas m 6 A signals within gene bodies were not markedly altered ( Supplementary Fig. S2H and I ). These findings demonstrate the specificity and reliability of m 6 Am sites identified by MeLACE-seq. We further compared m 6 A-containing genes identified by MeLACE-seq with those reported in published GLORI 3.0 and miCLIP datasets [ 37 , 38 ]. Notably, 87.9% and 93.1% of m 6 A-containing genes detected by MeLACE-seq overlapped with those identified in the GLORI 3.0 and miCLIP datasets, respectively ( Supplementary Fig. S2J and K ), indicating a high level of confidence in MeLACE-seq-derived m 6 A calls. In summary, these results demonstrate that MeLACE-seq is a highly sensitive and robust approach for the simultaneous profiling of m 6 A and m 6 Am modifications.
Next, we applied MeLACE-seq to profile m 6 A and m 6 Am modifications in human and mouse oocytes and pre-implantation embryos. For human samples, we analyzed oocytes at GV, metaphase I (MI), and metaphase II (MII) stages, as well as early embryos at the zygote (1C), four-cell (4C), eight-cell (8C), and blastocyst (BL) stages. For mouse samples, we included oocytes at the GV stage and early embryos at the zygote, early two-cell (E2C), late two-cell (L2C), 8C, morula (MO), and blastocyst stages (Fig. 1B and Supplementary Table S1 ). Each stage was analyzed with two biological replicates, which exhibited a high correlation across all developmental stages in both human and mouse species ( Supplementary Fig. S3A and B ). As expected, both human and mouse samples significantly enriched the m 6 A and m 6 Am signals near stop codons and TSSs across gene bodies (Fig. 1C ), with substantial signals also observed in exon regions ( Supplementary Fig. S3C ). The distribution patterns of these modifications were highly consistent across all stages of oocytes and embryos in both species. Specifically, internal m 6 A peaks were predominantly localized within the conserved DRACH motif, and cap m 6 Am peaks were enriched within the canonical BCA motif ( Supplementary Fig. S3D and E ), aligning with previously reported in vitro RNA m 6 A and m 6 Am enrichment motifs [ 36 , 37 ].
We identified a total of 15 560–75 811 m 6 A peaks and 4203–15 443 cap m 6 Am peaks across seven developmental stages in human samples, corresponding to 5311–12 494 genes and 2055–4908 genes, respectively. Similarly, for mouse samples, we detected 28 625–62 857 m 6 A peaks and 1894–10 978 cap m 6 Am peaks across seven developmental stages, mapping to 7691–9760 genes and 1089–3692 genes, respectively. Taken together, these results highlight the dramatic dynamic changes in both m 6 A and m 6 Am methylomes during human and mouse OETs.
We observed that the m 6 A and m 6 Am methylomes exhibited dynamic changes from mouse GV oocytes to blastocysts, clustering into three distinct groups in chronological order ( Supplementary Fig. S3F ). A significant transition occurred during the E2C and L2C stages, coinciding with the mouse ZGA. In contrast, the m 6 A and m 6 Am methylome in human GV oocytes exhibited a distinctive distribution ( Supplementary Fig. S3F ), clustering with the 4C and 8C samples, coinciding with the onset of human ZGA [ 39 ]. Notably, the cap m 6 Am methylome in both species followed a concordant chronological order from oocytes to blastocysts (Fig. 1D ). However, the human m 6 A methylome showed distinct clustering patterns, especially for GV oocytes (Fig. 1D ), underscoring species-specific differences in methylome dynamics.
To gain deeper insights into m 6 A and m 6 Am methylome dynamics, we analyzed the m 6 A- and m 6 Am-marked transcripts and peaks at each developmental stage. In human samples, the number of m 6 Am- and m 6 A-marked transcripts and peaks showed two waves of growth, with the second low point occurring around ZGA (Fig. 1E and Supplementary Fig. S3G ). However, the highest abundance of m 6 Am-marked transcripts and peaks occurred at the MII stage, while those of m 6 A-marked transcripts and peaks exhibited the highest abundance at the zygote stage. In contrast, in mouse samples, the number of m 6 Am-marked transcripts and peaks gradually declined until the L2C stage, after which they increased notably (Fig. 1E and Supplementary Fig. S3H ). Meanwhile, the number of m 6 A-marked transcripts and peaks remained relatively stable throughout early development, except for a pronounced decrease at the L2C stage. These findings reveal distinct dynamics of the m6A and m6Am methylomes during human and mouse OETs, identifying the ZGA stage as a critical regulatory node at which m 6 A and m 6 Am modification landscapes undergo dramatic, species-specific remodeling.
Then, we investigated whether expression changes were associated with the gain or loss of m 6 A/m 6 Am-marked transcripts across consecutive stages of human and mouse oocyte maturation and early embryonic development by integrating MeLACE-seq datasets with previously reported transcriptome datasets [ 31 , 32 ]. The analysis revealed that >79.6% of m 6 A- or m 6 Am-marked transcripts in human oocytes and early embryos remained unchanged between consecutive stages despite dynamic gains or losses of m 6 A/m 6 Am marks (Fig. 1F ). Similarly, in mouse oocytes and early embryos, the majority (>52.8%) of m 6 A- or m 6 Am-marked transcripts exhibited no significant expression changes across consecutive stages, despite pronounced alterations in m 6 A/m 6 Am modification (Fig. 1F ). These findings indicate that the gain or loss of m 6 A/m 6 Am peaks on specific mRNAs is largely independent of changes in transcript abundance, although a subset of transcripts shows a correlation between methylation dynamics and expression levels at specific developmental stages. Notably, this observation helps explain the apparent paradox that m 6 Am+/m 6 A + transcripts progressively increase during human oocyte maturation, while overall transcript abundance declines (Fig. 1E ).
Next, we identified differentially m 6 A- and m 6 Am-marked genes between each pair of consecutive developmental stages. This analysis revealed stage-specific marker genes and functional enrichments associated with key processes in human early development ( Supplementary Fig. S4A and B ). For instance, m 6 Am gain genes during the MII-to-1C transition were significantly enriched in functions related to chromosome segregation, including genes such as INO80, AXIN2 , and SEH1L ( Supplementary Fig. S4A ). In contrast, genes with m 6 A gain at multiple transitions, such as MRPL18, OGDH, RPUSD3, MTRF1L, METTL4, PUS1 , and HSD17B10 , were associated with mitochondrial activity ( Supplementary Fig. S4B ). These genes play crucial roles in modulating Ca 2+ signaling, ATP production, reactive oxygen species, and intermediary metabolites during early development [ 40 ].
Based on the dynamics of the cap m 6 Am and m 6 A methylomes during OET, we classified genes into four major groups: (i) ‘Oocyte-specific’ genes are marked by m 6 Am or m 6 A at the GV stage but lose these marks during OET and remain unmarked after ZGA (Fig. 2A and B ). Examples of m 6 Am-marked genes in this group include F-box genes [ 41 ], part of the E3 ligase complex, such as FBXO25 and Fbxo43 in both human and mouse species (Fig. 2C ). m 6 A-marked genes in this group include those involved in mRNA degradation and oocyte development, like BTG4 and Gdf9 (Fig. 2D ). (ii) ‘OET gain’ genes are initially unmarked by m 6 Am and m 6 A at the GV stage but gain these marks during OET. Genes involved in chromosome segregation and stage-specific regulation, such as NUF2 and Zscan4f , belong to this group (Fig. 2C ). (iii) ‘OET loss’ genes are marked by m 6 Am or m 6 A at the GV stage, lose these marks during OET, and are gained later. This group includes genes involved in transcription regulation (e.g. MAX ), transcription (e.g. Polr2e ), mRNA processing (e.g. WDR33 ), and microtubule motor activity (e.g. Kif1c ). (iv) ‘Embryonic’ genes are marked by m 6 Am or m 6 A after ZGA and are involved in transcription and embryonic development, such as FGFR4, Tmem159, GATA6 , and Cdx2 ( Supplementary Fig. S5A–E and Supplementary Table S4 ).
Global dynamics of m 6 A and m 6 Am methylome in human and mouse oocytes and early embryos. (A) Heat maps showing the clustering results based on the Cap m 6 Am status for all expressed RNA with m 6 Am status change from human and mouse consecutive stage transitions. The number of each cluster gene is listed. (B) Heat maps showing the clustering results based on the Cap m 6 A status for all expressed RNA with m 6 A status change from human and mouse consecutive stage transitions. The number of each cluster gene is listed. (C) Integrative Genomics Viewer (IGV) Genome browser views of representative m 6 Am-marked gene transcripts from human and mouse showing different clustering patterns. (D) IGV Genome browser views of representative m 6 A-marked gene transcripts from human and mouse showing different clustering patterns. (E) Box plots showing the expression levels and TE of m 6 Am marked ‘OET gain’ and ‘OET loss’ gene transcripts during human and mouse OETs. The transcriptome and translatome datasets of human and mouse oocytes and early embryos are derived from previous publications [ 31 , 32 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × interquartile range (IQR). (F) Box plots showing the expression levels and TE of m 6 A-marked ‘OET gain’ and ‘OET loss’ gene transcripts during human and mouse OETs. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. P -values in panels (E) and (F) were determined by the Kolmogorov–Smirnov test.
To investigate whether the dynamics of m 6 A and m 6 Am methylomes are coupled with changes in the transcriptome and translatome during the OET, we integrated our m 6 A and m 6 Am methylome data with published transcriptome and translatome datasets from human and mouse OET [ 31 , 32 ]. These published R2-lite datasets revealed the global transcriptomic and translational dynamics in human and mouse oocytes and early embryos. By comparing the expression levels and TE of m 6 A- or m 6 Am-marked ‘OET gain’ genes with those of ‘OET loss’ genes across developmental stages, we observed that the expression level and TE of m 6 Am-marked ‘OET loss’ genes were generally higher than those of ‘OET gain’ genes at each stage of human and mouse oocytes and embryos, except for the inner cell mass (ICM) isolated from mouse blastocysts (Fig. 2E ). Similarly, the expression level and TE of m 6 A-marked ‘OET loss’ genes were consistently higher than those of ‘OET gain’ genes across all stages in both species (Fig. 2F ). Collectively, these findings suggest that the dynamic m 6 A and m 6 Am methylomes potentially play a regulatory role in modulating the expression and translation of their target genes during OET.
To assess the evolutionary conservation of m 6 Am and m 6 A methylomes between human and mouse OET, we compared the dynamics of these modifications in homologous genes. For m 6 Am-marked genes, we found that ‘OET gain’ genes accounted for a majority of the analyzed human genes, which are primarily corresponding to the ‘OET gain’ and ‘OET loss’ groups in the mouse (Fig. 3A and Supplementary Fig. S6A ). Similarly, for m 6 A-marked genes, a large proportion of them belongs to ‘OET gain’ group, dominantly corresponding to the ‘OET gain’ and ‘OET loss’ groups in the mouse (Fig. 3B and Supplementary Fig. S6B ). Thus, the corresponding combination of human ‘OET gain’—mouse ‘OET loss’ genes marked by both m 6 A and m 6 Am constituted a great number of the methylation dynamic specificity for the two species, and the ‘OET gain’ genes marked by both m 6 A and m 6 Am exhibited a large number for methylation dynamic conservation across species compared to the other three groups, highlighting their potential importance in gain or loss m 6 A/m 6 Am modification during OET.
Conservation and divergence of m 6 Am and m 6 A methylomes between human and mouse oocytes and early embryos. (A) Alluvial plot showing the divergence of Cap m 6 Am-marked transcripts between human–mouse homologous genes. The representative genes are listed. (B) Alluvial plot showing the divergence of m 6 A-marked transcripts between human–mouse homologous genes. The representative genes are listed. (C) Alluvial diagrams depicting four groups of m 6 Am-marked human homologous genes, corresponding to their mouse homologs, with continuous m 6 Am marking or lack thereof during mouse OET (left), or vice versa (right). The representative genes are listed. (D) Alluvial diagrams depicting four groups of m 6 A-marked human homologous genes, corresponding to their mouse homologs, with continuous m 6 A marking or lack thereof during mouse OET (left), or vice versa (right). The representative genes are listed. (E) Alluvial diagrams depicting four groups of m 6 Am-marked genes in human but absent in mouse (left), or vice versa (right). The representative genes are listed. (F) Alluvial diagrams depicting four groups of m 6 A-marked genes in human but absent in mouse (left), or vice versa (right). The representative genes are listed. (G) Box plots showing the expression levels and TE of m 6 Am marked ‘OET gain’ and ‘OET loss’ gene transcripts from human–mouse homologous genes during human and mouse OETs. The transcriptome and translatome datasets of human and mouse oocytes and early embryos are derived from previous publications [ 31 , 32 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. (H) Box plots showing the expression levels and TE of m6A-marked ‘OET gain’ and ‘OET loss’ gene transcripts from human–mouse homologous genes during human and mouse OETs. Boxes represent the 25th−75th percentile (line at the median), with whiskers at 1.5 × IQR. P -values in panels (G) and (H) were determined by the Kolmogorov–Smirnov test.
In addition to the four groups of m 6 Am and m 6 A methylome dynamics in human homologous genes that aligned with corresponding groups in mice, we observed that a substantial proportion of human homologous genes from these groups were associated with mouse homologous genes that either retained or lacked methylation across all developmental stages. Specifically, the majority of these human genes (3050 out of 3176) lacked of m 6 Am modification throughout mouse OET, with a substantial portion (2022 out of 3050) categorized as human ‘OET gain’ genes (Fig. 3C ). Conversely, a similarly high percentage of mouse homologous genes (1767 out of 2020) continuously lacked of m 6 Am mark during human OET, with the majority classified as mouse ‘OET gain’ genes (722 out of 1767) (Fig. 3C ). For the m 6 A methylome, over half of the human homologous genes (1993 out of 3878) were continuously marked by m 6 A during mouse OET (Fig. 3D ). In contrast, a higher percentage of mouse homologous genes (1633 out of 2328) were continuously unmarked by m 6 A during human OET, with mouse ‘OET gain’ genes (895 out of 1633) constituting the largest subset (Fig. 3D ). Furthermore, we analyzed the m 6 Am and m 6 A methylomes of species-specific genes in human and mouse OETs. We found that the ‘OET gain’ genes marked by m 6 Am and m 6 A contributed predominantly to the species-specific methylome in both species (Fig. 3E and F ). Collectively, these results suggest that the m 6 Am and m 6 A methylomes associated with ‘OET gain’ genes exhibit less evolutionary conservation, underscoring the species-specific regulatory mechanisms during early development.
We next investigated whether the conserved m 6 A and m 6 Am methylomes in homologous genes are coupled with the expression level and TE of ‘OET gain’ and ‘OET loss’ conserved genes. Our results showed that the conserved m 6 A methylome in homologous genes generally maintains their relevance to the expression level and TE of ‘OET loss’ genes. However, the conserved m 6 Am dynamics during human OET do not significantly correlate with the TE of ‘OET loss’ genes compared to ‘OET gain’ genes (Fig. 3G and H ), suggesting that m 6 Am modifications may specifically enhance translation for a subset of these genes.
We classified transcripts into four groups at each developmental stage to investigate the individual effects of cap m 6 Am and internal m 6 A on mRNA expression and translation during human and mouse OETs. The first group consisted of transcripts marked solely by m 6 Am, the second group included transcripts marked exclusively by m 6 A, the third group contained transcripts marked by both modifications, and the remaining transcripts were used as the control group for comparison. Based on this classification, we found that the global RNA levels and TE of transcripts containing m 6 Am, m 6 A, or both modifications tended to be higher than those of the unmarked counterparts during both human and mouse OETs (Fig. 4A and B and Supplementary Fig. S7A ). Notably, transcripts marked with m 6 Am showed a more pronounced effect on expression and translation compared to those marked only with m 6 A. These results suggest that m 6 Am and m 6 A methylations are strongly associated with elevated expression levels and TE in human and mouse OET transcripts.
Cap m 6 Am and m 6 A relate to mRNA expression level and TE during human and mouse OETs. (A) Box plots showing the expression levels (left) and TE (right) of only m 6 A, m 6 Am, and both m 6 A and m 6 Am-marked transcripts during human OET. The transcriptome and translatome datasets of human oocytes and early embryos are derived from a previous publication [ 31 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. (B) Box plots showing the expression levels (left) and TE (right) of only m 6 A, m 6 Am, and both m 6 A and m 6 Am-marked transcripts during mouse OET. The transcriptome and translatome datasets of mouse oocytes and early embryos are derived from a previous publication [ 32 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. (C) Line charts showing maternally decayed and ZGA gene number of only m 6 A-marked, only cap m 6 Am-marked, both m 6 A and cap m 6 Am-marked, and the rest of unmarked transcripts across human and mouse OETs. (D) Box plots showing the expression levels and TE of maternally decayed RNAs with only m 6 A, m 6 Am, or both m 6 A and m 6 Am-marked during human OET. The transcriptome and translatome datasets of human oocytes and early embryos are derived from a previous publication [ 31 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. (E) Box plots showing the expression levels and TE of maternally decayed RNAs with only m 6 A, m 6 Am, or both m 6 A and m 6 Am-marked during mouse OET. The transcriptome and translatome datasets of mouse oocytes and early embryos are derived from a previous publication [ 32 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. (F) Box plots showing the expression levels and TE of ZGA RNAs with only m 6 A, m 6 Am, or both m 6 A and m 6 Am-marked during human OET. The transcriptome and translatome datasets of human oocytes and early embryos are derived from a previous publication [ 31 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. (G) , Box plots showing the expression levels and TE of ZGA RNAs with only m 6 A, m 6 Am, or both m 6 A and m 6 Am-marked during mouse OET. The transcriptome and translatome datasets of mouse oocytes and early embryos are derived from a previous publication [ 32 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. P -values in panels (A), (B), (D), (E), (F), and (G) were determined by the Kolmogorov–Smirnov test. (H) The dynamics of m 6 Am and m 6 A modifications of representative master regulators for lineage segregation in human and mouse oocytes and early embryos.
The OET marks the first major developmental shift in both human and mouse species, during which maternal mRNAs are degraded, and zygotic transcription is initiated. We sought to explore the individual roles of m 6 Am and m 6 A in affecting the expression and translation of maternally decayed RNAs and ZGA transcripts. For maternally decayed RNA transcripts, we observed that in human samples, the RNA expression levels and TE of transcripts marked only by m 6 Am were generally higher than those of unmarked transcripts before ZGA (Fig. 4C – E ). In contrast, the RNA expression levels of transcripts marked only by m 6 A showed no significant difference compared to unmarked transcripts, although their TE was slightly elevated during oocyte maturation. In mouse samples, the RNA expression levels of transcripts marked only by m 6 Am or both m 6 Am and m 6 A were consistently higher than those of unmarked transcripts during OET. However, the TE of these transcripts showed no significant difference compared to unmarked transcripts (Fig. 4D and E ).
Given that previous studies have shown that maternally encoded transcript decay occurs either before or after ZGA [ 18 , 42 ], we categorized these maternally decayed RNAs into three subsets based on their degradation patterns: M-decay, Z-decay, and C-decay ( Supplementary Fig. S7B ). We further investigated the individual effects of m 6 Am and m 6 A on the expression and translation of these subclasses of maternally decayed RNAs. Our results revealed that the impact of m 6 Am marking on TE of human maternally decayed RNAs is primarily attributable to Z-decay transcripts, rather than M-decay or C-decay transcripts, as m 6 Am-marked Z-decay transcripts generally exhibited higher expression levels and TE compared to unmarked Z-decay transcripts. While all three groups of maternally decayed RNAs marked by m 6 Am apparently affect the expression level of those transcripts moderately ( Supplementary Fig. S7C and D ). In contrast, mouse m 6 Am marks affect expression levels of all three groups of maternally decayed RNAs, while mouse m 6 A marks preferentially affect Z-decay transcripts at special stages, consistent with previous reports that m 6 A marks preferentially act on Z-decay transcripts with elevated expression levels [ 17 ]. Of note, neither modification significantly impacts the TE of any group of maternally decayed RNAs, except that m 6 A-marked C-decay transcripts exhibit reduced TE at the L2C stage ( Supplementary Fig. S7E ).
We also investigated the relevance of m 6 Am and m 6 A to the expression and translation of ZGA transcripts (Fig. 4C , F , and G). In human samples, the RNA expression levels and TE of transcripts marked only by m 6 A were generally higher than those of unmarked transcripts in oocytes and early embryos. In contrast, the RNA expression levels of transcripts marked only by m 6 Am showed no significant difference compared to unmarked transcripts around ZGA (4C–8C stages), although their TE was slightly elevated. In mouse samples, the RNA expression levels of transcripts marked only by m 6 Am, only by m 6 A, or both m 6 Am and m 6 A were consistently higher than those of unmarked transcripts around ZGA. Additionally, the TE of these transcripts was higher than that of unmarked transcripts at the L2C stage. Given that ZGA encompasses minor ZGA and major ZGA, we further investigated the individual effects of m 6 Am and m 6 A on the expression and translation of minor ZGA and major ZGA transcripts ( Supplementary Fig. S7B ). Our results showed that the minor ZGA transcripts in humans with m 6 A modification correlate with higher RNA levels and TE, while human major ZGA transcripts with m 6 A modification correlate with higher TE. In mouse early embryos, both m 6 A and m 6 Am-marked major ZGA transcripts correlate with higher RNA levels and TE ( Supplementary Fig. S8A and B ).
Considering the essential roles of master transcription regulators in ZGA and embryonic development [ 43 ] and their extensive modification by m 6 A (∼89% of expressed transcription factors are marked by m 6 A at least one developmental stage in mouse oocytes and early embryos [ 18 ]), we examined the status of m 6 A and m 6 Am modifications on transcription factor mRNAs in both human and mouse oocytes and early embryos. On average, around 70% of expressed transcription factor mRNAs were marked by m 6 A at each human developmental stage, compared to ∼79% in mice ( Supplementary Fig. S8C ). By contrast, ∼35% of human and 25% of mouse transcription factor mRNAs were marked by m 6 Am at each developmental stage. In human embryos, development proceeds into three distinct, concurrently forming lineages: epiblast (EPI), primitive endoderm (PrE), and trophectoderm (TE) [ 44 ]. The lineage specification marker genes, such as trophectoderm markers ( TEAD1/3 ), epiblast markers ( POU5F1/SOX2 ), and primitive endoderm markers ( GATA4/6 ), were widely marked by m 6 A at least one developmental stage but rarely marked by m 6 Am across all stages (Fig. 4H ). During mouse preimplantation embryo development, three distinct cell lineages are formed by a two-step process whereby outer TE cells are first segregated from ICM, followed by ICM refinement into either the PrE or EPI [ 45 ]. Similarly, in mouse embryos, marker genes for the first and second lineage specification events were extensively marked by m 6 A, consistent with previous reports [ 18 ]. However, marker genes for the primitive endoderm in the second lineage specification, such as Gata4/6, Pdgfra , and Sox7/17 , lacked m 6 Am marks across all developmental stages (Fig. 4H ). In sum, these findings suggest that, compared to m 6 A marks, m 6 Am marks exhibit greater variability and may play a more nuanced role in regulating the expression of lineage-specific transcription factors during early development.
Approximately half of the mammalian genome is derived from transposon elements, including three major classes of retrotransposons: long terminal repeat (LTR) elements, long interspersed nuclear elements (LINEs), and short interspersed nuclear elements (SINEs) [ 46 ]. These retrotransposons are widely recognized as key drivers of genome evolution due to their ability to rewire gene regulatory networks [ 47 ]. Some evolutionarily young retrotransposons are expressed and marked by m 6 A modifications in oocytes and early embryos [ 17 , 18 ]. These m 6 A-marked retrotransposons exhibit diverse regulatory functions in mouse embryonic stem cells (ESCs) and early embryos [ 48 – 51 ]. However, the role of m 6 Am-marked retrotransposons remains less explored, despite the potential for RNA polymerase II-transcribed retrotransposons like LINEs and LTRs to bear m 6 Am marks. To address this, we investigated the dynamics of m 6 A and m 6 Am marks on retrotransposons across human and mouse development, from oocyte to early embryo. Overall, we found that LINE1 RNAs were heavily marked by m 6 A around the ZGA in both human and mouse embryos, while MaLR (Multiple Active LTRs) RNAs were notably m 6 A-marked before ZGA in both species, in line with their expression patterns during OET. The mouse 2C-specific retrotransposon ERVL RNAs were heavily marked by m 6 A at the 2C stage, consistent with previous reports [ 18 ]. In contrast, human ERVL RNAs exhibited strong m 6 A enrichment prior to the ZGA stage (Fig. 5A – D ). m 6 A enrichment and expression of Alu RNAs, which belong to SINE retrotransposons, were most prominent at the human blastocyst stage, while mouse SINE retrotransposons, such as B1, B2, and B4, were heavily m 6 A-marked at the morula stage (Fig. 5A – D and Supplementary Fig. S9A and B ). Additionally, we identified m 6 Am marks at the TSS regions of retrotransposon loci, similar to what we observed in mRNA transcripts ( Supplementary Fig. S9C and D ). However, we found that m 6 Am marks are sparse on the main families of LTRs and LINEs, with <6% of each retrotransposon family bearing m 6 Am marks in human and mouse oocytes and embryos ( Supplementary Fig. S9E and F ).
Cap m 6 Am and m 6 A relate to retrotransposon expression during human and mouse OETs. (A) Heatmap plots showing the dynamics of m 6 A intensity (m 6 A signal reads/genomic copies) and expression values (RPKM > 0) for human representative retrotransposon subfamilies. The transcriptome datasets of human oocytes and early embryos are derived from a previous publication [ 31 ]. (B) IGV Genome browser views of m 6 A-marked retrotransposon subfamilies during human OET. (C) Heatmap plots showing the dynamics of m 6 A intensity (m 6 A signal reads/genomic copies) and expression values (RPKM > 0) for mouse representative retrotransposon subfamilies. The transcriptome datasets of mouse oocytes and early embryos are derived from a previous publication [ 32 ]. (D) IGV Genome browser views of m 6 A-marked retrotransposon subfamilies during mouse OET. (E) Box plots showing the expression levels of retrotransposon subfamilies with only m 6 A, only m 6 Am, or both m 6 A and m 6 Am-marked during human OET. The transcriptome datasets of human oocytes and early embryos are derived from a previous publication [ 31 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. (F) Box plots showing the expression levels of retrotransposon subfamilies with only m 6 A, only m 6 Am, or both m 6 A and m 6 Am-marked during mouse OET. The transcriptome datasets of mouse oocytes and early embryos are derived from a previous publication [ 32 ]. Boxes represent the 25th–75th percentile (line at the median), with whiskers at 1.5 × IQR. P -values in panels (E) and (F) were determined by the Kolmogorov–Smirnov test.
Similar to our analysis of gene transcripts, we also assessed the expression changes associated with the gain or loss of m6A/m6Am-marked retrotransposon RNAs across consecutive stages of human and mouse oocyte maturation and early embryonic development. The results showed that >90.6% of m 6 A- or m 6 Am-marked retrotransposon RNAs in human oocytes and early embryos remained unchanged between consecutive stages despite dynamic gains or losses of m 6 A/m 6 Am marks ( Supplementary Fig. S9G ). Likewise, in mouse oocytes and early embryos, >75.1% of m 6 A- or m 6 Am-marked retrotransposon RNAs showed no significant expression changes across consecutive stages, despite changes in methylation status ( Supplementary Fig. S9H ). These findings indicate that the gain or loss of m 6 A/m 6 Am marks on retrotransposon RNAs is largely independent of changes in transcript abundance, suggesting that these methylation dynamics are primarily governed by epitranscriptomic regulation.
Furthermore, we investigated the relevance of m 6 Am and m 6 A modifications to retrotransposon RNA expression. Overall, the RNA expression levels of the main families of LTR and LINE RNAs marked only by m 6 A were higher than their unmarked counterparts before ZGA, while this trend reversed around ZGA. In contrast, although m 6 Am-marked LTR and LINE RNAs showed no significant difference in RNA expression levels compared to unmethylated transcripts in human early development, the trend of these m 6 Am-marked RNAs was similar to that of m 6 A-marked retrotransposons (Fig. 5E ). For mouse retrotransposon RNAs, the RNA expression levels of LTR and LINE RNAs marked only by m 6 A tended to be higher than those of the unmarked counterparts throughout the entire developmental process from oocytes to early embryos. Similarly, m 6 Am-marked LTR and LINE RNAs generally exhibited higher expression levels than their unmarked counterparts during early development, except for m 6 Am-marked retrotransposons around the ZGA stage (Fig. 5F ). Thus, these findings suggest that m 6 A and m 6 Am methylation on human LTR and LINE RNAs probably shifts from stabilizing transcripts before ZGA to promoting their decay around ZGA. In contrast, in mice, LTR and LINE RNAs marked by m 6 A consistently exhibit higher expression levels than their unmarked counterparts throughout the OET. Meanwhile, m 6 Am displays a functional switch: prior to ZGA, it is positively correlated with the expression of its marked LTR and LINE RNAs, but around ZGA, this correlation becomes negative (Fig. 5F ).
Materials|Methods
This study was approved by the Medical Study Ethics Committee of The Third Affiliated Hospital of Guangzhou Medical University (2024-066) and The Affiliated Guangdong Second Provincial General Hospital of Jinan University (2022-SZ-KY-009–01), China, following the measures of the People’s Republic of China on the administration of Human Assisted Reproductive Technology, the ethical principles of the Human Assisted Reproductive Technology and the Human Sperm Bank as well as the Helsinki declaration.
All human gametes and embryos were clinically discarded and collected from volunteers between 24 and 35 years old after signing informed consent. Women diagnosed with chromosomal abnormalities, polycystic ovary syndrome, or endometriosis were excluded. For oocyte collection, controlled ovarian stimulation was carried out, and transvaginal ultrasound-guided oocyte retrieval was scheduled for 36 h after human chorionic gonadotrophin (hCG) administration. Cumulus-oocyte complexes were exposed to hyaluronidase and mechanically denuded at 38 h after hCG administration and cultured in G-IVF PLUS medium (Vitrolife). Mature MII were clinically used for intracytoplasmic sperm injection (ICSI) treatment at 40 h after hCG administration. The spare immature oocytes were clinically discarded following ICSI treatments and donated with informed consent of these infertile couples with male factors. These in vitro matured oocytes were fertilized by ICSI with donated sperm, which were used just for research. Human one-cell, four-cell, eight-cell, and blastocyst embryos were vitrified and cryopreserved at the appropriate time after in vitro fertilization. Clinically discarded embryos with morphologically high quality from tripronuclear (3PN) zygotes were also collected. Following insemination, embryos were cultured in G1-PLUS medium (Vitrolife) to obtain cleavage embryos in a humidified atmosphere with 6% CO 2 and 5% O 2 . On Day 3, embryos were transferred and cultured in G2-PLUS medium (Vitrolife) to obtain blastocysts. The thawed oocytes or embryos of high quality were selected randomly for the experimental groups. Zona pellucida was removed from the oocytes or embryos by acidic Tyrode’s solution (Sigma–Aldrich, T1788). The number of embryos at each developmental stage for MeLACE-seq was shown in Supplementary Table S1 .
All animal maintenance and experimental procedures used in this study were approved by the Affiliated Guangdong Second Provincial General Hospital of Jinan University, Guangzhou, China. Mettl3 conditional knockout (cKO) mice in oocytes ( Mettl3 flox/flox ; Zp3 -cre, referred to as Mettl3 cKO) were referenced from the previous report [ 24 ]. The Mettl3 flox/flox female mice were used as the control group (referred to as control). The control and Mettl3 cKO oocytes were collected as described before. Briefly, 4- to 6-week-old females of both control and Mettl3 cKO mice were injected with 5 International Unit (IU) of pregnant mare’s serum gonadotropin (PMSG), and germinal vesicle (GV) oocytes were collected from the ovaries of control and Mettl3 cKO female mice at 48 h post PMSG injection.
To collect GV oocytes from 8-week-old wild-type mice, the whole ovaries of C57BL/6 female mice were clipped mechanically with a razor blade. The zona pellucida was gently removed from oocytes by treatment with Tyrode’s solution. To obtain pre-implantation embryos, 6- to 8-week-old C57BL/6 female mice were super-ovulated by injection with 5 IU of PMSG, followed by human chorionic gonadotropin (hCG, 5 IU) 48 h later. Pre-implantation embryos were collected from the super-ovulated C57BL/6 female mice mated with DBA/2 male mice. Mouse zygote, early two-cell, late two-cell, eight-cell, morula, and blastocyst stage embryos were collected at 18, 35, 48, 65, 80, and 96 h post hCG administration, respectively. Embryos were cultured in KSOM (Millipore, MR-107-D) medium. The zona pellucida was gently removed by treatment with Tyrode’s solution (Sigma, T1788). The number of oocytes and embryos used for MeLACE-seq is shown in Supplementary Table S1 .
HEK293 cells were grown in Dulbecco’s modified Eagle’s medium (Thermo Scientific, 11995-065) with 10% newborn bovine serum plus 100 U penicillin/streptomycin (Thermo Scientific, 15240062) at 37°C in 5% CO 2 . Lipofectamine RNAiMAX (Thermo Scientific, 13778030) was used for small interfering RNA (siRNA) transfection, according to the manufacturer’s instructions.
The SELECT assay was performed as previously described [ 25 ]. Briefly, the zona pellucida was gently removed using Tyrode’s solution (Sigma, T1788). Oocytes were washed three times with 0.1% bovine serum albumin (BSA) in phosphate buffered saline (PBS) and subsequently lysed in 17 μl of lysis buffer containing 1 mM dNTPs, 1 μM forward primer, 1 μM reverse primer, 10 × CutSmart buffer, and nuclease-free H₂O. Lysis was carried out using the following temperature program: 90°C for 1 min; 80°C for 1 min; 70°C for 1 min; 60°C for 1 min; 50°C for 1 min; 40°C for 6 min. After lysis, a reaction mixture containing Bst 2.0 DNA polymerase (New England BioLabs, M0537S), SplintR ligase (New England BioLabs, M0375S), and ATP (New England BioLabs, P0756S) was added. The reaction was incubated at 40°C for 20 min, followed by enzyme inactivation at 80°C for 20 min. The resulting cDNA was quantified by real-time polymerase chain reaction (PCR) using a LightCycler 480 system (Roche) with Power SYBR Green Master Mix (Vazyme, Q111). All Reverse Transcription-quantitative PCR primers used in this study are listed in Supplementary Table S2 .
MeLACE-seq was adapted from the LACE-seq [ 23 ]. Samples were lysed in 50 µl ice-cold lysis buffer (1 × PBS, 0.1% SDS, 0.5% NP-40, 0.5% sodium deoxycholate) for 10 min on ice. The lysate was treated with 1 μl of RNase inhibitor (Thermo Scientific, EO0381) and 4 μl of RQ1 DNase (Promega, M6101) and incubated at 37°C for 3 min with gentle mixing. Then, 0.5–1 μg of anti-m 6 A antibody (Synaptic Systems, 202003) was added to the lysate and rotated for 2 h at 4°C. The solution was cross-linked twice with 0.4 J cm −2 UV light (254 nm) in a Hybrilinker.
For each sample, 10 μl protein A/G magnetic beads (Thermo Scientific, 26162) were blocked with block buffer (1 × PBS, 0.2 mg/ml glycogen, 0.2 mg/ml BSA) at room temperature (RT) for 2 h, then washed once with 0.1 M Na-phosphate buffer (93.2 mM Na 2 HPO 4 , 6.8 mM NaH 2 PO 4 , 0.05% Tween 20, pH 8.0). The blocked beads were resuspended in 50 μl of 0.1 M Na-phosphate buffer and transferred to the cross-linked solution. The mixture was incubated with rotation for 2 h at 4°C. The immunoprecipitated RNAs were fragmented on beads with 1 × 10 −8 U micrococcal nuclease (MNase, New England BioLabs, M0247S) for 3 min at 37°C. The fragmented RNA 3′ ends were dephosphorylated using FastAP alkaline phosphatase (Thermo Scientific, EF0651) at 37°C for 10 min. The 3′ linker was ligated to fragmented RNA with T4 RNA ligase 2 (New England 879 Biolabs, M0242) at RT for 2.5 h. After converting the RNA fragment to cDNA using a biotin-modified primer, the cDNA was released from Protein A/G beads by treating it with RNase H (Thermo Scientific, EN0202) and then captured using streptavidin C1 beads (Thermo Scientific, 650002). A cDNA 3′ linker was ligated to the 3′ end of cDNA with T4 RNA ligase 1 (New England BioLabs, M0437) overnight at RT, then the cDNA products were subjected to pre-PCR using KAPA HiFi HotStart Ready Mix (KAPA Biosystems, KK2601). The PCR products were purified using Ampure XP beads (Beckman Coulter, A63881) and subjected to in vitro transcription with T7 RNA Polymerase (New England BioLabs, M0251) at 37°C for 24 h. DNA template was removed by treating it with TURBO DNase (Thermo Scientific, AM2238) at 37°C for 30 min, and the RNA was purified using Agencourt RNA Clean beads (Beckman Coulter, A63987). Complementary ribosomal RNA (crRNA) fragments were depleted by hybridizing the crRNA probes to themselves and treating them with RNase H (Thermo Scientific, EN0202) at 37°C for 30 min. crRNA probes were digested with TURBO DNase (Thermo Scientific, AM2238) at 37°C for 30 min and the RNA was purified using Agencourt RNA Clean beads (Beckman Coulter, A63987). After performing reverse transcription and indexed PCR, the PCR products were subjected to size selection using a 2% agarose gel. Regions ranging from 250 to 500 bp were purified using a Gel Extraction Kit (Qiagen, 28604). Barcoded libraries were pooled and sequenced on the Illumina platforms with 150 bp paired-end reads.
The adapter sequences at both ends of the raw reads were removed using the Cutadapt [ 26 ] program (v1.15) with the following parameters: -m 18 -j 8 –max-n 4 –trim-n –times 2 -e 0.1 -O 3. Paired-end reads were merged into single reads using fastp [ 27 ] (version 0.21.0) if there was an overlap of >30 nt. After extracting the unique molecular identifier (UMI) sequence, the clean reads were initially aligned to the mouse pre-ribosomal RNA using Bowtie2 [ 28 ] software (version 2.5.1). The remaining unmapped reads were subsequently aligned to the mouse (mm9) reference genome using STAR [ 29 ] (version 2.5.2b) with the following parameters: –outSJfilterReads Unique –alignEndsType Extend5pOfRead1 –outFilterMismatchNoverLmax 0.04 –outFilterMismatchNmax 999. UMI sequences were used to remove PCR duplicates, and the retained reads were utilized for peak identification. For motif analysis, LACE-seq peaks were first extended by 30 nt upstream, and overrepresented motifs in the extended sequences were identified using the findMotifsGenome.pl function in Homer ( http://homer.ucsd.edu/homer/ ). The Pearson correlation coefficient between MeLACE-seq replicates was calculated using deepTools multiBigwigSummary (‘bins’ mode and window size = 10 kb) and plotCorrelation.
To define an m 6 Am peak, we followed previously published procedures. Specifically, A putative peak was identified as a position with a stack of at least five reads, and the peak was selected if it contained an adenosine (±1) and was located within 50 bp of the TSS [ 30 ]. The gene was defined as an m 6 Am + gene if any of its gene transcripts or isoforms overlapped with at least one m 6 Am peak; otherwise, the gene was considered an m 6 Am- gene.
m 6 A peak and gene were defined as follows: the 5′ ends of reads were stacked, and positions with fewer than 10 reads were filtered out. The remaining positions within 100 bp were merged into clusters. For each cluster, the highest site was selected iteratively, with a minimum distance of 50 bp between selected sites. Sites within 50 bp of an m 6 Am peak were excluded from being defined as m 6 A peaks. A gene was defined as an m 6 A + gene if its gene transcripts or isoforms overlapped with at least one m 6 A peak.
We reanalyzed the RNA-seq and translatome data from human and mouse oocytes and early embryos, following the instructions of the previous report [ 31 , 32 ]. For retrotransposon RNA analysis, multi-aligned reads were retained and subjected to TEtranscripts to quantify retrotransposon RNA families and conduct differential expression analysis. For published translatome data, the translation efficiency (TE) of different developmental stages for each gene was calculated by the ratio of Ribo-lite and mRNA-seq [FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) + 1/FPKM + 1].
To analyze the conservation and specificity of m6A and m6Am methylome between human and mouse, the human and mouse homologous genes were downloaded from Mouse Genome Database (MGD; http://www.informatics.jax.org/homology.shtml ) [ 33 ]. We chose the homolog gene pairs with only one homolog gene in each species, and those with multiple homolog gene pairs in either species were discarded.
The four major categories of m 6 A or m 6 Am marked genes, OET gain, oocyte-specific, OET loss, and embryonic, were identified as follows the cutoffs based on m 6 A and m 6 Am peaks:
Oocyte-specific: m 6 A or m 6 Am was deposited on RNAs only at the GV stage or from the GV stage and occurred in at least two consecutive stages, but not including the blastocyst.
OET gain: RNAs were unmarked by m 6 A or m 6 Am at the GV stage, but gained m 6 A or m 6 Am at any stage after the GV stage, except for blastocyst.
OET loss: m 6 A or m 6 Am was deposited on RNAs at the GV stage and lost m 6 A or m 6 Am marks at any stage after the GV stage. This category did not include oocyte-specific genes.
Embryonic: m 6 A or m 6 Am was deposited on RNAs only at the blastocyst stage or on RNAs from ZGA (L2C for mouse and 8C for human) and till the blastocyst stage.
Gene Ontology enrichment analyses were performed using DAVID [ 34 ].
As described in a previous study, we defined maternally decayed genes, including M-decay, Z-decay, and C-decay genes.
Maternally decayed genes: the RNA FPKM ≥ 5 at the GV or MII stage and the expression fold change (mouse GV/L2C or human GV/8C) ≥ 3.
M-decay genes: the expression fold change (GV/1C) ≥ 2 and the expression fold change (mouse GV/L2C or human GV/8C) < 2.
Z-decay genes: the expression fold change (GV/1C) ≥ 1 and the expression fold change (GV/1C) < 2 and the expression fold change (mouse 1C/L2C or human 1C/8C) ≥ 2.
C-decay genes: the expression fold change (GV/1C) ≥ 2 and the expression fold change (mouse 1C/L2C or human 1C/8C) ≥ 2.
Maternally decayed genes were the combination of M-decay, Z-decay, and C-decay RNAs.
Minor and major ZGA genes were identified as follows:
Mouse minor ZGA genes: the RNA FPKM < 5 at the GV and MII stage, the expression fold change (E2C/1C) ≥ 5;
Human minor ZGA genes: the RNA FPKM < 5 at the GV and MII stage, the expression fold change (4C/1C) ≥ 5.
Mouse major ZGA genes: the RNA FPKM < 5 at the GV and MII stage, the expression fold change (L2C/1C) ≥ 5;
Human major ZGA genes: the RNA FPKM < 5 at the GV and MII stage, the expression fold change (8C/1C) ≥ 5.
ZGA genes were a combination of minor ZGA and major ZGA genes.
We collected human and mouse transcription factors from the database AnimalTFDB4 [ 35 ]. Only the transcription factors with the FPKM ≥ 1 at any developmental stage were used to analyze their m 6 A and m 6 Am marking status.
All experiments were independently repeated at least twice, and no inconsistent results were observed. Statistical analyses were carried out using the GraphPad software or R Studio. The boxplot borders represent upper and lower quartiles (25th and 75th percentiles), and the center line represents the median. The statistical tests and P -values are indicated in the figure legends. P -values <0.05 were considered significant. All data are reproducible, and the details of replicates are stated in the figure legends.