RNA N6-methyladenosine (m6A, m6A) modification is a critical regulator for a range of physiological processes. However, the dynamic m6A profiles within human preimplantation embryos remain uncharacterized. Here, we present the first RNA m6A landscape of single human oocytes and early embryos. Comparative analyses with mouse data reveal an intriguing divergence during the window of zygotic genome activation. m6A-modified genes are involved in regulation of gene transcription, while unmodified genes are mainly associated with basic metabolic processes. Maternal decay mRNAs exhibit a propensity for m6A modifications, and these genes are targeted by miRNAs. m6A modified genes that are constantly expressed across all stages demonstrate higher translation efficiency. Moreover, we observe frequent m6A enrichment on stage-specifically expressed retrotransposons, particularly within young subfamilies. m6A inhibitor leads to m6A erasure on massive retrotransposons. In summary, this study provides a resource to broaden our understanding about the regulatory roles of m6A during early human embryo development.
Synopsis
This study presents the landscape of RNA m6A modification in human oocytes and early embryos, uncovering dynamic changes during zygotic genome activation and retrotransposon regulation. Functional analyses reveal key roles for m6A in regulating RNA stability and translation.
•
Dynamic m6A atlas of human oocytes and preimplantation embryos.
•
m6A marks and miRNAs co-regulate mRNA degradation, and m6A enhances the translation efficiency of stably expressed genes.
•
Dynamic retrotransposons, especially young subfamilies, are frequently m6A-modified; m6A inhibition erases these m6A marks.
•
m6A-modified genes are involved in transcriptional regulation, while unmodified genes drive basic metabolic processes.
Introduction
Human infertility is a growing global concern, with one in six individuals worldwide affected by this condition (Cox et al, 2022). Understanding the processes of human oocyte maturation and early embryonic development is crucial to deciphering its etiology. RNA modifications have emerged as pivotal regulators of gene expression programs (Shi et al, 2019; Yang et al, 2018). In mammals, the most prevalent internal mRNA modification is N6-methyladenosine (m6A, m6A) (Frye et al, 2018; Roundtree et al, 2017). m6A modification plays crucial roles in various biological processes related to reproduction and fertility (Klungland et al, 2017; Zheng et al, 2013). The key role of m6A dynamics in meiosis seems to be well conserved and was identified early in pioneering studies in yeast (Schwartz et al, 2013; Shah and Clancy, 1992). We and others have characterized the dynamic m6A modification in mouse oocytes and early embryos, uncovering its important functions in early embryonic development (Wang et al, 2023; Wu et al, 2022; Yao et al, 2023; Zhu et al, 2023). However, no study has yet investigated this crucial RNA modification in human early embryos. This is primarily due to two major challenges: first, human embryos are precious and difficult to obtain as research materials, as their use is strictly regulated by ethical constraints; second, various methods have been developed to identify transcriptome-wide m6A profiles, significantly advancing our understanding of m6A (Chen et al, 2015a; Dominissini et al, 2012; Garcia-Campos et al, 2019; Hu et al, 2022; Linder et al, 2015; Liu et al, 2023; Meyer, 2019a; Meyer et al, 2012; Pratanwanich et al, 2021; Shu et al, 2020; Wang et al, 2020; Xiao et al, 2023; Zhang et al, 2019), however, given the scarcity of human early embryos, none of these methods can be directly applied to this context, as they would require thousands of embryos. Recently, two unique droplet-based methods, scDART-seq and m6A-CT, have been developed to enable m6A profiling of single cells, however, both require thousands of cells as starting materials (Hamashima et al, 2023; Tegowski et al, 2022). Another method, scm6A-seq, which barcodes individual embryos, enabled the first m6A profiling of early mouse embryos by using multiple embryos as the starting material for each developmental stage (Yao et al, 2023). Recently, we developed picoMeRIP-seq to study m6A using single oocyte or embryo as starting material (Li et al, 2024). picoMeRIP-seq opens the door for m6A mapping in vivo in single cells and scarce cell types. Earlier studies have shown that reprogramming of DNA methylation is highly similar between mouse and human early embryos (Guo et al, 2014; Smith et al, 2014). Understanding whether the m6A landscapes in early embryonic RNAs are conserved between humans and mice is essential for comprehensive insights on RNA m6A methylation. Here, we provide a valuable resource on the RNA m6A landscape during the human oocyte-to-embryo transition. To the best of our knowledge, this is the first study to explore transcriptome-wide RNA modifications in human oocytes and early embryos, deepening our understanding of m6A’s regulatory roles in early human embryo development.
Results
We recently reported a highly sensitive picogram-scale m6A RNA immunoprecipitation and sequencing (picoMeRIP-seq) method that can start out with a single mouse oocyte, or a single mouse or zebrafish early embryo (Li et al, 2024). Here, we applied picoMeRIP-seq to human oocytes and early embryos to profile the m6A landscape and to identify potential regulators during human oocyte-to-embryo transition (OET) (Fig. EV1A; “Methods”).
Figure EV1. Performance evaluation of picoMeRIP-seq in human samples using hESCs and m6A profiling in human oocytes and preimplantation embryos; along with the RNA expression levels of m6A-related regulatory factors.
(A) Illustration of picoMeRIP-seq procedure performed on human oocytes and early embryos. The figure was created in BioRender. (B) Genome browser snapshots of two genes with m6A enrichment in hESCs. (C) Transcriptome-wide correlation analyses of picoMeRIP-seq experiments (IP samples) among samples from 1000, 100 and 10 cells in hESCs. (D) Heatmap showing the fraction of m6A-modified genes (identified in samples of 10 cells) over those identified in samples of 1000 and 100 cells, as well as over those from a previous study in hESCs (Batista et al, 2014). (E) Metagene profiles showing the enrichment of m6A peaks along protein-coding genes in hESCs. (F) Consensus motifs identified within m6A peaks in hESCs. (G) Principal component analysis (PCA) of transcriptome-wide m6A signals. Multiple biological replicates for each stage were generated. (H) Number of m6A peaks. For each stage, the averaged (mean) number across biological replicates 1 and 2 are shown. (I) Fraction of m6A+ genes under different expression cutoffs. (J) Boxplot showing mRNA levels of m6A associated proteins in human oocytes and early embryos. The data were downloaded from a previously published database EmAtlas (Zheng et al, 2023). The boxplots are arranged from top to bottom as follows: the maximal value, the 75th percentile, the median, the 25th percentile, and the minimal value. TB trophoblast.
First, we demonstrated that picoMeRIP-seq could acquire high-quality data from as few as 10 human embryonic stem cells (hESCs) (Fig. EV1B–F). As shown in genome browser, a high level of reproducibility for selected gene loci was observed (Fig. EV1B). Transcriptome-wide profiles of picoMeRIP-seq data from 1000 cells to 10 cells showed strong correlations between replicates and across different input cell numbers (Pearson correlation coefficients >0.8; Fig. EV1C). A high percentage of genes with m6A identified in samples of 10 cells showed an overlap with those identified in samples of 1000 cells (>85%), and with those reported in previously published m6A-seq data (>75%) (Batista et al, 2014) (Fig. EV1D). In agreement with previous studies (Batista et al, 2014; Dominissini et al, 2012; Meyer et al, 2012), m6A peaks were predominantly enriched around stop codons (Fig. EV1E), and the most significant consensus motif within m6A peaks was GGACU (Fig. EV1F). These results highlight the efficacy of picoMeRIP-seq in accurately profiling the m6A methylome from low-input human samples, demonstrating its strong sensitivity and reproducibility.
Next, we performed picoMeRIP-seq on single human oocytes at the germinal vesicle (GV) and metaphase II (MII) stages, and single human embryos at the zygote (1C), two-cell (2C), eight-cell (8C) and blastocyst (BLT) stages, with ≥2 biological replicates for each stage (Table EV1). Principal component analysis (Fig. EV1G) and hierarchical clustering analysis (Fig. 1A) of the transcriptome-wide m6A profiles showed the precise clustering of single oocytes and embryos by cell identity. For the 2C and 8C stages, we had only two biological replicates available. To maintain statistical rigor and ensure comparability between stages, we selected two biological replicates from each stage for downstream analyses. The number of m6A peaks identified ranged from 5401 to 9601 across the various stages (Fig. EV1H). We identified 5126 (GV), 4497 (MII), 3666 (1C), 5328 (2C), 4820 (8C) and 5519 (BLT) m6A-modified (m6A+) genes (Fig. 1B; Dataset EV1). Among highly expressed genes (TPM ≥ 10), the percentage of m6A+ genes ranged from 32 to 44%, aligning with the overall trend observed in the number of m6A+ genes across the six stages (Fig. EV1I). Overall, the numbers and the percentage of m6A+ genes decreased until fertilization and then increased after the zygote stage, the minimal number of m6A-modified genes at the 1C stage could result from significant maternal RNA degradation before substantial zygotic mRNA synthesis begins, which suggests that m6A is dynamically regulated during the human preimplantation embryos. This could be credited to the dynamic expression of the m6A methyltransferases METTL3 and METTL14, the m6A demethylases FTO and ALKBH5, and m6A readers at both RNA and protein levels (Yang et al, 2018; Zheng et al, 2023; Zhu et al, 2025) (Figs. EV1J and EV2A). Among these m6A+ genes, the majority were protein-coding genes (82–92%), but there were also long non-coding RNAs (lncRNAs, 5–12%), pseudogenes (2–6%) and a small fraction of other types of RNA species (<1%) (Fig. EV2B,C). The non-coding m6A + RNAs identified at each stage in this study may play potential regulatory roles during human oocyte and early embryonic development, as demonstrated by previous study that established the m6A-associated functions of lncRNAs in human fetal tissues (Xiao et al, 2019). For instance, m6A has been demonstrated to be required for the function of XIST (Patil et al, 2016), which is a lncRNA that mediates gene silencing on the X chromosome during mammalian development and was marked by m6A at the BLT stage (Fig. 1C). The m6A profiles were clearly detectable as peaks, as illustrated by genome browser snapshots of representative genes, such as oocyte-secreted factor GDF9 (Paulini and Melo, 2011), the 8C embryo marker gene LEUTX (Lewin et al, 2023), blastocyst-specifically activated gene ARID3A (Rhee et al, 2014), and an imprinted gene ZDBF2 that was constantly expressed across all stages (Duffie et al, 2014) (Fig. 1C). The m6A peaks typically exhibited significant enrichment near the stop codons (Figs. 1D,E and EV2D), and displayed the consensus motif GGACU (Figs. 1F and EV2E) across all stages.
Figure 1. Dynamic RNA m6A landscape in human oocytes and preimplantation embryos.
(A) Hierarchical clustering analysis of transcriptome-wide m6A profiles. (B) Number of m6A+ genes. (C) Representative genome browser view of picoMeRIP-seq signals. (D) Bubble plot representing the relative enrichment ratios in different genomic features. For each feature, the enrichment score was calculated as the log2 ratio of the observed peak numbers to the expected peak numbers. (E) Metagene profiles of m6A peak distribution. (F) Consensus motifs on m6A peaks identified in biological replicate 1.
Figure EV2. The distribution of m6A peaks across different developmental stages; along with the protein expression levels of m6A-related regulatory factors.
(A) Bar charts showing protein levels of m6A associated proteins in human oocytes and early embryos. The data were downloaded from a previously published database (Zhu et al, 2025). The mean values were shown as the bar. (B) Number of m6A+ genes for different types of genes. (C) Fraction of m6A+ genes among different types of RNAs. (D) Top, pie charts show genomic annotation of m6A peaks. Bottom, bar plots show the peak enrichment score, which was calculated as the log2 ratio of the observed over expected peak numbers. (E) Consensus motifs identified within m6A peaks from biological replicate 2.
We further explored the extent to which differences in the m6A methylome during OET are conserved between human and mouse by comparing their transcriptomes and m6A methylomes from equivalent stages (Wang et al, 2023) (Fig. 2A). To mitigate the influence of RNA abundance on m6A characterization in the two species, we specifically targeted genes with high expression levels (TPM ≥ 50) at each stage in either species, noting that m6A+ genes had median TPM around 50 (Fig. EV3A). For the non-homologous genes between the two species, an average of 39% were m6A+ across all stages in human, compared to 49% in mouse, with a notably low fraction of m6A modification at the BLT stage for both human (20%) and mouse (35%) (Fig. EV3B).
Figure 2. Conservation and divergence of m6A methylation during human and mouse OET.
(A) Schematic of human and mouse oocytes and embryos used for stage comparison. Human ZGA mainly occurs at the 8C stage, while mouse ZGA occurs at the 2C stage. (B) Heatmaps showing the expression levels (left) of three categories of homologous genes in humans and mice for each developmental stage; and the ratios of genes with different m6A statuses (right). (C, D) Comparison of gene lengths (C) and RRACH counts (D) between m6A modified and unmodified homologous genes at the human 8C and mouse 2C stages. The P values were calculated by the two-sided Wilcoxon rank-sum test. The boxplots are generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot. N number of genes used in the data analysis. (E) GO analysis for m6A modified and unmodified human-specifically expressed genes. Fisher’s exact test was used to calculate the one-sided P values.
Figure EV3. Comparison of gene length and RRACH motif between m6A modified and unmodified genes in human and mouse.
(A) Statistics of expression level of m6A+ genes in human and mouse data. TPM, transcripts per million. (B) The ratios of non-homologous genes (TPM ≥ 50) with m6A modification in human and mouse. (C, D) Comparison of gene lengths (C) and RRACH counts (D) between m6A modified and unmodified genes for human- (left) and mouse- (right) specifically expressed genes, as well as human-mouse co-expressed genes (middle). The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of genes used in the data analysis. (E) The average RRACH motif density after normalized per 1 kb length between m6A modified and unmodified genes for human- (left) and mouse- (right) specifically expressed genes, as well as human-mouse co-expressed genes (middle). The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of genes used in the data analysis. The boxplots in (C–E) were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot.
Next, we focused on homologous genes and categorized them into three groups based on their expression at each stage: human-mouse co-expressed (expressed in both human and mouse), human-specifically expressed, and mouse-specifically expressed genes (Fig. 2B; Dataset EV2; “Methods”). Among human-mouse co-expressed genes, approximately 35% were m6A+ in both species across various stages, with <11% uniquely modified in human while >17% uniquely modified in mouse (Fig. 2B, middle). Notably, 36% of co-expressed genes at the 8C/2C stage were exclusively modified at the mouse 2C stage versus only 7% at the human 8C stage. Human major zygotic genome activation (ZGA) occurs at the eight-cell stage, whereas in mice it takes place at the two-cell stage. This observed divergence could be attributed to differences in developmental timing or the distinct profiles of major zygotic gene expression, comparisons of human and mouse transcriptomes during the major ZGA show that only half of the transcriptomes overlap, suggesting that the homology between human and mouse zygotic genes is relatively low (Sha et al, 2020). This disparity is further reflected in the genes that are specifically expressed in either human or mouse. Regarding human-specifically expressed genes, over 50% were m6A+ at the GV, MII and BLT stages, whereas this decreased to 35% at the 8C stage (Fig. 2B, top). In contrast, in mouse-specifically expressed genes, >54% were m6A+ at every stage, with a peak of 65% at the 2C stage (Fig. 2B, bottom). These data reveal distinct RNA m6A methylome states during the human major ZGA, likely reflecting evolutionary divergence.
Recent advances in understanding m6A deposition mechanisms indicate that m6A acts like a default hardware feature (He et al, 2023; Uzonyi et al, 2023; Yang et al, 2022). These studies found that longer genes tend to exhibit higher m6A levels, with the exon junction complex (EJC) acting as an m6A suppressor to shape the m6A epitranscriptome, protecting RNA near exon junctions within coding sequences from methylation. Also, longer genes typically have more exons, which increases the availability of RRACH (R = G or A; H = A, C or U) motifs in exonic regions. We aim to determine if this finding represents a universal mechanism or context-dependent differences in human early embryonic development. For both human- and mouse-specifically expressed genes, m6A+ genes exhibited significantly longer gene lengths than unmodified genes across all stages (Figs. 2C and EV3C). In human-mouse co-expressed genes, m6A+ genes in both species were substantially longer than others; and interestingly, human homologous m6A+ genes were significantly shorter than mouse m6A+ homologous genes specifically at the 2C and 8C stage, respectively (Fig. 2C), a distinction not observed at other stages (Fig. EV3C). Again, this could be attributed to the unequal cell stage comparison between the mouse 2C and human 8C stages. m6A+ genes contained markedly more RRACH motifs than unmodified genes across all stages; and for human genes co-expressed in mouse, the genes were uniquely expressed in m6A + mRNAs in human had lower counts of RRACH motifs than those uniquely expressed in m6A + mRNAs in mouse, specifically in human 8C embryo, likely due to the unequal cell stage comparison (Figs. 2D and EV3D). It is important to note that, after normalizing for gene length, there is no significant difference in the density of RRACH motifs (Fig. EV3E). Instead, this may be more closely related to gene length and splicing events, as supported by recent studies (He et al, 2023; Uzonyi et al, 2023; Yang et al, 2022). These analyses indicate that our data quality is consistent with the latest findings in the field, demonstrating that picoMeRIP is a highly accurate method for RNA m6A profiling.
Gene ontology (GO) analysis provided functional insights into different groups of genes: of human-specifically expressed genes at the GV, MII, 1C and 8C stages, m6A+ genes were largely linked to transcription regulation while unmodified genes were mainly enriched in basic metabolic processes; at the BLT stage, human-specifically expressed m6A+ genes contributed to cell differentiation (Fig. 2E). Additionally, for mouse-specifically expressed m6A+ genes, divergence in function by stage was observed: they were implicated in the regulation of cell proliferation and multicellular organism development from GV to the 1C stage, and in transcription regulation from GV to the 2C stage (Fig. EV4A). Across all stages, for human-mouse co-expressed genes, m6A+ genes were largely associated with transcription, embryonic development as well as development-related signal pathways, while unmodified genes were frequently involved in diverse metabolic pathways (Fig. EV4B,C). Taken together, these results underscore the conserved characteristics of m6A modification between the two species. The observed divergences between human 8C and mouse 2C embryos suggest the possibility of species-specific roles for m6A during the critical window of major ZGA.
Figure EV4. GO analysis of m6A modified and unmodified genes for mouse-specifically expressed genes and human-mouse co-expressed genes.
(A–C) GO analysis for m6A + /− genes for mouse specifically expressed genes (A), for human genes that were co-expressed in mouse (B), and for mouse genes that were co-expressed in human (C). Fisher’s exact test was used to calculate the one-sided P values.
We sought to further focus on m6A profiling of the two critical events during OET: maternal RNA degradation (MD) and ZGA (Fig. 2A). We classified MD genes into two subsets based on their distinct decay patterns: M-decay, degraded during oocyte maturation, exemplified by VCAN that is associated with human oocyte developmental competence (Shen et al, 2020); and Z-decay, degrading after fertilization, such as BNC1 that is required for human oocyte meiosis (Zhang et al, 2018) (Fig. 3A,B; Dataset EV3; “Methods”). In total, we identified 1339 MD genes (M-decay, 396; Z-decay, 943) and 550 ZGA genes (Fig. 3A). The smaller number of m6A-marked genes detected in this study compared to mouse studies is likely influenced by the larger number of embryos (at least hundreds of embryos at each stage) used in mouse studies, which allows for greater enrichment of low-abundance RNAs (Wang et al, 2023; Wu et al, 2022; Zhu et al, 2023). Near half of ZGA genes were methylated at the 8C stage (Fig. 3C), including TFAP2C, a known marker of human ground-state naive pluripotency (Pastor et al, 2018) (Fig. 3B). Notably, some MD genes were consistently tagged with m6A throughout development. We hypothesize that although the original maternal mRNAs are degraded, the zygotic genome re-expressed new mRNAs from the same genes, ensuring continued support for embryonic development (Sha et al, 2020; Yan et al, 2013), consistent with the results in mouse (Wang et al, 2023; Wu et al, 2022). In contrast to M-decay genes (<19%), Z-decay genes (>30%) displayed significantly higher m6A+ proportions at the GV, MII and 1C stages (Fig. 3C), which can be partially explained by the relatively lower expression of M-decay genes (Fig. EV5A). For M-decay genes at the GV stage, Z-decay genes at the GV, MII and 1C stages, as well as ZGA genes at the 8C stage, those with m6A modification had significantly longer gene lengths (Figs. 3D and EV5B) and more RRACH motifs (Figs. 3E and EV5C) than unmodified genes. Unmodified MD genes were predominantly associated with metabolic pathways, apoptotic process and protein transport, whereas unmodified ZGA genes were linked to RNA processing and nucleosome assembly (Fig. 3F). Of note, m6A+ genes from both MD and ZGA functioned in regulation of cell proliferation and pluripotency, likely mediated by key factors essential for transcription (Fig. 3F). Given that m6A is believed to regulate cell fate in embryonic stem cells by modulating its abundance on transcription factor coding genes (Batista et al, 2014; Geula et al, 2015; Jin et al, 2021), we examined the expression and m6A modification dynamics of genes encoding transcription factors (Fig. 3G) and transcription cofactors (Fig. EV5D) throughout all six stages. Among the highly expressed genes (TPM ≥ 50), an average of 68% of transcription factors and 48% of transcription cofactors (Fig. EV5E) were m6A+ across all stages, with the highest methylation ratios at the BLT stage. The mature human blastocyst is comprised of the first three cell lineages of the embryo: trophectoderm, epiblast and primitive endoderm (Alberio, 2020; Blakeley et al, 2015; Stirparo et al, 2018). We further examined the m6A status on the lineage-specific marker genes across developmental stages and observed m6A occupancy on some master transcription regulators at the 8C and/or BLT stages, such as GATA2, GATA3 and TEAD3 (markers of trophectoderm), SOX2 and NANOG (markers of epiblast), and PDGFRA and GATA6 (markers of primitive endoderm) (Fig. EV5F). The prevalence of m6A deposition on maternal RNAs and the onset of new m6A methylation during ZGA, particularly on those genes involved in transcription regulation, hint at the involvement of m6A in these two critical events during human OET.
Figure 3. Enrichment of m6A in MD and ZGA genes, as well as in transcription factor coding genes.
(A) Expression (left) and m6A modification (right) dynamics of M-decay, Z-decay and ZGA genes. (B) Genome browser views showing the m6A profiles of representative M-decay, Z-decay and ZGA genes. (C) Percentage of m6A+ genes for each event at each specified stage. The P values are calculated by the one-sided Fisher’s exact test. (D, E) Comparison of gene lengths (D) and RRACH counts (E) between m6A modified and unmodified genes for M-decay, Z-decay and ZGA genes. The P values were calculated by the two-sided Wilcoxon rank-sum test. The boxplots were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot. N number of genes used in the data analysis. (F) GO analysis for m6A modified and unmodified genes for MD and ZGA genes. Fisher’s exact test was used to calculate the one-sided P values. (G) Expression (left) and m6A status (right) dynamics of transcription factor coding genes.
Figure EV5. Comparison of m6A-modified and unmodified genes in terms of M-decay, Z-decay, and ZGA categories; along with the expression levels and m6A modification status of transcription factors and cofactors.
(A–C) Comparison of expression levels (A), gene lengths (B), and RRACH counts (C) between m6A modified and unmodified genes for M-decay, Z-decay and ZGA genes. The P values were calculated by the two-sided Wilcoxon rank-sum test. The boxplots were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot. N number of genes used in the data analysis. (D) Expression (left) and m6A status (right) dynamics of transcription cofactor coding genes. (E) Proportion of m6A+ transcription factors (left) and transcription cofactors (right) coding genes at each stage with different expression thresholds. (F) Expression and m6A marking status of genes essential for the lineage specification events in human early embryos. Top: gene expression, down: m6A status.
MicroRNAs (miRNAs), along with m6A, are another key regulator of RNA stability and translation, and are essential for mammalian development (DeVeale et al, 2021). We and others have observed that m6A modifications are particularly enriched at miRNA targeting sites in various human and mouse tissues (Liu et al, 2020) as well as in mouse pluripotent stem cells and early embryos (Chen et al, 2015b; Wang et al, 2023). Next, we examined the co-occupancy between m6A modification and miRNAs on human MD and ZGA genes (Dataset EV4; “Methods”). miRNAs had a significant tendency to target m6A + Z-decay genes as compared with unmodified genes at the GV, MII and 1C stages (Figs. 4A and EV6A,B). However, unlike our observations in mouse data (Chen et al, 2015b; Wang et al, 2023), we did not discern a distinct preference of miRNA targeting between modified and unmodified M-decay genes at the GV stage and ZGA genes at the 8C stage (Figs. 4A and EV6A,B). Previous studies have demonstrated significant differences in gene expression and miRNA profiles between humans and mice during oocyte maturation and major ZGA. We reasoned that mice and human may have divergent mechanisms for miRNA regulation, which likely explains the distinct miRNA targeting phenomenon observed here (Paloviita et al, 2021; Sha et al, 2020; Yang et al, 2016; Zhao et al, 2020). This discrepancy may suggest species-specific differences in the coordination of m6A modification and miRNA targeting during OET between humans and mice. However, when examining genes that were constantly expressed across all stages, the preference of miRNA targeting on m6A+ genes was substantially higher than on unmodified genes (Figs. 4B and EV6C,D).
Figure 4. Association of m6A modification with miRNA targeting and translation efficiency.
(A) Comparison of the proportion of mRNA genes targeted by miRNAs between m6A modified and unmodified genes for M-decay, Z-decay and ZGA genes. The P values are calculated by the one-sided Fisher’s exact test. (B) Comparison of the count of miRNAs targeting m6A modified versus unmodified genes for constantly expressed mRNA genes. The P values are calculated by the two-sided Wilcoxon rank-sum test. N number of genes used in the data analysis. (C, D) Comparison of translation efficiency between m6A modified and unmodified genes for M-decay, Z-decay and ZGA mRNA genes (C), as well as for constantly expressed mRNA genes (D). The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of genes used in the data analysis. (E) Comparison of translation efficiency in constantly expressed mRNA genes with m6A across different genic regions. N number of genes used in the data analysis. The boxplots in (B–E) were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot.
Figure EV6. Comparison of miRNA targeting and translation efficiency between m6A modified and unmodified genes for M-decay, Z-decay and ZGA genes, as well as for constantly expressed genes; along with the analysis of m6A enrichment on retrotransposon-derived RNAs.
(A) Comparison of the count of miRNAs targeting m6A modified versus unmodified mRNA genes for M-decay, Z-decay and ZGA genes. The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of genes used in the data analysis. (B) Comparison of the normalized count of miRNAs targeting m6A modified versus unmodified mRNA genes for M-decay, Z-decay and ZGA genes. See “Methods” for the calculation of normalized mRNA counts. The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of genes used in the data analysis. (C) Comparison of the proportion of mRNA genes targeted by miRNAs between m6A modified and unmodified genes for constantly expressed genes. The P values were calculated by the one-sided Fisher’s exact test. (D) Comparison of the normalized count of miRNAs targeting m6A modified versus unmodified mRNA genes for constantly expressed genes. See “Methods” for the calculation of normalized mRNA counts. The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of genes used in the data analysis. (E) Distribution of m6A signals along the full-length sequences of the representative retrotransposon subfamilies, based on unique assignment strategy. The mean (represented by the central black line) and standard deviation (represented by the band around the black line) across all loci were plotted. These plots are based on uniquely aligned reads. (F) Distribution of m6A signals along the full-length sequences of the representative retrotransposon subfamilies, based on a random assignment strategy. The mean (represented by the central black line) and standard deviation (represented by the band around the black line) across all loci were plotted. These plots are based on a random assignment strategy for those multiply aligned reads. (G) Heatmaps showing the number of expressed loci (left), the number of those with m6A among the expressed loci (middle), and the percentage of loci with m6A (right) for the representative retrotransposon subfamilies. These results are based on a random assignment strategy for those multiply aligned reads. (H) Comparison of the number of processed sequencing reads between control (DMSO) and treatment (STM2457) groups in human ESCs. (I) Comparison of the number of m6A peaks between control (DMSO) and treatment (STM2457) groups in human ESCs. The boxplots in (A, B, D) were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot.
m6A is known to regulate RNA translation efficiency (Jain et al, 2023; Shan et al, 2023; Wang et al, 2015). By correlating previous translatome data in human oocytes and early embryos (Zou et al, 2022) with our m6A methylome data, we found that m6A + Z-decay genes had relatively higher translation efficiency than unmodified genes at the MII and 1C stages, while there was no significant association between m6A modification and mRNA translation efficiency for M-decay and ZGA genes (Fig. 4C; Dataset EV5; “Methods”). For constantly expressed genes, the median translation efficiencies of m6A+ genes were higher than unmodified genes across different stages (Fig. 4D). The locations of m6A deposition along mRNAs are closely related to its distinct roles in translation control (Meyer, 2019b; Shan et al, 2023). We categorized constantly expressed m6A+ genes based on their modification sites into 5’UTR, coding sequence (CDS), stop codon and 3’UTR. Despite previous reports of a negative correlation between CDS m6A and translation in other research models (Meyer, 2019b; Shan et al, 2023), our observations showed that genes with CDS m6A residing had relatively higher translation efficiencies in comparison to those with m6A at other locations, especially in oocytes (Fig. 4E). The role of CDS m6A in either promoting or inhibiting protein production remains elusive (Meyer, 2019b). The variations in m6A associated translation efficiency between MD and ZGA events, as well as the divergent impacts of CDS m6A on translation between our observation and previous studies, suggest that the effect of m6A on translation is multifaceted and complex. It possibly involves different m6A readers to modulate translation positively or negatively under different conditions (Meyer, 2019b; Shan et al, 2023).
Recent findings showed that the m6A modified transcripts derived from retrotransposons could act both in cis and, in trans to regulate the chromatin environment, hence, the m6A modification has been linked to cell fate determination through modulating retrotransposon-derived RNAs (Chelmicki et al, 2021; Chen et al, 2021; Liu et al, 2021; Sun et al, 2023; Xu et al, 2021). Here, we found that m6A was largely enriched on retrotransposon-derived RNAs during human OET (Fig. 5A; Dataset EV6; “Methods”). The most recently integrated human endogenous retrovirus H (HERVH) were consistently activated across all stages, and >51% of the expressed HERVH loci/copies were modified by m6A from the GV to 2C stages (Fig. 5A). HERVK showed marked m6A modification at the BLT stage. In addition to full-length endogenous retroviruses (ERVs), solo long terminal repeats (LTRs), such as LTR12C and THE1C, were also heavily occupied by m6A prior to ZGA (Fig. 5A). Certain evolutionarily young subfamilies of LINE-1, including L1HS, L1PA2, and L1PA3, demonstrated their highest m6A frequencies at the BLT stage. SVA elements also exhibited a high degree of m6A modification concurrent with ZGA. There was no frequent enrichment of m6A on Alu elements. Notably, HERVH displayed a broad distribution of m6A across the internal sequences, while solo THE1C exhibited a clear positional preference for m6A occupancy (Figs. 5B and EV6E). Further examination of other types of retrotransposons supports this divergent distribution of m6A along their entire sequences (Fig. EV6E). Due to their highly repetitive nature, one of the technical challenges in studying retrotransposons is the difficulty in assigning multi-mapped sequencing reads to unique genomic copies/loci (Lanciano and Cristofari, 2020). Our random assessment strategy for multi-mapped sequencing reads yielded consistent conclusions, ruling out computational biases (Fig. EV6F,G). These data reveal that m6A marks retrotransposon-derived RNAs and their stage-specific expression during the OET process, suggesting a potential role in the post-transcriptional regulation of both maternal and ZGA retrotransposon-derived RNAs. It is likely that the widespread distribution of m6A marks on retrotransposons is largely ‘hard-wired’ by the genomic architecture (He et al, 2023; Uzonyi et al, 2023; Yang et al, 2022).
Figure 5. m6A deposition on human retrotransposon-derived RNAs.
(A) Heatmaps showing the number of expressed loci (left), the number of those with m6A among the expressed loci (middle), and the percentage of loci with m6A (right) for the representative retrotransposon subfamilies. These results are based on uniquely aligned reads. (B) Distribution of m6A signals along the full-length HERVH and solo THE1C sequences (based on uniquely aligned reads from biological replicate 1). The mean (represented by the central black line) and standard deviation (represented by the band around the black line) across all loci were plotted. (C) Comparison of the number of retrotransposon-associated m6A peaks between DMSO control (Ctrl) and STM2457 treatment groups in human ESCs. (D) Comparison of the m6A signal on representative retrotransposon subfamilies between DMSO and STM2457 groups in hESCs (based on uniquely aligned reads from biological replicate 1). The P values were calculated by the two-sided Wilcoxon rank-sum test. The boxplots were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot. N number of repeat loci used in the data analysis.
Human reproduction is regulated by retrotransposons (Macfarlan et al, 2012; Xiang et al, 2022) and studies on the potential intervention and influence of m6A status on retrotransposons in human embryo development is intriguing. Due to the limited availability of human early embryos and ethical considerations, investigating m6A status and retrotransposon activity in human early embryos is currently not feasible. Instead, to observe the possible effect of inhibiting m6A levels of retrotransposon-derived RNAs, we resorted to use hESCs as a model. We treated hESCs with STM2457, a highly potent and selective catalytic inhibitor of the m6A writer METTL3 (Yankova et al, 2021). As expected, at the transcriptome level, STM2457 treatment led to a considerable reduction of uniquely mapped sequencing reads (>50%) (Fig. EV6H), and a substantial loss of ~73% in the number of m6A peaks (Fig. EV6I). In the context of retrotransposons, this treatment resulted in an around 47% reduction in m6A peaks associated with retrotransposons (Fig. 5C). The observed variation in reduction rates could be due to differences in transcript structure and the inherent transcriptional regulation of retrotransposons. Retrotransposons, as repetitive elements, may have sequence-specific features that make them less sensitive to m6A modification machinery inhibitors. Of note, m6A levels on retrotransposons, including those from LTR (e.g., HERVH, HERVK, HERVL, HERV9NC, MLR1B, MSTA, THE1C) and LINE1 (e.g., L1HS, L1PA2, L1PA3), were significantly diminished after STM2457 treatment (Figs. 5D and EV7A–C). Following STM2457 treatment, more than half of the genes involved in transcription and metabolic pathways lost m6A modification, yet their expression levels remained largely unchanged (Fig. EV7D–F). The retrotransposon stage-specific expression patterns during human OET, combining with the analysis of METTL3 inhibitor experiment collectively highlight the critical role of m6A in labeling active retrotransposon transcripts.
Figure EV7. Analysis of m6A profile changes after STM2457 treatment in hESC.
(A) Comparison of the m6A signal on representative retrotransposon subfamilies between DMSO and STM2457 groups in hESCs (based on uniquely aligned reads from biological replicate 2). The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of repeat loci used in the data analysis. (B, C) Comparison of the m6A signal on representative retrotransposon subfamilies between DMSO and STM2457 groups in hESCs (based on a random assignment strategy for those multiply aligned reads; see “Methods”). The (B) is biological replicate 1; and (C) is biological replicate 2. The P values were calculated by the two-sided Wilcoxon rank-sum test. N number of repeat loci used in the data analysis. (D–F) The expression level and fold change (left), m6A status (middle) and percentage (right) of transcription factor (D), transcription cofactors (E) and metabolic pathways (F) mRNAs. The boxplots in (A–C) were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile; and the mean value indicated by a dot.
To the best of our knowledge, no studies to date have explored RNA modifications transcriptome-wide in human oocytes and embryos, leaving the roles of RNA modifications in human preimplantation embryos unknown. Here, we first explored this long-standing question of how conserved RNA m6A reprogramming is between mouse and human early embryos. We observed the highly dynamic landscape of the m6A methylome on both coding and non-coding RNAs (Fig. 6). Recent advances in understanding the mechanism of m6A deposition suggest that m6A functions akin to a hardware feature that is added by default (He et al, 2023; Uzonyi et al, 2023; Yang et al, 2022). Its specificity is globally regulated by suppressors that prevent m6A deposition in unmethylated transcriptome regions. EJC have been identified as an m6A suppressor that protect exon junction–proximal RNA within coding sequences from being methylated. Longer genes tend to have more exons and splicing events, increasing the availability of RRACH motifs in exonic regions. Our findings are consistent with this unified model; indeed, even during the transition from oocyte to embryo, m6A-modified mRNAs consistently exhibit longer gene lengths and a greater number of RRACH motifs. We also revealed distinct m6A methylome patterns during human major ZGA compared to mice. These differences likely reflect evolutionary adaptations to species-specific needs. We explored the potential functions of m6A in the degradation of M-decay and Z-decay maternally expressed genes (Fig. 6). m6A-modified transcripts were significantly targeted by miRNAs, pointing to a possible collaborative mechanism in RNA decay. We also observed that m6A-modified genes, which were consistently expressed across all developmental stages, exhibited higher translational efficiency (Fig. 6). These observations prompted several important questions for future research, such as the precise mechanisms by which m6A promotes RNA decay and translation. We found m6A modifications are frequently enriched in stage-specifically expressed retrotransposons, especially within young subfamilies such as HERVH, L1HS, and SVA_D (Fig. 6). It will be interesting to determine whether co-transcriptional incorporation of m6A in retrotransposons is related with ZGA. We also evidenced the retrotransposon m6A levels was downregulated by METTL3 inhibitor STM2457 in hESCs (Fig. 6). Due to the limited availability of human embryos, we were unable to conduct large-scale single-embryo analyses to explore heterogeneity or experimentally validate these intriguing questions. By utilizing picoMeRIP-seq to explore the RNA m6A landscape in human embryos, we aim to inspire further research into RNA modifications and their roles in early human embryonic development.
Figure 6. Summary of the RNA m6A methylome landscape during mammalian OET.
In summary, our study offers a comprehensive overview of the dynamic m6A methylome in human oocytes and embryos, highlighting both conserved and human-specific features compared to mice (Fig. 6). The mechanisms by which m6A influences mRNA stability and translation require further elucidation to fully understand the epitranscriptomic regulatory programs governing early mammalian embryonic development. Future research should aim to unravel these intricate processes to provide deeper insights into the role of m6A in developmental biology and its implications for human reproduction health.
To comply with ethical permissions and the European GDPR regulations, human embryo work approval was obtained from regional committees for medical and healthcare research ethics (Regionale komiteer for medisinsk og helsefaglig forskningsetikk, reference ID 2016/1263). hESCs work adheres to the International Society for Stem Cell Research (ISSCR) guidelines. Human oocytes and sperm for the study were sourced from patients undergoing fertility treatment at the Department of Reproductive Medicine, Oslo University Hospital (OUH). These materials were collected only after obtaining voluntary informed consent from the patients and their partners. Patients were given the option of continued storage, disposal, or donation of embryos to this research (RNA modifications during early human embryo development). Patients received no compensation and were offered counseling. They could withdraw their participation at any time until the embryo was used for research. This project did not affect the patients’ IVF cycles or normal treatment processes. After embryo evaluation and selection for patient treatment by individual embryologists, only suboptimal quality oocytes and embryos, deemed unsuitable for clinical treatment, were collected for this project. Embryos were not cultured beyond day 5 post-fertilization. Suboptimal oocytes and embryos were collected and transported immediately to the laboratory for processing.
Human ESC, oocyte and early embryo collection
Patients underwent controlled ovarian stimulation with recombinant or urine-derived follicle stimulating hormone (FSH). Spontaneous ovulation was prevented by either co-administration of GnRH antagonist or suppression with GnRH agonist from the mid-luteal phase of the preceding cycle. Final follicular maturation was induced by administering human chorionic gonadotropin (hCG) upon observation of ≥3 follicles with a diameter of ≥17 mm via ultrasound monitoring. Oocyte retrieval was performed 36 h after hCG administration under transvaginal ultrasound guidance. Retrieved embryos were cultured in Vitrolife G1, G2, or G-TL media under controlled conditions of 6% CO2 and 37 °C, in accordance with standard clinical protocols at the Department of Reproductive Medicine, OUH. Healthy sperm samples were cryopreserved in liquid nitrogen using cryoprotectants. Immature oocytes, identified by the presence of a germinal vesicle (GV), were selected for analysis. Single matured (MII) oocytes from fertile, healthy donors were included for single-cell analysis. These oocytes underwent vitrification and warming using RapidVit and RapidWarm oocyte kits (Vitrolife, Sweden). Intra cytoplasmic sperm injection (ICSI) was employed for fertilization. Fertilization was assessed post-ICSI, and embryo cleavage was monitored every 12 h. Suboptimal grade embryos, deemed unsuitable for cryopreservation, were donated to the project with patient consent. Single human oocytes and embryos were manually sorted into 12 μl of 1× lysis buffer (Takara). The samples were snap-frozen in liquid nitrogen and stored at −80 °C until further analysis.
Human embryonic stem cells (H1-hESC cell line) obtained from ATCC were cultured in feeder-free conditions. Regular mycoplasma testing was conducted, confirming the absence of mycoplasma contamination in all cell lines. For METTL3 inhibition, STM2457 (MCE, HY-134836) was added into the medium at a final concentration of 5 μM for 2 days (Yankova et al, 2021). Cells were sorted into 1× lysis buffer (Takara) in an 8-strip PCR tube using a BD FACSMelody cell sorter (BD Biosciences) following the manufacturer’s instructions. Plates were sealed and stored at −80 °C until further processing.
picoMeRIP-seq experiments
The whole procedures were conducted within a UV decontaminated LAF bench.
rRNA and DNA depletion
To achieve single-tube rRNA and DNA depletion, we employed the NEBNext rRNA depletion kit (NEB) with slight modifications to the manufacturer’s protocol. Initially, RNA/probe master mix (3 µl) was added to a 12-µl sample. Subsequently, the sample underwent a precise temperature ramp from 95 °C to 22 °C, controlled at a rate of −0.1 °C s–1, followed by a 5-min incubation period at 22 °C. The RNase H reaction mix (5 µl) was then added to the samples, which were incubated at a precisely controlled temperature of 37 °C for 30 min. Following this, the samples were treated with 30 µl of DNase I digestion mix and incubated at 37 °C for an additional 30 min. The resulting samples were purified using 2.2× volume of RNAClean XP beads, washed twice with 80% freshly prepared ethanol, and finally eluted with 78 μl of nuclease-free water. To prevent RNA degradation, 2 μl RiboLock RNase inhibitor (40 U µl–1) was added to each sample, ensuring a final sample volume of 80 μl.
RNA fragmentation by sonication
Sonication of the samples was carried out using a UP100H Ultrasonic Processor (Hielscher) equipped with a 2-mm probe. The samples underwent 2 × 30 s sonication cycles, precisely alternating between 30 s of sonication and 30 s on ice for each cycle. After sonication, 20% RNA fractions from several oocytes or embryos were pooled as input controls for each developmental stage. Add appropriate nuclease-free water to ensure a final volume of 80 μl. After that, 20 μl 5× IP buffer was added to the samples to achieve a final volume of 100 µl for subsequent processing steps.
Antibody–bead incubation
Prior to use, Dynabeads (Invitrogen) underwent washing to ensure purity. Specifically, 20 μl of beads was washed twice with 1× IP buffer, with each washing step involving vortexing, rapid centrifugation on a MiniGalaxy, and placement on a magnetic rack to remove the supernatant. Following this, the antibody solution was prepared by diluting 4 µl of Millipore anti-m6A (Cat num: ABE572; Lot num: 3194595) antibody in 16 µl of 5× IP buffer (50 mM Tris-HCl (pH 7.5), 750 mM NaCl, 0.5% (vol/vol) NP-40, 5 U/µl RiboLock) and 60 µl of nuclease-free water. The antibody-containing solution (80 µl) was added to the washed beads, and the antibody–bead mixture was incubated overnight with head-over-tail rotation on a HulaMixer at 4 °C (40 r.p.m.).
Immunoprecipitation (IP) and washes
The prepared antibody-coated beads were captured on the tube wall using a magnetic rack, and the supernatant from the antibody–bead incubation was carefully discarded. The antibody-coated beads then underwent washing to remove unbound antibodies that could otherwise compete for binding to the epitope. The beads were washed twice with 200 μl of 1× IP buffer, with each wash involving vortexing (four times for 5 s each) to ensure thorough washing. After vortexing during the second wash, the antibody-coated beads were transferred to 0.2-ml PCR tubes, and a volume of 10 μl of homogenously antibody-coated bead solution was transferred to each PCR tube. The tubes were rapidly centrifuged on a MiniGalaxy and placed on a magnetic rack for at least 2 min or until the solution became clear. Following removal of the supernatant, 100 μl of the sonicated RNA sample was added to each antibody–bead-containing tube, and the samples were incubated with head-over-tail rotation on a HulaMixer at maintained 4 °C for 2 h (40 r.p.m.). Subsequently, the tubes were swiftly centrifuged on a MiniGalaxy and placed on a magnetic rack. The supernatant was meticulously removed, and the RNA–antibody–bead complexes underwent a series of washes to eliminate non-specific binding. The RNA–antibody–bead complexes were washed four times in the following solutions, quickly spun and placed in a magnetic rack in between washes: washed once with ice-cold medium-stringency RIPA buffer (10 mM Tris-HCl (pH 8.0), 300 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 1% (vol/vol) Triton X-100, 0.2% (vol/vol) SDS and 0.1% (vol/vol) sodium deoxycholate), washed twice with ice-cold high-stringency RIPA buffer (10 mM Tris-HCl (pH 8.0), 350 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 1% (vol/vol) Triton X-100, 0.23% (vol/vol) SDS and 0.1% (vol/vol) sodium deoxycholate) and washed once with ice-cold medium-stringency RIPA buffer. Following the four wash steps, the tubes were centrifuged and placed on a magnetic rack, and the supernatant was carefully discarded. The RNA–antibody–bead complexes were then resuspended in 100 μl of 1× IP buffer and incubated for 5 min. Following this, the samples were swiftly centrifuged and placed on a magnetic rack, and the supernatant was removed. The RNA–antibody–bead complexes were resuspended in 147.9 μl of elution buffer (5 mM Tris-HCl (pH 7.5), 1 mM EDTA, 0.05% (vol/vol) SDS, and 1 U µl–1 RiboLock RNase inhibitor). Proteinase K (2.1 μl; NEB) was added to each tube, and the tubes were then incubated on a Thermomixer at 1200 r.p.m. and 55 °C for 1.5 h. Following incubation, the tubes were briefly centrifuged and incubated further on a Thermomixer at 80 °C for 20 min to inactivate the Proteinase K. Subsequently, the samples were placed on a magnetic rack for 2–3 min, and the supernatant containing the m6A-immunoprecipitated RNA was carefully transferred to a new 1.5-ml low-binding tube. The remaining beads were resuspended again in 147.9 μl of elution buffer, and 2.1 μl of Proteinase K was added. The samples were immediately placed on a Thermomixer at 1200 r.p.m. and 55 °C for 5 min, followed by inactivation of Proteinase K on a Thermomixer at 80 °C for 20 min. The tubes were then placed back on a magnetic rack for 2–3 min, and the supernatant was collected and pooled with the first supernatant in the same 1.5-ml low-binding tube to recover as much of the m6A-immunoprecipitated RNA as possible, resulting in a total volume of ~300 μl.
Ethanol precipitation
For both input and immunoprecipitated RNA samples, nuclease-free water was added to each tube to achieve a final volume of 400 μl. Subsequently, 40 μl of 3 M sodium acetate (pH 5.2; Thermo Fisher Scientific) and 10 μl of linear acrylamide 5 mg μl–1 (Thermo Fisher Scientific) were added to each sample, followed by 1000 μl of ice-cold 100% ethanol. The samples were vigorously vortexed without centrifugation or spinning and immediately placed at –80 °C for at least 2 h or overnight until completely frozen. Once the samples were recovered from –80 °C, they were allowed to briefly thaw on ice, ensuring visual confirmation that all samples had thawed before commencing centrifugation. The samples were then centrifuged at 20,000× g at 4 °C for 15 min, and the supernatant was removed without disturbing the visible pellet. The pellet was carefully washed twice with 1 ml of ice-cold 75% ethanol. For each wash, 75% ethanol was added, centrifugation was then repeated with 10 min at 4 °C. Following the last wash, as much as possible of the supernatant was removed, and the tube lid was left open until all ethanol had evaporated, resulting in a dried pellet. The dried pellet was resuspended in 7 μl of nuclease-free water for library preparation.
Library preparation and sequencing
With modifications to the manufacturer’s protocol, as described below, the SMART-Seq stranded kit (Takara, 634442) was used to construct sequencing libraries. For the fragmented input or immunoprecipitated RNA, the protocol was performed without the fragmentation step. After the first PCR amplification and following AMPure bead purification, the beads were resuspended by adding 46.5 µl of nuclease-free water, and the ribosomal cDNA depletion protocol was skipped. The samples were incubated at room temperature for 5 min for rehydration, and 46 µl of supernatant was recovered from each sample. Subsequently, the protocol was followed until completion. The libraries were assessed for quantity using KAPA library quantification kits (Roche), and size distribution was evaluated using TapeStation D1000 ScreenTape (Agilent Technologies). These assessments provided a good estimation of pooling at equimolar ratios. The pooled libraries were sequenced on a NovaSeq system (Illumina).
picoMeRIP-seq data processing and m6A peak calling
Quality control, alignment, and processing of sequencing reads
The quality of the raw sequencing reads for both Input and IP samples was evaluated using FastQC (v0.11.8) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and sequencing adapters were trimmed using Cutadapt (v1.8.1) (Marcel, 2011) with settings “-m 20 --max-n 0.01 --trim-n”. The trimmed reads were aligned to the human reference genome (hg38) using HISAT2 (v2.1.0) (Kim et al, 2015) with the parameter “−5 8 --no-mixed --no-discordant”. We retained only those read pairs aligning uniquely to one genomic location, as indicated by HISAT2. PCR duplicates were removed using SAMtools (v1.9) (Li et al, 2009) fixmate & markdup, and the read pairs aligning to ribosomal RNAs (per GENCODE v39 annotation) were filtered out using BEDTools (v2.28.0) (Quinlan and Hall, 2010) intersect function.
Gene and transcript quantification from the Input samples
Based on GENCODE (v39) gene annotation, the expression abundance (TPM, transcript per million) of genes and transcript/isoforms were estimated by StringTie (v1.3.5) (Pertea et al, 2015) with the parameter “-e -A”. Further, the normalization of expression abundances among samples was performed by the R package NormExpression (v0.1.0) (Wu et al, 2019) with the normalization factor (method = ”DESeq”).
Definition of m6A signal
The genome coverage bigWig files with a resolution of 10 bp (i.e., bin size = 10 bp) and normalization according to RPKM (reads per kilobase per million reads) were generated by deepTools (v3.2.0) (Ramirez et al, 2014) bamCoverage with the parameter “-bs 10 --normalizeUsing RPKM”. The visualization of read density along exonic regions of genes was based on these bigWig files. For each gene, the exonic regions of the transcript/isoform with highest expression value was used to visualize the read density. The transcripts/isoforms used for genome browser visualization are ENST00000687138.1 (GDF9), ENST00000638280.2 (LEUTX), ENST00000263620.8 (ARID3A), ENST00000429829.6 (XIST), ENST00000649650.1 (ZDBF2), ENST00000513960.5 (VCAN), ENST00000345382.7 (BNC1), and ENST00000201031.3 (TFAP2C).
For each bin (size = 10 bp), the m6A signal was calculated as the log2 transformed ratio of (the IP sample’s RPKM value + 1) over (the Input sample’s RPKM value + 1).
PCA, Pearson correlation, and hierarchical clustering analyses
By deepTools multiBigwigSummary (“bins” mode and window size = 1 kb) & plotPCA & plotCorrelation, we performed principal component analysis (PCA), and calculated Pearson correlation coefficients (r) along with hierarchical clustering of the samples. In the hierarchical clustering, the dissimilarity distance between any two samples was defined as 1.0 − r.
m6A peak calling
MACS2 (v2.1.2) (Zhang et al, 2008) callpeak with the parameter “-g 421727543 --keep-dup all -B --nomodel --call-summits” was used identify m6A peaks (q value < 0.05), with the IP sample and the corresponding Input sample.
Genomic annotation, metagene profiling, and motif analysis of m6A peaks
Peak annotation was carried out using BEDTools based on the GENCODE (v39) annotation library. To avoid assigning a peak to multiple genomic features, the hierarchy of genomic features prioritized as follows: stop codon (ranging from 200 bp upstream to 200 bp downstream of the annotated stop codon in GENCODE), 3’ UTR, 5′ UTR, CDS (coding sequencing), exon, intron, and intergenic. The enrichment score at each genomic region was calculated as the log2 ratio of the observed over expected peak numbers. The expected number of peaks is calculated based on the length proportion of different genomic regions (i.e., 5’UTR, CDS, 3’UTR, stop codon, exon, intron, intergenic) relative to the total length of reference genome hg38.
Metagene profiles of m6A peak summits along protein coding genes were generated using MetaPlotR (Olarerin-George and Jaffrey, 2017). For each gene, only the concatenated exonic regions of the transcript/isoform with the highest expression value were used for plotting metagene profiles.
We searched for consensus m6A motifs using the peak regions as detected by MACS2 and the regions (within a 100 bp, 150 bp, 200 bp, 300 bp, and 400 bp windows centered on peak summits as detected by MACS2), and by Homer (v4.11.1) (Heinz et al, 2010) findMotifsGenome.pl with the parameter “-rna -len 5,6,7,8”.
Definition of m6A gene
We defined a gene as m6A+ if ≥1 of its transcript/isoform intersected with ≥1 m6A peak; otherwise, the gene was considered as m6A − . Considering that we had two IP samples for each developmental stage, we opted for a union approach: considering a gene as m6A+ if it was m6A+ in any of two IP samples.
For mouse samples (oocytes at the GV and MII stages, and embryos at the zygote/1C, 2-cell/2C and blastocyst/BLT stages), we derived their expression abundance and m6A profiles from our previous study (Wang et al, 2023).
Defining co-expressed and species-specifically expressed genes: within a given developmental stage that was comparable between human and mouse, a gene having an expression level of TPM ≥ 50 in both species was classified as a co-expressed genes; a gene with TPM ≥ 50 in human but TPM ≤ 10 in mouse was denoted as a human-specifically expressed gene; a gene with TPM ≥ 50 in mouse but TPM ≤ 10 in human was denoted as a mouse-specifically expressed gene.
Definition of M-decay, Z-decay, ZGA, and constantly expressed genes
M-decay gene: (1) TPM was >10 in GV oocyte; (2) TPM in GV oocyte was greater than in MII oocyte; (3) the expression fold change of GV oocyte versus 8-cell was >2; (4) the expression fold change of GV versus zygote was >2; and (5) the expression fold change of GV versus MII was <2.
Z-decay gene: (1) the TPM was >10 in GV oocyte; (2) the expression fold change of zygote versus GV was <1.5; (3) the expression fold change of GV versus 8-cell was >2; (4) the expression fold change of GV versus zygote was <2; and (5) the expression fold change of zygote versus 8-cell was >2.
MD genes were the combination of M-decay and Z-decay genes.
ZGA gene: (1) TPM < 1 in both GV and MII oocytes; (2) TPM > 10 at the 8-cell stage; and (3) TPM at the 8-cell is greater than at the zygote stage.
Constantly expressed genes: Based on gene expression values across all Input samples, we first calculated the entropy specificity by the R package BioQC (v1.10.0) (Zhang et al, 2017). A gene with <0.5 entropy specificity scores and with TPM ≥ 10 across all samples was defined as constantly expressed one.
Annotation of transcription factor, transcription cofactor and imprinted genes
The annotation of human transcription factors and transcription cofactors was downloaded from the database AnimalTFDB4 (Hu et al, 2019) on August 30, 2023. Human imprinted genes were obtained from Geneimprint database (https://www.geneimprint.com/site/genes-by-species) on August 30, 2023.
miRNA targeting analyses
The miRNA expression data (read counts) from four developmental stages (GV, MII, zygote/1C and 8-cell/8C) were obtained from the Dataset EV2 of a previous study (Paloviita et al, 2021). We further normalized the expression of miRNA to reads per million (RPM) within each sample. Given that there were multiple replicates for each developmental stage in the study (Paloviita et al, 2021), we computed the mean RPM value for each miRNA by averaging across all replicates of a given developmental stage.
To reduce the likelihood of false-positive miRNA target predictions, we focused on those miRNAs that were from conserved miRNA families and are highly expressed (RPM ≥ 10). Only the conserved gene targets of the conserved miRNAs were considered based on the miRNA and target gene annotations from TargetScanHuman (v8.0) (Agarwal et al, 2015). A gene can be targeted by multiple different miRNAs. A gene was considered miRNA-targeted at a particular stage if it was predicted to be a target of ≥1 miRNA; and we then counted the number of regulated miRNAs and the weighted number (i.e., the sum of log10 (RPM) of all targeted miRNAs) of regulated miRNAs.
Of note, only the protein coding genes (mRNAs) were used for miRNA analysis.
Analyses of human translatome data
We obtained translatome (Ribo-lite, single-end) and corresponding transcriptome (mRNA-seq, paired-end) data from a previous study (Zou et al, 2022), available at the NCBI GEO under accession number GSE197265, with two replicates for each developmental stage. The quality of sequencing reads was evaluated by FastQC, and the sequencing adapters were trimmed using Cutadapt with the parameter “cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -u 3 -m 18 --max-n 0.01 --trim-n <input.file> | cutadapt -a AAAAAAAAAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA -m 18 -o <output.file > -” for Riob-lite data, and with the parameter “cutadapt -a AAGCAGTGGTATCAACGCAGAGTACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA -A AAGCAGTGGTATCAACGCAGAGTACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT -G AAGCAGTGGTATCAACGCAGAGTAC -g AAGCAGTGGTATCAACGCAGAGTAC -m 20 --max-n 0.01 --trim-n” for mRNA-seq data. Using the same procedures as picoMeRIP-seq data (see above), both Ribo-lite and mRNA-seq data underwent alignment (HISAT2 with the default parameter for Ribo-lite, and with the parameter “−5 8 --no-mixed --no-discordant” for mRNA-seq), extraction of uniquely aligned read pairs, removal of PCR duplicates and rRNA-derived read pairs, and estimation of gene abundance (FPKM, Fragments Per Kilobase of transcript per Million mapped reads). We averaged the FPKM values over two replicates for each developmental stage, and focused on highly expressed mRNA genes (FPKM ≥ 10, estimated by mRNA-seq data) to calculate translation efficiency as the fold change of (Ribo-lite’s FPKM + 1) versus (mRNA-seq’s FPKM + 1).
Of note, only the protein coding genes (mRNAs) were used for translation analysis.
Retrotransposon analyses
Retrotransposon annotation
For LINE, SINE and SVA retrotransposons, we extracted the annotation using UCSC Table Browser with the setting “clade = Mammal, genome = Human, assembly = GRCh38/hg38, group = Repeats, track = RepeatMasker, table=rmsk” on March 8, 2020. For LTR retrotransposon, we downloaded the annotation from HERVd database (version 3) (Paces et al, 2004) on January 2, 2024. To minimize the effect of highly-fragmented/truncated retrotransposon copies/loci which lack transcription and transposition activities (Cordaux and Batzer, 2009), (1) for LINE, SINE, SVA and solo LTR, we imposed a completeness threshold whereby only genomic loci/copies with ≥90% completeness were included in the analyses; completeness being defined by the ratio of the annotated retrotransposon sequence length in the human reference genome to that of the transposable element reference sequence’s full length. Further, for the full-length ERV but not for solo LTR, only loci/copies exhibiting the complete architecture (comprising 5’LTR, internal sequence and 3’LTR) were considered. To avoid the effects of the expression of regular genes on retrotransposon analyses, we excluded the retrotransposon loci/copies that intersected with exonic regions of regular genes (annotated by GENCODE v39).
Expression and m6A enrichment of retrotransposons
For each retrotransposon locus/copy defined above, we took into account two metrics: expression value and m6A signal value. The expression value was computed as the maximum RPKM observed across all bins (size =10 bp) overlapping the locus, using Input sample. The m6A signal value was calculated as the maximum m6A signal value across all bins (size =10 bp) overlapping the locus (see “Definition of m6A signal” under the section “picoMeRIP-seq data processing and m6A peak calling”). For each developmental stage, the mean values (expression or m6A signal) across two biological replicates were assigned to this locus. A locus/copy was deemed “expressed” if its expression value was >0; and it was considered as “m6A modified” if its m6A signal value was >0.
Visualization of boxplots
All boxplots in this study (except Fig. EV1J) were generated with five percentile-based whiskers, arranged from top to bottom as follows: the 95th percentile, the 75th percentile, the median, the 25th percentile, and the 5th percentile. Additionally, each boxplot includes a single point/dot representing the mean value. The boxplot of Fig. EV1J was directly downloaded from the database EmAtlas, with arranged from top to bottom as follows: the maximal value, the 75th percentile, the median, the 25th percentile, and the minimal value.
GO analysis
We performed GO analysis with the online tool DAVID (v6.8) (Huang et al, 2009) using all human and mouse genes as background, respectively.
Data availability
Both raw sequencing data and processed data for hESCs, and processed data for human oocytes and embryos generated in this study are available in GEO under accession number GSE268645. The public m6A-seq data for hESCs were obtained from the study (Batista et al, 2014) with the SRA accession numbers SRR1035221, SRR1035222, SRR1035223 and SRR1035224. The expression values of miRNAs at the human GV, MII, 1C and 8C stages were obtained from the Dataset EV2 of the study (Paloviita et al, 2021). The translatome (Ribo-lite) and the corresponding transcriptome (mRNA-seq) data were obtained from the study (Zou et al, 2022) with the GEO accession number GSE197265.
Yanjiao Li: Conceptualization; Data curation; Formal analysis; Validation; Investigation; Visualization; Methodology; Writing—original draft; Writing—review and editing. Yunhao Wang: Conceptualization; Data curation; Software; Formal analysis; Validation; Investigation; Visualization; Methodology; Writing—original draft; Writing—review and editing. Aylin Cengiz: Resources; Investigation; Writing—review and editing. Kang-Xuan Jin: Resources; Investigation; Methodology. Blanca Corral Castroviejo: Resources. Xiaolin Lin: Resources. Marie Indahl: Resources. Rujuan Zuo: Resources. Trine Skuland: Writing—review and editing. Madeleine Fosslie: Writing—review and editing. Maria Biba: Resources; Supervision. Xuechen Wu: Resources. Peter Fedorcsak: Resources; Supervision. Magnar Bjørås: Resources; Supervision. Adam Filipczyk: Resources. John Arne Dahl: Methodology; Writing—review and editing. Gareth D Greggains: Resources; Supervision; Writing—review and editing. Kin Fai Au: Conceptualization; Supervision; Funding acquisition; Writing—review and editing. Arne Klungland: Conceptualization; Supervision; Funding acquisition; Project administration; Writing—review and editing.
Source data underlying figure panels in this paper may have individual authorship assigned. Where available, figure panel/source data authorship is listed in the following database record: biostudies:S-SCDT-10_1038-S44318-025-00474-5.
Disclosure and competing interests statement
The authors declare no competing interests.
Acknowledgements
This work was supported by South-Eastern Norway Regional Health Authority, Grant 2018086 (to AK); Research Council of Norway, FRIPRO Researcher Project 275286 (to AK), and partially by the Research Council of Norway through its Centres of Excellence scheme, project number 332713, CRESCO; National Institutes of Health grants R01HG011469 and R01GM136886 (to KFA and YW); Institutional fund from the Gilbert S Omenn Department of Computational Medicine and Bioinformatics, University of Michigan (to KFA and YW); South-Eastern Norway Regional Health Authority Early Career Grant 2016058 and Grants 2018063; 2023098 (to JAD); Research Council of Norway, FRIPRO Grant 289467 (to JAD); Oslo University Hospital Strategic Research Environment Grant (to JAD and GDG). The sequencing service was provided by the Norwegian Sequencing Centre (www.sequencing.uio.no), a national technology platform hosted by the University of Oslo and supported by the “Functional Genomics” and “Infrastructure” programs of the Research Council of Norway and the Southeastern Regional Health Authorities”.
Blakeley P, Fogarty NM, del Valle I, Wamaitha SE, Hu TX, Elder K, Snell P, Christie L, Robson P, and Niakan KK Defining the three cell lineages of the human blastocyst by single-cell RNA-seq Development 2015 142 3151-3165
Chelmicki T, Roger E, Teissandier A, Dura M, Bonneville L, Rucli S, Dossin F, Fouassier C, Lameiras S, and Bourc’his D m(6)A RNA methylation regulates the fate of endogenous retroviruses Nature 2021 591 312-316
Chen C, Liu W, Guo J, Liu Y, Liu X, Liu J, Dou X, Le R, Huang Y, Li C, et al. Nuclear m(6)A reader YTHDC1 regulates the scaffold function of LINE1 RNA in mouse ESCs and early embryos Protein Cell 2021 12 455-474
Chen K, Lu Z, Wang X, Fu Y, Luo GZ, Liu N, Han D, Dominissini D, Dai Q, Pan T, et al. High-resolution N(6) -methyladenosine (m(6) A) map using photo-crosslinking-assisted m(6) A sequencing Angew Chem Int Ed Engl 2015 54 1587-1590
Chen T, Hao YJ, Zhang Y, Li MM, Wang M, Han WF, Wu YS, Lv Y, Hao J, Wang LB, et al. m6A RNA methylation is regulated by microRNAs and promotes reprogramming to pluripotency Cell Stem Cell 2015 16 289-301
Cox CM, Thoma ME, Tchangalova N, Mburu G, Bornstein MJ, Johnson CL, and Kiarie J Infertility prevalence and the methods of estimation from 1990 to 2021: a systematic review and meta-analysis Hum Reprod Open 2022 2022 hoac051
Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, Cesarkas K, Jacob-Hirsch J, Amariglio N, Kupiec M, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq Nature 2012 485 201-206
Duffie R, Ajjan S, Greenberg MV, Zamudio N, Escamilla del Arenal M, Iranzo J, Okamoto I, Barbaux S, Fauque P, and Bourc’his D The Gpr1/Zdbf2 locus provides new paradigms for transient and dynamic genomic imprinting in mammals Genes Dev 2014 28 463-478
Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon-Divon M, Hershkovitz V, Peer E, Mor N, Manor YS, et al. Stem cells. m6A mRNA methylation facilitates resolution of naive pluripotency toward differentiation Science 2015 347 1002-1006
Guo H, Zhu P, Yan L, Li R, Hu B, Lian Y, Yan J, Ren X, Lin S, Li J, et al. The DNA methylation landscape of human early embryos Nature 2014 511 606-610
Hamashima K, Wong KW, Sam TW, Teo JHJ, Taneja R, Le MTN, Li QJ, Hanna JH, Li H, Loh YH (2023) Single-nucleus multiomic mapping of m(6)A methylomes and transcriptomes in native populations of cells with sn-m6A-CT. Mol Cell
Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, and Glass CK Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities Mol Cell 2010 38 576-589
Hu H, Miao YR, Jia LH, Yu QY, Zhang Q, and Guo AY AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors Nucleic Acids Res 2019 47 D33-D38
Hu L, Liu S, Peng Y, Ge R, Su R, Senevirathne C, Harada BT, Dai Q, Wei J, Zhang L, et al. m(6)A RNA modifications are measured at single-base resolution across the mammalian transcriptome Nat Biotechnol 2022 40 1210-1219
Huang DW, Sherman BT, and Lempicki RA Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources Nat Protoc 2009 4 44-57
Jain S, Koziej L, Poulis P, Kaczmarczyk I, Gaik M, Rawski M, Ranjan N, Glatt S, and Rodnina MV Modulation of translational decoding by m6A modification of mRNA Nat Commun 2023 14 4784
Jin KX, Zuo R, Anastassiadis K, Klungland A, Marr C, and Filipczyk A N6-methyladenosine (m(6)A) depletion regulates pluripotency exit by activating signaling pathways in embryonic stem cells Proc Natl Acad Sci USA 2021 118 e2105192118
Lewin TD, Blagrove JR, and Holland PWH Rapid evolution of the embryonically expressed homeobox gene LEUTX within primates Genome Biol Evol 2023 15 evad097
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, and Genome Project Data Processing S The Sequence Alignment/Map format and SAMtools Bioinformatics 2009 25 2078-2079
Li Y, Wang Y, Vera-Rodriguez M, Lindeman LC, Skuggen LE, Rasmussen EMK, Jermstad I, Khan S, Fosslie M, Skuland T, et al. Single-cell m(6)A mapping in vivo using picoMeRIP-seq Nat Biotechnol 2024 42 591-596
Linder B, Grozhik AV, Olarerin-George AO, Meydan C, Mason CE, and Jaffrey SR Single-nucleotide-resolution mapping of m6A and m6Am throughout the transcriptome Nat Methods 2015 12 767-772
Liu C, Sun H, Yi Y, Shen W, Li K, Xiao Y, Li F, Li Y, Hou Y, Lu B, et al. Absolute quantification of single-base m(6)A methylation in the mammalian transcriptome using GLORI Nat Biotechnol 2023 41 355-366
Liu J, Gao M, He J, Wu K, Lin S, Jin L, Chen Y, Liu H, Shi J, Wang X, et al. The RNA m(6)A reader YTHDC1 silences retrotransposons and guards ES cell identity Nature 2021 591 322-326
Liu J, Li K, Cai J, Zhang M, Zhang X, Xiong X, Meng H, Xu X, Huang Z, Peng J, et al. Landscape and regulation of m(6)A and m(6)Am methylome across human and mouse tissues Mol Cell 2020 77 426-440.e426
Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, and Jaffrey SR Comprehensive analysis of mRNA methylation reveals enrichment in 3’ UTRs and near stop codons Cell 2012 149 1635-1646
Olarerin-George AO and Jaffrey SR MetaPlotR: a Perl/R pipeline for plotting metagenes of nucleotide modifications and other transcriptomic sites Bioinformatics 2017 33 1563-1564
Paces J, Pavlicek A, Zika R, Kapitonov VV, Jurka J, and Paces V HERVd: the Human Endogenous RetroViruses Database: update Nucleic Acids Res 2004 32 D50
Paloviita P, Hyden-Granskog C, Yohannes DA, Paluoja P, Kere J, Tapanainen JS, Krjutskov K, Tuuri T, Vosa U, and Vuoristo S Small RNA expression and miRNA modification dynamics in human oocytes and early embryos Genome Res 2021 31 1474-1485
Pastor WA, Liu W, Chen D, Ho J, Kim R, Hunt TJ, Lukianchikov A, Liu X, Polo JM, Jacobsen SE, et al. TFAP2C regulates transcription in human naive pluripotency by opening enhancers Nat Cell Biol 2018 20 553-564
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, and Salzberg SL StringTie enables improved reconstruction of a transcriptome from RNA-seq reads Nat Biotechnol 2015 33 290-295
Pratanwanich PN, Yao F, Chen Y, Koh CWQ, Wan YK, Hendra C, Poon P, Goh YT, Yap PML, Chooi JY, et al. Identification of differential RNA modifications from nanopore direct RNA sequencing with xPore Nat Biotechnol 2021 39 1394-1402
Ramirez F, Dundar F, Diehl S, Gruning BA, and Manke T deepTools: a flexible platform for exploring deep-sequencing data Nucleic Acids Res 2014 42 W187-W191
Rhee C, Lee BK, Beck S, Anjum A, Cook KR, Popowski M, Tucker HO, and Kim J Arid3a is essential to execution of the first cell fate decision via direct embryonic and extraembryonic transcriptional regulation Genes Dev 2014 28 2219-2232
Sha QQ, Zheng W, Wu YW, Li S, Guo L, Zhang S, Lin G, Ou XH, and Fan HY Dynamics and clinical relevance of maternal mRNA clearance during the oocyte-to-embryo transition in humans Nat Commun 2020 11 4917
Shan T, Lib F, Wen MM, Chen ZG, Li SP, Wang YF, Cheng H, and Zhou Y m6A modification negatively regulates translation by switching mRNA from polysome to P-body via IGF2BP3 Mol Cell 2023 83 4494-4508.e6
Shen QZ, Chen M, Zhao X, Liu Y, Ren XL, and Zhang L Versican expression level in cumulus cells is associated with human oocyte developmental competence Syst Biol Reprod Med 2020 66 176-184
Shu X, Cao J, Cheng M, Xiang S, Gao M, Li T, Ying X, Wang F, Yue Y, Lu Z, et al. A metabolic labeling method detects m(6)A transcriptome-wide at single base resolution Nat Chem Biol 2020 16 887-895
Smith ZD, Chan MM, Humm KC, Karnik R, Mekhoubad S, Regev A, Eggan K, and Meissner A DNA methylation dynamics of the human preimplantation embryo Nature 2014 511 611-615
Stirparo GG, Boroviak T, Guo G, Nichols J, Smith A, and Bertone P Integrated analysis of single-cell embryo data yields a unified transcriptome signature for the human pre-implantation epiblast Development 2018 145 dev158501
Sun T, Xu Y, Xiang Y, Ou J, Soderblom EJ, and Diao Y Crosstalk between RNA m(6)A and DNA methylation regulates transposable element chromatin activation and cell fate in human pluripotent stem cells Nat Genet 2023 55 1324-1335
Tegowski M, Flamand MN, and Meyer KD scDART-seq reveals distinct m(6)A signatures and mRNA methylation heterogeneity in single cells Mol Cell 2022 82 868-878.e810
Uzonyi A, Dierks D, Nir R, Kwon OS, Toth U, Barbosa I, Burel C, Brandis A, Rossmanith W, Le Hir H, et al. Exclusion of m6A from splice-site proximal regions by the exon junction complex dictates m6A topologies and mRNA stability Mol Cell 2023 83 237-251.e237
Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, Weng X, Chen K, Shi H, and He C N(6)-methyladenosine modulates messenger RNA translation efficiency Cell 2015 161 1388-1399
Wang Y, Li Y, Skuland T, Zhou C, Li A, Hashim A, Jermstad I, Khan S, Dalen KT, Greggains GD, et al. The RNA m(6)A landscape of mouse oocytes and preimplantation embryos Nat Struct Mol Biol 2023 30 703-709
Wang Y, Xiao Y, Dong S, Yu Q, and Jia G Antibody-free enzyme-assisted chemical approach for detection of N(6)-methyladenosine Nat Chem Biol 2020 16 896-903
Wu Y, Xu X, Qi M, Chen C, Li M, Yan R, Kou X, Zhao Y, Liu W, Li Y, et al. N(6)-methyladenosine regulates maternal RNA maintenance in oocytes and timely RNA decay during mouse maternal-to-zygotic transition Nat Cell Biol 2022 24 917-927
Wu Z, Liu W, Jin X, Ji H, Wang H, Glusman G, Robinson M, Liu L, Ruan J, and Gao S NormExpression: an R package to normalize gene expression data using evaluated methods Front Genet 2019 10 400
Xiang X, Tao Y, DiRusso J, Hsu FM, Zhang J, Xue Z, Pontis J, Trono D, Liu W, and Clark AT Human reproduction is regulated by retrotransposons derived from ancient Hominidae-specific viral infections Nat Commun 2022 13 463
Xiao S, Cao S, Huang Q, Xia L, Deng M, Yang M, Jia G, Liu X, Shi J, Wang W, et al. The RNA N(6)-methyladenosine modification landscape of human fetal tissues Nat Cell Biol 2019 21 651-661
Xiao YL, Liu S, Ge R, Wu Y, He C, Chen M, and Tang W Transcriptome-wide profiling and quantification of N(6)-methyladenosine by enzyme-assisted adenosine deamination Nat Biotechnol 2023 41 993-1003
Xu W, Li J, He C, Wen J, Ma H, Rong B, Diao J, Wang L, Wang J, Wu F, et al. METTL3 regulates heterochromatin in mouse embryonic stem cells Nature 2021 591 317-321
Yan L, Yang M, Guo H, Yang L, Wu J, Li R, Liu P, Lian Y, Zheng X, Yan J, et al. Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells Nat Struct Mol Biol 2013 20 1131-1139
Yang Q, Lin J, Liu M, Li R, Tian B, Zhang X, Xu B, Liu M, Zhang X, Li Y, et al. Highly sensitive sequencing reveals dynamic modifications and activities of small RNAs in mouse oocytes and early embryos Sci Adv 2016 2 e1501482
Yang Y, Hsu PJ, Chen YS, and Yang YG Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism Cell Res 2018 28 616-624
Yankova E, Blackaby W, Albertella M, Rak J, De Braekeleer E, Tsagkogeorga G, Pilka ES, Aspris D, Leggate D, Hendrick AG, et al. Small-molecule inhibition of METTL3 as a strategy against myeloid leukaemia Nature 2021 593 597
Yao H, Gao CC, Zhang D, Xu J, Song G, Fan X, Liang DB, Chen YS, Li Q, Guo Y, et al. scm(6)A-seq reveals single-cell landscapes of the dynamic m(6)A during oocyte maturation and early embryonic development Nat Commun 2023 14 315
Zhang D, Liu Y, Zhang Z, Lv P, Liu Y, Li J, Wu Y, Zhang R, Huang Y, Xu G, et al. Basonuclin 1 deficiency is a cause of primary ovarian insufficiency Hum Mol Genet 2018 27 3787-3800
Zhang JD, Hatje K, Sturm G, Broger C, Ebeling M, Burtin M, Terzi F, Pomposiello SI, and Badi L Detect tissue heterogeneity in gene expression data with BioQC BMC Genomics 2017 18 277
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al. Model-based analysis of ChIP-Seq (MACS) Genome Biol 2008 9 R137
Zhang Z, Chen LQ, Zhao YL, Yang CG, Roundtree IA, Zhang Z, Ren J, Xie W, He C, and Luo GZ Single-base mapping of m(6)A by an antibody-independent method Sci Adv 2019 5 eaax0250
Zhao ZH, Meng TG, Li A, Schatten H, Wang ZB, and Sun QY RNA-Seq transcriptome reveals different molecular responses during human and mouse oocyte maturation and fertilization BMC Genomics 2020 21 475
Zheng GQ, Dahl JA, Niu YM, Fedorcsak P, Huang CM, Li CJ, Vågbo CB, Shi Y, Wang WL, Song SH, et al. ALKBH5 is a mammalian RNA demethylase that impacts RNA metabolism and mouse fertility Mol Cell 2013 49 18-29
Zheng L, Liang P, Long C, Li H, Li H, Liang Y, He X, Xi Q, Xing Y, and Zuo Y EmAtlas: a comprehensive atlas for exploring spatiotemporal activation in mammalian embryogenesis Nucleic Acids Res 2023 51 D924-D932
Zhu W, Ding Y, Meng J, Gu L, Liu W, Li L, Chen H, Wang Y, Li Z, Li C, et al. Reading and writing of mRNA m(6)A modification orchestrate maternal-to-zygotic transition in mice Genome Biol 2023 24 67
Zhu W, Meng J, Li Y, Gu L, Liu W, Li Z, Shen Y, Shen X, Wang Z, Wu Y, et al. Comparative proteomic landscapes elucidate human preimplantation development and failure Cell 2025 188 814-831.e821
Zou ZN, Zhang CX, Wang QY, Hou ZZ, Xiong ZQ, Kong F, Wang QJ, Song JZ, Liu BY, Liu BF, et al. Translatome and transcriptome co-profiling reveals a role of TPRXs in human zygotic genome activation Science 2022 378 41-44
This cover highlights the article Early coordination of cell migration and cardiac fate determination during mammalian gastrulation by Shayma Abukar, Kenzo Ivanovitch and colleagues. Using long-term live imaging of whole mouse embryos with light-sheet microscopy, this study tracks single cells to map how mesodermal cells coordinate their migration and fate decisions during gastrulation to drive early heart formation. The cover image shows the cardiac crescent (blue) and a late T/Bra-lineage-positive mesodermal subpopulation (red) in a DAPI-stained E7.5 mouse embryo (white). The image illustrates that late mesodermal cells contributing to the atria do not populate the cardiac crescent.
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Select your manager software from the list below and click Download.
