High-stoichiometry m6A sites are evolutionarily conserved

  1. Schraga Schwartz
  1. Department of Molecular Genetics, Weizmann Institute of Science, Rehovot 7610001, Israel
  1. Corresponding author: schwartz{at}weizmann.ac.il

  1. Handling editor: Eric Phizicky

 
Next Section

Abstract

N6-methyladenosine (m6A) is the most prevalent internal mRNA modification in eukaryotes, yet whether m6A sites are functionally important or represent neutral byproducts remains unclear. Previous evolutionary analyses failed to detect consistent conservation signatures at m6A sites, and report conflicting patterns of conservation across genic regions, such as the coding sequence (CDS) and untranslated regions (UTRs). To reconcile these inconsistencies and definitively determine whether m6A sites are under selection, we developed novel motif-level conservation metrics that incorporate knowledge of m6A biogenesis to distinguish m6A-specific selection from other confounding sources. We analyzed ∼500,000 candidate sites with quantitative, single-nucleotide resolution m6A measurements across a phylogeny spanning 447 mammalian species. After controlling for proximity to exon-junctions, we observed a clear, dose-dependent relationship between m6A stoichiometry and evolutionary conservation in both CDS and UTRs. Highly methylated sites (>60%) exhibited significantly increased conservation compared to lowly methylated sites—with an effect size approximately one-third of the typical CDS–UTR difference—providing definitive evidence of purifying selection and supporting a model where highly modified sites contribute functionally to gene regulation. We established a methodological framework for evolutionary analysis of RNA modifications, highlighting the necessity of quantitative measurements, comprehensive phylogenetic sampling, and careful consideration of modification biogenesis.

Keywords

Previous SectionNext Section

INTRODUCTION

N6-methyladenosine (m6A) is the most prevalent internal modification of messenger RNA in eukaryotes, deposited at varying levels at >200,000 sites across the transcriptome (Liu et al. 2022). The installation of m6A predominantly occurs at DRACH (D = A/G/U; R = A/G; H = A/C/U) motifs. From a functional standpoint, in mammalian cells, m6A has been mostly implicated in inducing mRNA decay (Wang et al. 2015; Du et al. 2016; Zaccara and Jaffrey 2020; Dierks et al. 2021; Uzonyi et al. 2023) but has also been implicated in regulating other steps of mRNA metabolism including export and translation, underscoring potentially diverse modes via which this modification can have functional consequences.

The m6A writer machinery is highly conserved across eukaryotes, from human, through yeast, to plants (Schwartz et al. 2013; Yue et al. 2019; Ensinck et al. 2023), pointing to a conserved function of this machinery. Consistently, disruption of the methylation machinery manifests in severe phenotypes in species sampled from across the eukaryotic lineage (Roignant and Soller 2017; Wang et al. 2018; Jiang et al. 2021; Xu et al. 2021; Perlegos et al. 2022), again pointing to a function of the machinery. This functional conservation could arise through distinct mechanisms: All m6A sites may contribute to function; only one or a few sites may be functional while others represent neutral “passengers”; or the methylation machinery itself may be conserved for methylation-independent roles—as demonstrated in yeast (Ensinck et al. 2023)—with no individual m6A sites being functionally relevant.

Diverse lines of evidence lend support to these competing hypotheses. The model where the role of the writer machinery is entirely methylation-independent is challenged by the observation that not only the writer machinery is highly conserved but also its target “DRACH”-like motif (Schwartz et al. 2013; Lence et al. 2016; Dagan et al. 2022; Ensinck et al. 2023; Brodersen and Arribas-Hernández 2024; Wang et al. 2024), suggesting that both the presence of the methylation machinery and its output are subject to selection. Further, studies in mammalian species point to a strong and causal association between the “load” of m6A across a gene and its stability (Dierks et al. 2021; Uzonyi et al. 2023), suggesting that individual m6A sites can impart a molecular outcome. There are also a growing number of cases where the functionality of individual m6A sites has been experimentally established (Porman et al. 2022; Xin et al. 2025), again supporting the functionality of individual methylation sites. On the other hand, there is also evidence to support the notion that m6A sites may to a considerable extent be functionally neutral. First, nearly all m6A sites are deposited in a nonstoichiometric manner, and the majority of methylated sites are deposited at low stoichiometries (Liu et al. 2022), indicating that at the single-transcript level, m6A deposition at specific sites is typically not essential. Second, surgical disruption of m6A often does not lead to a phenotype (Schwartz et al. 2013), and natural variation at m6A sites also often fails to result in substantial effects on gene expression (Zhang et al. 2020b; Xiong et al. 2021). Third, the permissive—rather than selective—deposition of m6A sites suggests that many sites may be incidental rather than purposeful marks. Initially thought to be a highly selective mark, m6A was recently established to be introduced by default at all DRACH consensus motifs in mammalian motifs, except within a ∼100 nt region surrounding splice junctions by the exon-junction complex (Yang et al. 2022; He et al. 2023; Luo et al. 2023; Uzonyi et al. 2023). This permissive modification mechanism suggests that the barrier on the emergence or loss of new m6A sites is low, and that point mutations arising from neutral genetic drift could substantially reshape methylation profiles by turnover of the DRACH consensus motifs. Indeed, substantial changes in methylation profiles were observed between two closely related yeast species, which could be attributed in large part to such turnover rather than changes in writer activity (Shachar et al. 2024). Potential precedent for a largely neutral model also comes from literature on modifications other than m6A, such as ADAR2-mediated adenosine deamination, which are catalyzed at thousands of sites transcriptome-wide, but where the lethality of an ADAR2 knockout is rescued by modification of only a single target (Higuchi et al. 2000).

Evolutionary sequence conservation can distinguish between these competing hypotheses: If individual m6A sites are functional, they should exhibit signatures of purifying selection. To date, there have been two key attempts to measure the evolutionary constraints acting on methylation target motifs. A first study compared the sequence conservation between mouse and human of the methylated adenosine nucleotide of m6A sites, using unmethylated sites as controls (Liu and Zhang 2018). This study observed only a subtle increase in conservation differences, considered by the authors to be too low to be supportive of functional relevance. A second study used a similar approach, and found evidence for reduced conservation of methylated As in the first and second codons, but subtly increased conservation for methylated adenosines at the third position. This study also reported a pronounced conservation for methylated sites in the 3′ UTR (Zhang et al. 2020a). These two studies thus failed to find consistent evidence for m6A conservation within coding sequences and reached conflicting conclusions about 3′-UTR conservation. However, the conclusions of these studies were undermined by severe technological and analytical limitations. First, both relied on only a few thousand m6A sites, reducing statistical power. Second, m6A detection lacked single-nucleotide precision, complicating direct site-by-site comparisons. Third, m6A levels could not be measured quantitatively, preventing correlation of evolutionary constraint with modification stoichiometry. Equally important, their evolutionary analyses were constrained by narrow phylogenetic sampling (often limited to two species) and by focusing solely on the methylated adenosine, rather than evaluating the full DRACH motif. These combined shortcomings may have undermined the ability of these studies to detect true selection acting on m6A.

Whether individual m6A sites are subject to selection, and whether such selection occurs predominantly in certain regions within genes (e.g., UTRs or CDSs), can also provide critical clues regarding the function of this mark. Earlier literature on the mark associated m6A levels with a strong, destabilizing effect independently of its relative localization within genes (Wang et al. 2014; Du et al. 2016; Dierks et al. 2021; Uzonyi et al. 2023), but more recent studies have reported that m6A primarily destabilizes transcripts when present in the CDS, but not in the 3′ UTR (Zhou et al. 2024; Linder et al. 2025; Murakami et al. 2025), raising the possibility that 3′-UTR methylation might be functionally neutral.

Here, we sought to overcome limitations of prior studies and definitively determine whether m6A sites are under selection. We develop novel motif-level conservation metrics that incorporate knowledge of m6A biogenesis, specifically the requirement of the DRACH consensus motif. This approach can distinguish m6A-specific selection from other overlapping selective constraints while also improving the statistical power to detect subtle evolutionary signals. Our analysis is further statistically empowered by use of recent advances in m6A mapping allowing comprehensive, single-nucleotide resolution quantitative mapping of this mark, and the availability of multiple genome alignments encompassing hundreds of mammal species. We identify exon–intron architecture as a critical confounder that constrained previous studies; which once controlled for, reveals a clear, dose-dependent increase in conservation as a function of methylation level. We demonstrate that individual highly modified methylated sites (>60%) within both the CDS and the UTRs are under purifying selection, whereas lowly modified sites are indistinguishable from background rates of conservation. Together, these findings provide strong evidence for the functional importance of highly modified sites. More broadly, our study establishes a methodological framework for the evolutionary analysis of other RNA modifications, emphasizing the importance of accounting for modification biogenesis; employing quantitative and single-nucleotide measurements; and ensuring comprehensive phylogenetic sampling.

Previous SectionNext Section

RESULTS

Assembly of m6A and evolutionary data sets

To assess the extent to which m6A was subject to evolutionary selection, we first analyzed the methylation levels of all expressed (coverage >15) DRACH motifs across the entire Hek293T transcriptome, on the basis of measurements conducted via GLORI, a quantitative and single-nucleotide-based approach for m6A detection (Liu et al. 2022). We defined two groups: “methylated” sites, defined as all sites with methylation levels exceeding 20% (n = 98,959), and “unmethylated controls” with levels below 5% (n = 321,289), controlling for any source of methylation-independent conservation acting on DRACH motifs.

We also established an additional control group of matched non-DRACH sites to control for any gene- or region-specific sources of selection that are independent of m6A or DRACH and may confound the analysis. To ensure we had not inadvertently selected a motif with an unrelated function, we considered a number of different consensus motifs to act as a control that had similar properties to DRACH (see Materials and Methods). After verifying that candidate motifs showed similar distribution and conservation patterns (Supplemental Fig. 1A–C), we ultimately chose GDBYC, the inverse DRACH motif. To assist in disentangling the source of constraints, we only included “matched GDBYC controls”: For each DRACH site, we identified the nearest GDBYC motif within 25 bp and within the same gene region (3′ UTR, etc.) (n = 218,225). In subsequent analyses, we stratified all sites, including these nonmethylatable GDBYC controls, based on methylation level by assigning each control the methylation level of its matched DRACH site.

All sites were additionally annotated by their gene location: 5′ untranslated region (5′ UTR), coding sequence (CDS), or 3′ untranslated region (3′ UTR).

To obtain homologous sequences for each site identified in the human transcriptome, we utilized the Cactus-447 multiple genome alignment. This alignment was constructed by combining the genomes of species collected as part of the Zoonomia project (Zoonomia Consortium 2020) with those from a survey of primate species (Kuderna et al. 2023), resulting in a comprehensive genome alignment comprising 447 placental mammal species. While this broad evolutionary sampling provided substantial statistical power for estimating motif conservation, the data set contains inevitable taxonomic bias (Fig. 1A). Primates account for 233 species, ∼50% of the data set, despite representing only a single mammalian order and a small fraction of mammalian diversity overall. This taxonomic oversampling, resulting from the nonrandom and incomplete species sampling, poses potential risks for phylogenetic inference without appropriate correction (Pybus and Harvey 2000; Cusimano and Renner 2010), as addressed later in this paper.

FIGURE 1.
View larger version:
FIGURE 1.

Establishing measures of motif conservation. (A) Phylogenetic tree from the Cactus-447 data set, comprising two constituent data sets: Zoonomia 241 (black branches) and Primate data set (blue branches). Branch lengths represent substitutions per site, with Homo sapiens as the reference species (red dot). (B,C) Simplified tree topologies illustrating how phylogenetic structure biases motif conservation estimates. Each tree shows parsimony-reconstructed nucleotide substitutions (black ticks) for different permutations of the same motif sequences, demonstrating how relying solely on motif frequency across taxa can lead to overestimation (B) or underestimation (C) of m6A consensus motif (DRACH) conservation. Branch colors indicate the reconstructed motif state along each branch: Red represents DRACH motifs, while blue represents non-DRACH motifs. (D) Comparison of consensus frequency scores (CFSs) between expressed DRACH and control motif sites. (E) Comparison of motif entropy scores (MESs) between expressed DRACH and control motif sites.

Homologous sequences corresponding to all identified m6A sites and control sites in the human transcriptome were extracted from this whole-genome alignment for each species. This process yielded a set of homologous sequences across all species for each site. In cases where homologous sequences were either absent or contained motif-disruption indels, the data were treated as missing, and those species were excluded from the analysis for the affected sites. In all analyses below, we restricted our analysis to sites with homologous sequences identified in at least 304 species, a threshold representing the top 75% of sites based on species coverage (Supplemental Fig. 1D). In total, this data set comprised sequence alignments for 474,226 expressed DRACH motifs including 78,154 unambiguously methylated sites (>20%) and 307,540 unmethylated sites (<5%).

Developing motif evolutionary conservation metrics

To measure evolutionary conservation of these sites, we developed two complementary metrics. Both metrics are calculated from the observed frequencies of the 5-mer sequence at each site to measure variation at the motif level, rather than examining only the central methylated adenosine, as done in the earlier studies (Liu and Zhang 2018).

The first conservation metric, consensus frequency score (CFS), quantifies the evolutionary constraint to retain the consensus motif (DRACH for m6A sites, or the respective motif for controls). For a given site, CFS is the frequency of the consensus motif, calculated as the sum of the frequencies for each motif that forms the consensus (see Materials and Methods).

The second metric, motif entropy score (MES), is based on Shannon entropy, previously used to quantify nucleotide and amino acid conservation (Shenkin et al. 1991; Durbin et al. 1998). Unlike CFS, MES measures overall motif variation at a site without incorporating prior knowledge of consensus preferences and captures the overall evolutionary rate of motif change. For a given site, MES is a negative of the Shannon entropy of the observed motifs across species, normalized such that 1.0 is the mean value within coding regions and 0.0 is the corresponding value across both UTRs (see Materials and Methods). For both CFS and MES, a higher value reflects greater conservation, and conversely, a lower value indicates reduced conservation.

Mitigating phylogenetic bias using branch-weighted frequency estimates

To mitigate bias in conservation metrics caused by the oversampling of primate taxa in the Cactus-447 data set, we weighted our estimates of the motif frequency by branch length, an established approach in assessing regulatory motif conservation (Kheradpour et al. 2007; Friedman et al. 2009). Estimation of motif frequencies from the multiple sequence alignment, without accounting for phylogenetic relationships among taxa, can lead to overestimation (Fig. 1B) or underestimation (Fig. 1C) of the true evolutionary time during which the consensus motif was conserved (which in this case is equivalent to the CFS).

To calculate branch-weighted motif frequencies, we first reconstructed ancestral motifs at internal nodes of the phylogeny for each site using a parsimony-based approach (see Materials and Methods). We then inferred the locations of nucleotide substitutions consistent with the reconstructed ancestral motifs, assuming substitutions occurred at branch midpoints as illustrated in the previous example of phylogenies (Fig. 1B,C). The frequency of a given motif was then calculated as the fraction of total tree length where that motif was present.

The evolutionary conservation metrics, CFS and MES, were then calculated using these estimated frequencies for each site. While these metrics are positively correlated (Supplemental Fig. 1E), the imperfect fit indicates they are capturing distinct aspects of the selective constraint on the motif, because of their contrasting emphases on consensus motif preservation versus overall constraint. Comparing the distributions of these metrics between DRACH sites and controls revealed no significant difference in overall conservation (Fig. 1D,E). Both metrics indicated that the sites are relatively highly conserved, with a median MES close to 1.0, consistent with the majority of sites being located in the CDS (Supplemental Fig. 1B).

Inconsistent differences in conservation metrics across methylated sites and unmethylated controls

We performed an initial set of analyses, aiming to explore whether m6A sites were subject to distinct evolutionary pressures. By separating the DRACH sites into methylated and unmethylated groups and comparing them with matching GDBYC controls, our first analyses revealed no significant differences among groups (Fig. 2A). Unexpectedly, the MES metric indicated that unmethylated sites appeared marginally more conserved than their methylated counterparts (Supplemental Fig. 2A).

FIGURE 2.
View larger version:
FIGURE 2.

Highly methylated m6A sites are the most evolutionarily conserved. (A) Comparison of evolutionary conservation between methylated (m6A+) and unmethylated (m6A−) DRACH sites to matched control (GDBYC) sites using consensus frequency score (CFS) and (B) separated by RNA region. (C) Comparison of CFS of DRACH and control motifs grouped by methylation frequency and (D) separated by RNA region.

Considering the possibility that m6A might be more functional (and hence more conserved) in specific regions within genes, we next stratified the sites by RNA regions (Fig. 2B; Supplemental Fig. 2B). While the distributions remained largely overlapping, the methylated sites in both the UTRs had increased conservation and began to be distinguishable from the unmethylated sites and controls. However, within the CDS, methylated sites appeared to be less conserved than their unmethylated counterparts. Similar to previous studies, the analyses thus far fail to identify a consistent conservation signature at methylated sites.

A major advantage of our data set is the availability of a quantitative assessment of m6A at each site, allowing us to move beyond binary classification between “methylated” and “unmethylated” sites, instead providing a continuous assessment of methylation levels. Accordingly, we next stratified all DRACH sites by their precise methylation levels (Fig. 2C; Supplemental Fig. 2C). This analysis revealed a somewhat complex picture. Among sites methylated at levels >5%, increased methylation levels displayed a clear, dose-dependent positive association with conservation, based both on the CFS and MES metrics (Supplemental Fig. 2C). Such a relationship was lacking among the matched GDBYC controls, where the annotated methylation levels were taken from their corresponding matched DRACH sites.

Nonetheless, perplexingly, unmethylated sites (with levels <5%) were clear outliers in this analysis, displaying unusually high conservation scores, comparable to the effect in the most highly modified sites. Further stratification by RNA region (Fig. 2D) showed that while the m6A-level effect on conservation was maintained across all regions, the anomalous levels of conservation at unmethylated sites were uniquely associated with m6A sites residing within the CDS, but absent from sites in the 5′ and 3′ UTRs. This phenomenon was more prominent in the entropy-based metric (Supplemental Fig. 2D), hinting that this pattern does not result from selective pressure to preserve the motif sequence. This suggests the presence of a confounding variable: an unidentified factor, acting primarily on CDS sites, that simultaneously promotes evolutionary conservation and creates conditions disfavoring methylation.

Distance from exon-junction confounds association between m6A and conservation

Given that similar conservation patterns between highly methylated sites and unmethylated DRACH sites occurred exclusively within the CDS, we hypothesized that proximity to splice junctions may act as a confounder. For the distance between an m6A site and a splice junction to act as a confounder, it would need to be independently associated both with m6A occurrence and with conservation. Indeed, m6A sites tend to be much more distant from splice sites than unmethylated controls (Fig. 3A), consistent with recent findings that m6A sites are depleted near splice junction via the exon junction complex (EJC) (Yang et al. 2022; He et al. 2023; Luo et al. 2023; Uzonyi et al. 2023). In parallel, considering only the GDBYC control sites (and thereby ruling out any effect acting on methylation motifs), we found that sites in close proximity to the splice junction tend to be considerably more conserved than more distant sites (Fig. 3B). These findings are consistent with prior findings and with the interpretation that these splice-proximal regions are enriched in functional elements regulating splicing (Goren et al. 2006). These results thus suggested that the increased conservation levels associated with unmethylated sites may have been an inadvertent consequence of the relative enrichment of these unmethylated sites near splice junctions.

FIGURE 3.
View larger version:
FIGURE 3.

Evolutionary conservation of methylation is distinguishable from proximity to splice junction effect. (A) Distance distributions from methylated (m6A+) and unmethylated (m6A−) DRACH sites to their nearest exon junction complex (EJC). (B) Median consensus frequency score (CFS) of matched control (GDBYC) sites grouped by distance to the nearest EJC. (C,D) Conservation heat maps showing CFS for (C) DRACH sites and (D) matched control sites grouped by methylation frequency and EJC proximity. Color intensity represents conservation strength, with darker red indicating higher conservation. Numbers within cells show the count of sites per group.

To control for this confounder, we jointly stratified sites based both on EJC distance and methylation level (Fig. 3C). This analysis highlighted a clear, dose-dependent positive association between proximity to the EJC and conservation levels independently of methylation levels and—conversely—a clear, dose-dependent positive association between methylation levels and conservation, independently of distance to the EJC. The previously observed anomalous high conservation effect in unmethylated sites was now eliminated, confirming that we had correctly identified the confounding variable. These independent effects were consistently detected using both conservation metrics we developed (Supplemental Fig. 2E). Moreover, applying this stratification approach only to the matched GDBYC motifs (which retain distances to splice junctions, but are not subject to methylation), we confirmed that the conservation effect associated with EJC proximity remained robust, whereas the associations with stoichiometry were abolished (Fig. 3D; Supplemental Fig. 2F).

As further validation, we examined evolutionary conservation in single-exon genes, which lack introns and therefore EJCs (Supplemental Fig. 3A,B). Despite the limited number of sites, we observed a consistent result: The methylation-dependent effect remained present, while the anomalous high conservation in unmethylated sites was absent. Together, these dose-dependent associations thus strongly suggest that m6A modification stoichiometry directly impacts sequence conservation patterns, rather than acting through secondary gene- or region-specific conservation differences.

How strong is the selection on m6A sites?

After accounting for EJC proximity, we observed that highly methylated sites (80%–100%) showed an ∼7% increase in CFS compared to unmethylated sites (<5%) with a corresponding 0.3 unit increase in MES. This conservation effect was diminished near EJCs, likely because these regions are already subject to strong selective constraints, limiting the dynamic range available for increased motif conservation. Given the scaling of the MES metric, we can interpret that the magnitude of this methylation-associated conservation is roughly one-third of the typical conservation difference observed between UTR and coding sequence regions.

Limitations in detecting conservation signatures of RNA modification

While we demonstrated a clear dose-dependent relationship between m6A and evolutionary conservation, this study also highlights the limitations of assessing conservation signatures at modification sites. In our study, we had access to quantitative measurements, to measurements across tens of thousands of sites, and we considered the full consensus motif. Prior studies lacked access to these and did not observe the consistent association between methylation and conservation that we observe. As a guide to future studies (either focused on m6A or on other modifications), we sought to understand which of these differences contributed most to our ability to observe the association between methylation and conservation.

We first investigated measurement inaccuracy—a common limitation in sequencing approaches, and particularly pertinent when single-nucleotide measurements are not available—by randomly reassigning methylation levels to subsets of sites (see Materials and Methods) and measuring how frequently we could still detect Spearman correlation coefficients between average motif conservation and methylation levels, excluding the unmethylated category (0.0–0.05), above specified thresholds. We intentionally avoided using our full data set because of the previously established EJC bias toward unmethylated and weakly methylated sites, which would obscure the relationship between higher methylation levels and increased evolutionary conservation. Notably, the methylation–conservation association persisted until methylation group misassignment exceeded 50% (Fig. 4A), suggesting that the underlying signal remains detectable even with substantial additional measurement noise. Overall, these results indicate that improved modification level measurements may be the least critical factor among the methodological advances that enabled detection of the positive association between m6A methylation and evolutionary conservation in this study.

FIGURE 4.
View larger version:
FIGURE 4.

Effects of limited statistical power on detecting correlation between site methylation frequency and evolutionary selection. (A,B) Detection frequency of significant Spearman correlations between methylation groups (>5% methylation) and DRACH conservation under different experimental conditions: (B) with varying proportions of sites having inaccurate methylation measurements and (C) with different sample sizes. Analysis used correlation thresholds of 0.5 and 0.8, with each condition resampled 1000 times to estimate detection probability. (C) Evolutionary conservation of the central nucleotide in DRACH sites binned by methylation frequency, compared to matched control (GDBYC) sites. For DRACH sites, this represents the central adenosine conservation; for control sites, whatever central nucleotide is present at the site in humans is used as the reference.

We next evaluated how limited site numbers affect detection of the conservation–methylation relationship (Fig. 4B). Using a similar approach, we randomly sampled subsets of methylated m6A sites and tested whether correlations above specified thresholds could still be detected across the reduced data sets (see Materials and Methods). We observed that detecting a correlation between evolutionary conservation and methylation levels comparable to our full data set results (ρ = 0.94) likely requires at least 10,000 modified sites. This sample size threshold explains both why prior studies that were limited by fewer known sites at the time failed to detect this evolutionary signal and suggests that this overall approach to identifying evolutionary constraints may be underpowered for other RNA modifications, such as m5C (Ma et al. 2022), whose transcriptome-wide distribution is far more limited.

Finally, we evaluated the impact of restricting the conservation metric to only the central adenosine residue, rather than considering the entire m6A consensus motif. This single-base approach resembles the approach utilized by Liu and Zhang (2018); however, rather than conducting a pairwise comparison between only mouse and human, we measured the fraction of adenosines conserved across all taxa in our data set (see Materials and Methods). We used the same set of m6A sites and neighboring control sites as before: those located within human DRACH and GDBYC motifs, respectively. Since the GDBYC control motif permits any nucleotide except adenosine at the central position, we quantified evolutionary conservation based on the specific nucleotide present at that site in the human reference. We found that considering only the methylation site reduced the association between methylation and conservation (Fig. 4C). Qualitatively, a similar trend was captured as when considering the entirety of the motif, but its dynamic range was compressed (compared to Fig. 2C).

Previous SectionNext Section

DISCUSSION

In this study, we quantify the evolutionary pressures on m6A sites to reveal whether the deep conservation of the m6A machinery reflects methylation-dependent versus methylation-independent functions. Furthermore, by contrasting conservation profiles across gene regions, we aim to pinpoint loci where m6A exerts stronger functional impact. Our findings identify a dose-dependent association between methylation levels and conservation across mammals, highly suggestive of functionality of at least a subset of methylated sites. The conservation signature is apparent only for highly methylated m6A sites (>60%) in both the CDS and in 3′ UTRs, with the degree of conservation being roughly one-third the magnitude of the conservation difference typically seen between UTR and coding sequence (CDS) regions. These results support a model where m6A deposition at target sites imparts a consistent and evolutionarily selected function.

Our results challenge those in previous studies (Liu and Zhang 2018; Zhang et al. 2020b) that either found no evidence for conservation or found mixed effects depending on the whereabouts of the modification. We attribute this discrepancy primarily to several methodological improvements: our larger sample size with quantitative m6A measurements; the substantially broader evolutionary tree we surveyed; our analysis of the entire DRACH motif rather than just the modified site; and our improved understanding of m6A biogenesis and its dependence on exon junction proximity, which allowed us to stratify methylation sites based on their distance to the nearest splice sites.

Recently, a series of manuscripts reported that m6A primarily triggers mRNA decay when present in the CDS, but not in the 3′ UTR (Zhou et al. 2024; Linder et al. 2025; Murakami et al. 2025), implying the UTR sites may be incidental. Here we observe clear selection acting on m6A sites across all regions. The methylation-associated effect in 3′ UTRs even appears stronger than the CDS, although this comparison requires cautious interpretation since the inherently high baseline conservation of coding regions may limit the dynamic range available to measure additional conservation effects. Together, these findings may hint that m6A plays distinct but equally important functional roles across CDS and UTR regions.

Despite our observation that more highly methylated sites tend to be, on average, more conserved, using conservation patterns in order to prioritize functional methylation sites remains a challenge. The evolutionary conservation level at any given site is predominantly determined by factors other than its methylation status, including the genetic code, splice signals, and binding sites to diverse RNA binding proteins. Consequently, discovering an individual highly conserved m6A site provides little information pertaining to the likelihood of the methylation playing a functional role. To identify the most functionally relevant m6A sites, future work will need to identify modifications that exhibit greater conservation than expected based on their local sequence context alone. This will require developing statistical models that predict baseline conservation scores for unmethylated sites in equivalent contexts, which can then serve as null expectations.

An additional limitation of our study is its anthropocentric focus. Specifically, the candidate sites analyzed were identified exclusively from human transcriptome data, and the evolutionary data set was moreover biased toward primate species. Access to high-quality m6A measurements in additional species will allow us to extend this exploration of m6A conservation.

Our study further has important bearings on evolutionary analyses directed at additional RNA modifications. While m6A is the most abundant internal mRNA modification and benefits from robust detection methods, most other mRNA modifications are far less frequent and are often detected with lower resolution, lower sensitivity, or less quantitative accuracy (Wiener and Schwartz 2021). Moreover, the sequence and structural preferences guiding the installation of these modifications remain less well-defined, limiting our ability to identify appropriate motif-based controls or correct for context-dependent confounders, as we did for m6A. These challenges can lead to false negatives in evolutionary studies, where lack of apparent conservation may result from inadequate resolution, insufficient sample size, or failure to account for positional biases—rather than indicating true evolutionary neutrality of the modification. Our findings thus underscore the need for caution when interpreting evolutionary analyses of RNA modifications. Specifically, they argue that evolutionary neutrality cannot be confidently inferred from the absence of conservation unless the analysis meets key criteria: quantitative measurement of modification stoichiometry, sufficient number of high-confidence modified sites, robust phylogenetic sampling, and consideration of sequence and structural determinants of modification placement. By establishing these requirements, our study offers a framework for future analyses of RNA modifications beyond m6A and suggests that with the right data sets and analytical approaches, it may be possible to uncover previously hidden signatures of selection acting on other marks as well.

In summary, our study provides strong evidence that m6A sites with high stoichiometry are under evolutionary selection, revealing a clear, dose-dependent relationship between methylation level and sequence conservation. This finding argues against models in which m6A is largely nonfunctional, and instead supports a view in which m6A deposition contributes to gene regulation at a broad scale, with functional consequences encoded at the level of individual sites. Together, our results highlight the power of evolutionary analysis as a tool for probing the functional relevance of epitranscriptomic marks—but also the necessity of integrating quantitative measurements, prior knowledge of modification biogenesis, and rigorous control for confounding variables.

Previous SectionNext Section

MATERIALS AND METHODS

Site discovery in human

Filtering and compiling GLORI-seq reads

To construct our data set of GLORI-treated raw RNA reads, we combined the two control replicates from the original GLORI study (Liu et al. 2022). These reads were aligned to the hg38 human reference genome using Hisat-3n (Zhang et al. 2021), which corrects for the A-to-G conversions induced by GLORI treatment. After removing all unmapped reads, we applied a stringent filtering step to all reads with more than three unconverted sites, as these likely reflect incomplete treatment. This unconverted site filter removed only ∼0.03% of mapped reads.

Identifying consensus motif positions

In addition to the m6A consensus motif DRACH, we analyzed seven comparable control consensus motifs. Like DRACH, each of these consensus motifs contains 18 possible variants. The controls include GDBYC (the inverse of DRACH), the reverse complement of each motif in DRACH (rc_DRACH), and five randomly selected consensus motifs (DCDAR, DTYCH, HGRTH, HTYCH, and RGHTD) that matched DRACH in key properties like GC content while maintaining minimal submotif overlap. All the controls (Supplemental Fig. 1C,D) showed similar overall MES conservation levels but differing CFS levels; we focused on GDBYC as our control since it most closely resembled DRACH, though the choice of control did not affect the overall results (data not shown).

For every consensus motif, we identified all genomic positions (relative to hg38) that showed expression levels of ≥15 reads in the GLORI RNA-seq data set outlined above.

Calling m6A sites

While identifying motif positions in the previous step, we simultaneously quantified methylation frequency for the m6A consensus motif DRACH. This frequency was calculated as the ratio of modified sites (detected via A-to-G conversions in GLORI-seq data) to the total read count at the central adenosine position of each DRACH motif. For pairwise analyses, we classified positions as modified (>20%) or unmodified (<5%).

Annotation of motif sites

We annotated each identified motif position using the RefSeq annotations (Pruitt et al. 2007), to include gene, region (3′ UTR, CDS, and 5′ UTR) as well as the nearest exon junction complex (EJC) boundary. Only motifs located in mRNA and nonintronic regions were considered in the analysis.

To determine the nearest EJC boundary for each motif, we scanned the gene isoform annotations for the closest upstream and downstream boundaries (excluding transcription start and end sites). The shortest distance between the motif's central nucleotide and the nearest boundary was then calculated. For genes with multiple isoforms exhibiting differing boundaries, we used the mean distance across all isoforms.

m6A matched controls

To account for gene-level confounding factors, we only included control sites that were paired with the nearest DRACH motif (within 25 bp) in the same gene region. Unmatched control motifs were excluded from further analysis. For negative control analyses, each control site, despite being nonmethylatable, was annotated with the methylation level of its matched DRACH site.

Building the phylogenetic data set

We extracted homologous motifs corresponding to the previously identified human sites from the Cactus-447 whole-genome alignment data set (multiple alignment format, [MAF]). This data set combines two sources: Shao et al. (2023) and Kuderna et al. (2023). For each site, we identified alignment columns containing homologs relative to the human reference sites and extracted all corresponding motifs. We excluded positions that spanned multiple alignment blocks in the MAF file, though such cases were rare as our analysis focused on relatively well-conserved genes that typically maintain continuous alignment blocks. Taxa lacking homologous sites were treated as missing data.

Minimal species coverage

We excluded homologous motifs with gaps or deletions relative to the human reference. We retained only sites that we covered at least 304 species (top 75% by species coverage) to minimize bias.

Estimating conservation metrics

Both conservation metrics used in this analysis are derived from the branch-weighted frequency of each motif across the phylogeny, equivalent to the relative evolutionary time that each motif is present across the tree topology. This makes the reasonable assumption that aligned sites have a single common ancestor. We used the phylogenetic tree provided with the Cactus-447, constructed as described in Kuderna et al. (2023).

Estimating branch-length weighted motif frequency

To estimate the motif frequency across each branch of the phylogenetic tree, it is necessary to reconstruct the ancestral motif at every internal node. This is achieved using a crude parsimony-based approach, which recursively infers the motif state of each internal node by tracing backward from the known motif sequences of the taxa at the tips of the tree toward the root. While parsimony-based approaches have acknowledged limitations for ancestral state reconstruction (Randall et al. 2016), we consider this method adequate for our purposes, as the inferred ancestral motifs need only be sufficiently accurate to enable the estimation of overall motif frequencies.

For each site, ancestral reconstruction is performed by independently reconstructing the nucleotide states at each position within the motif. The motif sequence at each internal node is then determined by concatenating the reconstructed nucleotide states from all positions within the motif.

Taxa with sequences containing gaps or lacking a homologous sequence at the position identified in humans are treated as missing data. Internal nodes where all descendants have missing sequences are also classified as missing, and no sequence reconstruction is attempted for such nodes.

The reconstruction of nucleotide states at each site begins at the tip taxa and proceeds recursively backward up the tree. Each internal node is visited only after the states of all its descendants have been inferred. Provided not all descendants are missing sequences, the nucleotide state of an internal node is assigned as the most frequently observed nucleotide among the tip taxa within the monophyletic group defined by that node. In the event of a tie, the nucleotide that is most common across the remaining taxa in the tree is selected; if the tie persists, the nucleotide is chosen at random. Note that the states of previously inferred internal nodes do not influence the determination of subsequent nucleotide states.

With the motif sequence inferred for each node in the tree, we can calculate the motif frequency across the tree by focusing on the set of branches, Bi, for each site i, where a complete sequence is present at both the ancestral and descendant nodes of the branch. Using Equation 1 we can calculate the frequency, Formula of motif m at site i, where lb is the length of branch b. Equation 1 utilized an indicator function (Eq. 2) to specify whether motif, m is the ancestral (Anc) or descendant state (Dec) of branch bi at site i.

When a branch exhibits different ancestral and descendant motif states, the transition is assumed to occur at the midpoint of the branch. As a result, half of the branch length is assigned to each of the observed motifs.Formula (1) Formula (2) Note that the denominator of Equation 1 will vary between sites because different combinations of taxa may have missing data. As a result, distinct subtrees will effectively be considered for each site.

Consensus frequency score (CFS)

The consensus frequency score (CFS) metric quantifies the relative frequency of the consensus motif, C, observed at site i across the phylogenetic tree, as described by Equation 3. Here, the set C represents the collection of motifs that form the consensus motif. For example, in the case of DRACH, the set C would be defined as { AAACT, GAACT … }.Formula (3)

Motif entropy score (MES)

The motif entropy score (MES) metric quantifies the motif variation observed at site i across the phylogenetic tree, independent of any specific consensus motif. It is derived from the Shannon motif entropy, Hi at a given locus (shown in Eq. 4) and transformed to assist interpretability. This transformation (shown in Eq. 5) involves taking the negative of the Shannon entropy at each site, ensuring that a site with low entropy receives a high score, indicating greater conservation. Additionally, the Shannon entropy is normalized by the median entropy of control GDBYC loci located in UTR, Formula and CDS regions, Formula. As a result, a score of 0 corresponds to the average conservation level in the UTR region, while a score of 1 represents the average conservation level in a coding region.Formula (4) Formula (5)

Central nucleotide conservation

When testing the case where the consensus motif is unknown, we used percentage nucleotide conservation. This is calculated in the same manner as CFS using branch-length weighted frequencies, but rather than reconstructing the motif, we simply take the frequency of the relevant nucleotide at the central position. For DRACH sites, this represents central adenosine conservation; for control sites, the central nucleotide present at the site in humans is used as the reference.

Power analysis with limited data

For measurement inaccuracy

We randomly reassigned methylation levels to different proportions of sites across the data set. For selected sites, the new “incorrect” methylation values were drawn randomly from the observed distribution of methylated sites (ranging from 0.05 to 1.0). We excluded the unmethylated site group (0.0–0.05) from correlation calculations and conducted 1000 independent iterations at each error rate. Following the same approach as the minimum site detection analysis, we calculated Spearman correlation coefficients between mean motif conservation scores and methylation levels and then determined the proportion of iterations that yielded correlations exceeding the specified thresholds (ρ > 0.5 and ρ > 0.8).

For minimum site detection

We randomly sampled subsets of methylated m6A sites, excluding unmethylated sites (<0.05), in groups of 100, 1000, 2500, 5000, 7500, and 10,000 sites, performing 1000 independent sampling iterations for each subset size. For each sampled subset, we calculated Spearman correlation coefficients between mean motif conservation scores and methylation levels and then determined the proportion of iterations that yielded correlations exceeding specified thresholds (ρ > 0.5 and ρ > 0.8).

Previous SectionNext Section

DATA DEPOSITION

The whole-genome alignment data underlying this article are publicly available from the Cactus-447 repository at https://hgdownload.soe.ucsc.edu/goldenPath/hg38/cactus447way/.

Previous SectionNext Section

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.

Previous SectionNext Section

ACKNOWLEDGMENTS

S.S. would like to acknowledge funding by the ISF (2291/25) and by DIP (512/1-1).

Previous SectionNext Section

Footnotes

  • Received November 11, 2025.
  • Accepted December 9, 2025.

This article, published in RNA, is available under a Creative Commons License (Attribution 4.0 International), as described at http://creativecommons.org/licenses/by/4.0/.

Previous SectionNext Section

REFERENCES

  1. ↵
    Brodersen P, Arribas-Hernández L. 2024. The m6A-YTH regulatory system in plants: a status. Curr Opin Plant Biol 82: 102650. doi:10.1016/j.pbi.2024.102650
    CrossRefMedlineGoogle Scholar
  2. ↵
    Cusimano N, Renner SS. 2010. Slowdowns in diversification rates from real phylogenies may not be real. Syst Biol 59: 458–464. doi:10.1093/sysbio/syq032
    CrossRefMedlineWeb of ScienceGoogle Scholar
  3. ↵
    Dagan Y, Yesharim Y, Bonneau AR, Frankovits T, Schwartz S, Reddien PW, Wurtzel O. 2022. m6A is required for resolving progenitor identity during planarian stem cell differentiation. EMBO J 41: e109895. doi:10.15252/embj.2021109895
    CrossRefMedlineGoogle Scholar
  4. ↵
    Dierks D, Garcia-Campos MA, Uzonyi A, Safra M, Edelheit S, Rossi A, Sideri T, Varier RA, Brandis A, Stelzer Y, et al. 2021. Multiplexed profiling facilitates robust m6A quantification at site, gene and sample resolution. Nat Methods 18: 1060–1067. doi:10.1038/s41592-021-01242-z
    CrossRefGoogle Scholar
  5. ↵
    Du H, Zhao Y, He J, Zhang Y, Xi H, Liu M, Ma J, Wu L. 2016. YTHDF2 destabilizes m6A-containing RNA through direct recruitment of the CCR4-NOT deadenylase complex. Nat Commun 7: 12626. doi:10.1038/ncomms12626
    CrossRefMedlineGoogle Scholar
  6. ↵
    Durbin R, Eddy SR, Krogh A, Mitchison G. 1998. Biological sequence analysis: probabilistic models of proteins and nucleic acids. Cambridge University Press, Cambridge.
    Google Scholar
  7. ↵
    Ensinck I, Maman A, Albihlal WS, Lassandro M, Salzano G, Sideri T, Howell SA, Calvani E, Patel H, Bushkin G, et al. 2023. The yeast RNA methylation complex consists of conserved yet reconfigured components with m6A-dependent and independent roles. eLife 12: RP87860. doi:10.7554/eLife.87860
    CrossRefMedlineGoogle Scholar
  8. ↵
    Friedman RC, Farh KK-H, Burge CB, Bartel DP. 2009. Most mammalian mRNAs are conserved targets of microRNAs. Genome Res 19: 92–105. doi:10.1101/gr.082701.108
    Abstract/FREE Full Text
  9. ↵
    Goren A, Ram O, Amit M, Keren H, Lev-Maor G, Vig I, Pupko T, Ast G. 2006. Comparative analysis identifies exonic splicing regulatory sequences—the complex definition of enhancers and silencers. Mol Cell 22: 769–781. doi:10.1016/j.molcel.2006.05.008
    CrossRefMedlineWeb of ScienceGoogle Scholar
  10. ↵
    He PC, Wei J, Dou X, Harada BT, Zhang Z, Ge R, Liu C, Zhang L-S, Yu X, Wang S, et al. 2023. Exon architecture controls mRNA m6A suppression and gene expression. Science 379: 677–682. doi:10.1126/science.abj9090
    CrossRefMedlineGoogle Scholar
  11. ↵
    Higuchi M, Maas S, Single FN, Hartner J, Rozov A, Burnashev N, Feldmeyer D, Sprengel R, Seeburg PH. 2000. Point mutation in an AMPA receptor gene rescues lethality in mice deficient in the RNA-editing enzyme ADAR2. Nature 406: 78–81. doi:10.1038/35017558
    CrossRefMedlineWeb of ScienceGoogle Scholar
  12. ↵
    Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, Yang C, Chen Y. 2021. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Ther 6: 74. doi:10.1038/s41392-020-00450-x
    CrossRefGoogle Scholar
  13. ↵
    Kheradpour P, Stark A, Roy S, Kellis M. 2007. Reliable prediction of regulator targets using 12 Drosophila genomes. Genome Res 17: 1919–1931. doi:10.1101/gr.7090407
    Abstract/FREE Full Text
  14. ↵
    Kuderna LFK, Gao H, Janiak MC, Kuhlwilm M, Orkin JD, Bataillon T, Manu S, Valenzuela A, Bergman J, Rousselle M, et al. 2023. A global catalog of whole-genome diversity from 233 primate species. Science 380: 906–913. doi:10.1126/science.abn7829
    CrossRefMedlineGoogle Scholar
  15. ↵
    Lence T, Akhtar J, Bayer M, Schmid K, Spindler L, Ho CH, Kreim N, Andrade-Navarro MA, Poeck B, Helm M, et al. 2016. m6A modulates neuronal functions and sex determination in Drosophila. Nature 540: 242–247. doi:10.1038/nature20568
    CrossRefMedlineGoogle Scholar
  16. ↵
    Linder B, Sharma P, Wu J, Birbaumer T, Eggers C, Murakami S, Ott RE, Fenzl K, Vorgerd H, Erhard F, et al. 2025. tRNA modifications tune m6A-dependent mRNA decay. Cell 188: P3715–P3727.e13. doi:10.1016/j.cell.2025.04.013
    CrossRefGoogle Scholar
  17. ↵
    Liu Z, Zhang J. 2018. Most m6A RNA modifications in protein-coding regions are evolutionarily unconserved and likely nonfunctional. Mol Biol Evol 35: 666–675. doi:10.1093/molbev/msx320
    CrossRefMedlineGoogle Scholar
  18. ↵
    Liu C, Sun H, Yi Y, Shen W, Li K, Xiao Y, Li F, Li Y, Hou Y, Lu B, et al. 2022. Absolute quantification of single-base m6A methylation in the mammalian transcriptome using GLORI. Nat Biotechnol 41: 355–366. doi:10.1038/s41587-022-01487-9
    CrossRefMedlineGoogle Scholar
  19. ↵
    Luo Z, Ma Q, Sun S, Li N, Wang H, Ying Z, Ke S. 2023. Exon-intron boundary inhibits m6A deposition, enabling m6A distribution hallmark, longer mRNA half-life and flexible protein coding. Nat Commun 14: 4172. doi:10.1038/s41467-023-39897-1
    CrossRefMedlineGoogle Scholar
  20. ↵
    Ma J, Song B, Wei Z, Huang D, Zhang Y, Su J, de Magalhães JP, Rigden DJ, Meng J, Chen K. 2022. m5C-Atlas: a comprehensive database for decoding and annotating the 5-methylcytosine (m5C) epitranscriptome. Nucleic Acids Res 50: D196–D203. doi:10.1093/nar/gkab1075
    CrossRefMedlineGoogle Scholar
  21. ↵
    Murakami S, Olarerin-George AO, Liu JF, Zaccara S, Hawley B, Jaffrey SR. 2025. m6A alters ribosome dynamics to initiate mRNA degradation. Cell 188: 3728–3743.e20. doi:10.1016/j.cell.2025.04.020
    CrossRefGoogle Scholar
  22. ↵
    Perlegos AE, Shields EJ, Shen H, Liu KF, Bonini NM. 2022. Mettl3-dependent m6A modification attenuates the brain stress response in Drosophila. Nat Commun 13: 5387. doi:10.1038/s41467-022-33085-3
    CrossRefMedlineGoogle Scholar
  23. ↵
    Porman AM, Roberts JT, Duncan ED, Chrupcala ML, Levine AA, Kennedy MA, Williams MM, Richer JK, Johnson AM. 2022. A single N6-methyladenosine site regulates lncRNA HOTAIR function in breast cancer cells. PLoS Biol 20: e3001885. doi:10.1371/journal.pbio.3001885
    CrossRefMedlineGoogle Scholar
  24. ↵
    Pruitt KD, Tatusova T, Maglott DR. 2007. NCBI reference sequences (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res 35: D61–D65. doi:10.1093/nar/gkl842
    CrossRefMedlineWeb of ScienceGoogle Scholar
  25. ↵
    Pybus OG, Harvey PH. 2000. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc Biol Sci 267: 2267–2272. doi:10.1098/rspb.2000.1278
    CrossRefMedlineWeb of ScienceGoogle Scholar
  26. ↵
    Randall RN, Radford CE, Roof KA, Natarajan DK, Gaucher EA. 2016. An experimental phylogeny to benchmark ancestral sequence reconstruction. Nat Commun 7: 12847. doi:10.1038/ncomms12847
    CrossRefMedlineGoogle Scholar
  27. ↵
    Roignant J-Y, Soller M. 2017. m6A in mRNA: an ancient mechanism for fine-tuning gene expression. Trends Genet 33: 380–390. doi:10.1016/j.tig.2017.04.003
    CrossRefMedlineGoogle Scholar
  28. ↵
    Schwartz S, Agarwala SD, Mumbach MR, Jovanovic M, Mertins P, Shishkin A, Tabach Y, Mikkelsen TS, Satija R, Ruvkun G, et al. 2013. High-resolution mapping reveals a conserved, widespread, dynamic mRNA methylation program in yeast meiosis. Cell 155: 1409–1421. doi:10.1016/j.cell.2013.10.047
    CrossRefMedlineWeb of ScienceGoogle Scholar
  29. ↵
    Shachar R, Dierks D, Garcia-Campos MA, Uzonyi A, Toth U, Rossmanith W, Schwartz S. 2024. Dissecting the sequence and structural determinants guiding m6A deposition and evolution via inter- and intra-species hybrids. Genome Biol 25: 48. doi:10.1186/s13059-024-03182-1
    CrossRefMedlineGoogle Scholar
  30. ↵
    Shao Y, Zhou L, Li F, Zhao L, Zhang B-L, Shao F, Chen J-W, Chen C-Y, Bi X, Zhuang X-L, et al. 2023. Phylogenomic analyses provide insights into primate evolution. Science 380: 913–924. doi:10.1126/science.abn6919
    CrossRefMedlineGoogle Scholar
  31. ↵
    Shenkin PS, Erman B, Mastrandrea LD. 1991. Information-theoretical entropy as a measure of sequence variability. Proteins 11: 297–313. doi:10.1002/prot.340110408
    CrossRefMedlineWeb of ScienceGoogle Scholar
  32. ↵
    Uzonyi A, Dierks D, Nir R, Kwon OS, Toth U, Barbosa I, Burel C, Brandis A, Rossmanith W, Le Hir H, et al. 2023. Exclusion of m6A from splice-site proximal regions by the exon junction complex dictates m6A topologies and mRNA stability. Mol Cell 83: 237–251.e7. doi:10.1016/j.molcel.2022.12.026
    CrossRefMedlineGoogle Scholar
  33. ↵
    Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, Fu Y, Parisien M, Dai Q, Jia G, et al. 2014. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature 505: 117–120. doi:10.1038/nature12730
    CrossRefMedlineWeb of ScienceGoogle Scholar
  34. ↵
    Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, Weng X, Chen K, Shi H, He C. 2015. N6-methyladenosine modulates messenger RNA translation efficiency. Cell 161: 1388–1399. doi:10.1016/j.cell.2015.05.014
    CrossRefMedlineGoogle Scholar
  35. ↵
    Wang C-X, Cui G-S, Liu X, Xu K, Wang M, Zhang X-X, Jiang L-Y, Li A, Yang Y, Lai W-Y, et al. 2018. METTL3-mediated m6A modification is required for cerebellar development. PLoS Biol 16: e2004880. doi:10.1371/journal.pbio.2004880
    CrossRefMedlineGoogle Scholar
  36. ↵
    Wang G, Li H, Ye C, He K, Liu S, Jiang B, Ge R, Gao B, Wei J, Zhao Y, et al. 2024. Quantitative profiling of m6A at single base resolution across the life cycle of rice and Arabidopsis. Nat Commun 15: 4881. doi:10.1038/s41467-024-48941-7
    CrossRefMedlineGoogle Scholar
  37. ↵
    Wiener D, Schwartz S. 2021. The epitranscriptome beyond m6A. Nat Rev Genet 22: 119–131. doi:10.1038/s41576-020-00295-8
    CrossRefGoogle Scholar
  38. ↵
    Xin T, Zhang Z, Zhang Y, Li X, Wang S, Wang G, Li H, Wang B, Zhang M, Li W, et al. 2025. Recessive epistasis of a synonymous mutation confers cucumber domestication through epitranscriptomic regulation. Cell 188: 4517–4529.e15. doi:10.1016/j.cell.2025.06.007
    CrossRefMedlineGoogle Scholar
  39. ↵
    Xiong X, Hou L, Park YP, Molinie B, GTEx Consortium, Gregory RI, Kellis M. 2021. Genetic drivers of m6A methylation in human brain, lung, heart and muscle. Nat Genet 53: 1156–1165. doi:10.1038/s41588-021-00890-3
    CrossRefMedlineGoogle Scholar
  40. ↵
    Xu W, Li J, He C, Wen J, Ma H, Rong B, Diao J, Wang L, Wang J, Wu F, et al. 2021. METTL3 regulates heterochromatin in mouse embryonic stem cells. Nature 591: 317–321. doi:10.1038/s41586-021-03210-1
    CrossRefMedlineGoogle Scholar
  41. ↵
    Yang X, Triboulet R, Liu Q, Sendinc E, Gregory RI. 2022. Exon junction complex shapes the m6A epitranscriptome. Nat Commun 13: 7904. doi:10.1038/s41467-022-35643-1
    CrossRefMedlineGoogle Scholar
  42. ↵
    Yue H, Nie X, Yan Z, Weining S. 2019. N6-methyladenosine regulatory machinery in plants: composition, function and evolution. Plant Biotechnol J 17: 1194–1208. doi:10.1111/pbi.13149
    CrossRefMedlineGoogle Scholar
  43. ↵
    Zaccara S, Jaffrey SR. 2020. A unified model for the function of YTHDF proteins in regulating m6A-modified mRNA. Cell 181: 1582–1595.e18. doi:10.1016/j.cell.2020.05.012
    CrossRefMedlineGoogle Scholar
  44. ↵
    Zhang H, Shi X, Huang T, Zhao X, Chen W, Gu N, Zhang R. 2020a. Dynamic landscape and evolution of m6A methylation in human. Nucleic Acids Res 48: 6251–6264. doi:10.1093/nar/gkaa347
    CrossRefMedlineGoogle Scholar
  45. ↵
    Zhang Z, Luo K, Zou Z, Qiu M, Tian J, Sieh L, Shi H, Zou Y, Wang G, Morrison J, et al. 2020b. Genetic analyses support the contribution of mRNA N6-methyladenosine (m6A) modification to human disease heritability. Nat Genet 52: 939–949. doi:10.1038/s41588-020-0644-z
    CrossRefMedlineGoogle Scholar
  46. ↵
    Zhang Y, Park C, Bennett C, Thornton M, Kim D. 2021. Rapid and accurate alignment of nucleotide conversion sequencing reads with HISAT-3N. Genome Res 31: 1290–1295. doi:10.1101/gr.275193.120
    Abstract/FREE Full Text
  47. ↵
    Zhou Y, Ćorović M, Hoch-Kraft P, Meiser N, Mesitov M, Körtel N, Back H, Naarmann-de Vries IS, Katti K, Obrdlík A, et al. 2024. m6A sites in the coding region trigger translation-dependent mRNA decay. Mol Cell 84: 4576–4593.e12. doi:10.1016/j.molcel.2024.10.033
    CrossRefMedlineGoogle Scholar
  48. ↵
    Zoonomia Consortium. 2020. A comparative genomics multitool for scientific discovery and conservation. Nature 587: 240–245. doi:10.1038/s41586-020-2876-6
    CrossRefMedlineGoogle Scholar
Previous Section
 

MEET THE FIRST AUTHOR

Hamish N.C. Pike
View larger version:
Hamish N.C. Pike

Meet the First Author(s) is an editorial feature within RNA, in which the first author(s) of research-based papers in each issue have the opportunity to introduce themselves and their work to readers of RNA and the RNA research community. Hamish N.C. Pike is the first author of this paper, “High-stoichiometry m6A sites are evolutionarily conserved.” Hamish is currently a postdoctoral researcher in Dr. Schraga Schwartz's lab, which investigates the RNA epitranscriptome and how various chemical base modifications regulate the RNA life cycle. Hamish completed his PhD with Dr. David Pollock at the University of Colorado Anschutz Medical Campus, where his research centered on understanding the mechanisms underlying protein evolution. Specifically, Hamish developed statistical methods to quantify patterns of amino acid substitution.

What are the major results described in your paper and how do they impact this branch of the field?

Specialized sequencing technologies developed over the past decade have enabled transcriptome-wide mapping of m6A—the methylation of adenosine in RNA—identifying hundreds of thousands of sites and revealing key factors determining its distribution: a specific sequence motif, and exclusion from the 100 nucleotides around exon junction complexes. Yet, despite being associated with numerous regulatory processes, particularly mRNA stability, deducing the function of individual m6A sites has proven elusive. It has even been proposed that individual m6A sites lack any function, with observed correlations representing mere coincidence: Artifacts of the transcriptional and splicing conditions present when m6A writers are active. Our work addresses this debate by showing that high-stoichiometry m6A sites are evolutionarily conserved, providing evidence that at least some m6A modifications are genuinely functional.

What led you to study RNA or this aspect of RNA science?

RNA biology is compelling as it offers technically challenging yet biologically relevant problems that I enjoy solving. RNA is not only the intermediate between DNA and protein in the central dogma, but also a molecule that combines the distinctive features of both. Like DNA, RNA is tractable: Its familiar sequence can be interrogated using the same high-throughput sequencing technologies to generate snapshots of cellular state. Yet unlike DNA, which serves an exclusively symbolic role, RNA can function more like proteins too: participating directly in the metabolic action and exhibiting diverse enzymatic and regulatory behaviors. As a computational biologist, it is hard to find problems that are simultaneously interesting, biologically relevant, and amenable to computational analysis using exclusively existing publicly available data sets. Fortunately, RNA's dual nature—being both experimentally tractable and mechanistically dynamic—means that there is a uniquely rich ecosystem of data sets spanning all aspects of its biology, providing ample opportunity to find curious problems with methodological and statistical challenges waiting to be solved.

If you were able to give one piece of advice to your younger self, what would that be?

There's no rush. The transition from undergraduate to graduate work was an exciting time for me. The shift in my experience of science—from passive observer consuming lectures and completing exams with prescribed answers, to active participant engaging with uncertainty and open discussion—was intoxicating. The rigidity of school had always felt stifling, but I had suddenly discovered a world where previously marginalized behaviors such as experimentation, dissent, and creativity were suddenly encouraged. I was fortunate that this transitional period coincided with my time in Richard Goldstein's group at UCL; the lab was generous in making me feel I had a seat at the table, even though I had no right to be. However, I got ahead of myself. I didn't fully appreciate how much of a beginner I was and would seek out questions and problems well beyond my capabilities at the time. In hindsight, I think I could have made things a lot easier for myself later if I had slowed down, and tried to focus on solving a small problem well before moving on to bigger things.

Are there specific individuals or groups who have influenced your philosophy or approach to science?

I want to give a shout-out to my PhD advisor, David Pollock, who demonstrated a certain kind of bravery and conviction, reflecting his optimistic faith in science as the pursuit of truth, rather than citation count. Sometimes, corners of the academic world can get caught up in dogma—even unconsciously—and overly fixated on superficial details while missing the key questions. It's an academic right of passage to feel frustrated at seeing papers gain favor more for aligning with prevailing views than for scientific merit, but I think the best (and maybe the only) response is to actively try and find a better approach. Challenging the consensus in a scientifically grounded manner requires significant courage and conviction, and David showed me what that could look like.

| Table of Contents
OPEN ACCESS ARTICLE

This Article

  1. Published in Advance December 18, 2025, doi: 10.1261/rna.080858.125 RNA 32: 328-340 © 2026 Pike and Schwartz.; Published by Cold Spring Harbor Laboratory Press for the RNA Society

Article Category

ORCID

  1. Profile for Pike, H. N.http://orcid.org/0009-0003-6786-5631
  2. Profile for Schwartz, S.http://orcid.org/0000-0002-3671-9709

Share

Navigate This Article

Current Issue

  1. July 2026, 32 (7)
  1. Current Issue
  1. Alert me to new issues of RNA