Introduction

N6-methyladenosine (m6A) is the most abundant internal modification in mammalian mRNA. Initially discovered five decades ago (Desrosiers et al., 1974; Perry and Kelley, 1974), research into this modification was catalyzed with the advent of transcriptome-wide detection methodologies, namely m6A-seq/m6A-meRIP (Dominissini et al., 2012; Meyer et al., 2012). These methods rely on the immunoprecipitation of m6A-containing RNA fragments, followed by high-throughput sequencing and facilitate the extent of m6A methylation stochiometry on the m6A-site, -gene (m6A-GI) and sample level (m6A-SI) (Dierks et al., 2021). m6A has been implicated in most steps of mRNA fate and the disruption of the m6A-methylation machinery has critical manifestations in diverse systems (Batista et al., 2014; Cui et al., 2017; Ke et al., 2017; Li et al., 2017; Liu et al., 2015; Roundtree et al., 2017; Shi et al., 2017; Wang et al., 2015, 2014). Nonetheless, the most consistently observed molecular outcome of m6A is cytoplasmic mRNA degradation (Du et al., 2016; Lasman et al., 2020; Lee et al., 2020; Li et al., 2018; Uzonyi et al., 2023; Zaccara and Jaffrey, 2020), mediated via cytoplasmic ‘m6A readers’ of the YTH family (Du et al., 2016; Wang et al., 2014). In mouse embryonic stem cells (mESCs), m6A appears to be the predominant driver of mRNA stability, explaining roughly 30% of the variation in mRNA half-lives (Dierks et al., 2021; Uzonyi et al., 2023).

A fundamental question pertaining to m6A is the extent to which it is dynamically modulated within cells (i.e., across subcellular compartments) and across biological stimuli (such as stress responses). In both cases, the literature to date supports the notion that, at first approximation, the m6A methylome is fixed, and yet that m6A levels may be subjected to changes across both dimensions. Specifically, m6A profiles were found to be qualitatively similar across highly different cell types, developmental stages and tissues (Liu et al., 2020; Schwartz et al., 2014). Nonetheless, changes in m6A levels - often at only a minority of sites - have been reported across diverse systems, ranging from cellular responses to diverse stressors, infection, disease, and during development (Anders et al., 2018; Knuckles et al., 2017; Lichinchi et al., 2016; Liu et al., 2022; Su et al., 2018; Zhou et al., 2015). Similarly, in considering m6A levels across different compartments within cells, studies have consistently found m6A levels across chromatin-associated, nucleoplasmic and cytoplasmic mRNA fractions to be very similar (Ke et al., 2017; Roundtree et al., 2017). Nonetheless, slightly higher levels within the chromatin-associated fraction in comparison to the nucleoplasmic and cytoplasmic fractions (Ke et al., 2017) were observed in one study, contrasting with a more recent study finding reduced m6A levels in chromatin-associated RNA.

m6A deposition is catalyzed via a megadalton ‘m6A-writer’ complex, in which METTL3 acts as the catalytic component. m6A is to a considerable extent hard-coded, governed by the presence of the m6A consensus motif (‘DRACH motif’) and distance from splice sites (Garcia-Campos et al., 2019; He et al., 2023; Luo et al., 2023; Uzonyi et al., 2023; Yang et al., 2022). In light of a ‘hard-coded’ deposition baseline, it is interesting to dissect how dynamics in m6A level across compartments or biological stimuli can be achieved. The vast majority of the literature has focused on ‘active’ dynamics, wherein changes in m6A levels were proposed to be a consequence of changes in deposition or removal of m6A at specific sites (Hess et al., 2013; Knuckles et al., 2017; Su et al., 2018; Zhao et al., 2014; Zhou et al., 2015). Yet, in the context of other types of marks, ‘passive’ dynamics have been established, in which changes in the level of a modification are not caused by direct deposition or removal but instead result from shifts in metabolism affecting the molecules that carry these modifications. For example, 5mC DNA methylation during early development is primarily thought to be regulated passively via cell division, whereby each division results in a two-fold dilution of the DNA 5mC content (Mayer et al., 2000; Rougier et al., 1998). It is thus intriguing to speculate that m6A deposition might, too, be shaped by passive mechanisms. Nonetheless, a complicating factor in the context of m6A (with respect to 5mC or other instances of classic passive dynamics) is that the methylation mark itself shapes the metabolism of the transcript harboring it. Considering the well-characterized role of m6A in triggering mRNA degradation, m6A levels, therefore, have the potential of being governed by what we term ‘semi-passive’ dynamics, wherein m6A actively triggers the decay of transcripts harboring it, thereby passively altering m6A levels. Under such a model, m6A levels, mRNA decay and subcellular RNA localization are intricately coupled.

We previously established an approach, m6A-seq2, allowing to capture m6A level either for a specific site (m6A-site score) or a metric capturing the overall methylation for a whole gene body (m6A gene index’, m6A-GI) (Dierks et al., 2021). Inspired by several instances in which we observed gene-level changes in m6A and where such relationships were manifested differently across different subcellular compartments, we explored whether changes in m6A levels can be attributed to semi-passive dynamics. We establish a model capturing key steps in the life cycle of an mRNA. Despite the simplicity of the model and the fact that it does not allow any site-specific active reshaping of m6A levels, it predicts ample modulation of m6A levels, and captures relationships between mRNA localization, mRNA half-lives and m6A levels that could be confirmed based on experimental data. We propose that m6A dynamics, when present, are attributable in part to semi-passive regulation, wherein changes in mRNA metabolism impact the susceptibility of a transcript to m6A-dependent decay, thereby altering m6A levels.

Results

Previous studies have explored active mechanisms via which m6A levels could be modulated. We sought to examine whether changes in m6A could be passive consequences of changes in mRNA metabolism. The motivation for exploring this model was the realization that m6A levels both reflect and drive mRNA intracellular metabolism. Specifically, we considered a simple model wherein (1) m6A promotes mRNA decay, and (2) such decay occurs exclusively in the cytoplasm. This model intuitively links m6A levels and mRNA localization (nucleus/cytoplasm) (Figure 1a). On the one hand, under this model, changes in m6A cause changes in intracellular mRNA distribution: An increase in cytosolic m6A level of a given gene will trigger selective decay of cytoplasmic m6A harboring transcripts causing increased nuclear:cytoplasmic ratios of the gene. On the other hand, under this model, changes in mRNA localization will lead to changes in m6A levels: An increase in relative nuclear abundance of an mRNA leads to decreased cytoplasmic m6A-dependent degradation, thereby increasing overall m6A levels (Figure 1a). Thus, under these model assumptions, m6A levels, cytoplasmic decay, and subcellular localization are deeply intertwined, with changes in m6A driving and reflecting mRNA metabolism.

Figure 1.

To model the complex interplay between m6A levels, mRNA localization, and mRNA stability, we established ‘m6ADyn’, a differential-equation-based model. In m6ADyn, a transcript is produced in the nucleus at a rate of ⍺, gets exported from the nucleus into the cytoplasm at a rate of β, where it is degraded at the rate of γ. To capture the difference in the metabolism of m6A-containing and non-containing mRNA, we established two sets of parameters: ⍺^A, β^A and γ^A representing production, export and decay of nonmethylated transcripts, respectively, and ⍺^m6A , β^m6A and γ^m6A representing these rates for methylated counterparts. Throughout our simulations, unless indicated otherwise, we assumed (1) that the export rates of methylated and unmethylated transcripts were identical (β^A = β^m6A), and we will therefore refer to these rates as β, and (2) that the fate of methylated mRNA differed from that of its unmethylated counterparts exclusively in increased susceptibility to cytoplasmic decay, which was implemented by multiplying the unmethylated decay rate (γ^A) by a constant (γ^m6A = S * γ^A, with S > 1). Under steady state conditions (e.g. under normal growth conditions), the model equations can be solved analytically to determine mRNA and m6A levels in the nucleus and in the cytoplasm (Figure 1b).

m6ADyn further enables us to determine the full dynamical response of mRNA and m6A levels under non steady state conditions. For instance, m6ADyn allows us to explore how the modulation of parameters governing export rate or cytoplasmic m6A-dependent decay impacts m6A levels. Increasing either m6A-dependent cytoplasmic degradation (γ^m6A) or export (β) leads to reduced m6A levels - the former due to depletion of m6A harboring transcripts, and the latter due to increased levels of m6A-harboring mRNAs in the cytoplasm, rendering them susceptible to m6A-dependent decay (Figure 1c, left). In parallel, m6ADyn also allows evaluation of how different aspects of metabolism are impacted by m6A levels. For example, increasing the m6A deposition rate (⍺^m6A) results in increased nuclear localization of a transcript, due to the increased cytoplasmic decay to which m6A-containing transcripts are subjected (Figure 1c, right). Thus, m6ADyn allows us to determine how m6A levels drive and reflect mRNA metabolism.

We note that the assumptions made by m6ADyn are reductionist and possibly over-simplistic, given that m6A has been tied, among others, to production, transcription termination, nuclear decay, and export (Hsu et al., 2019; Ke et al., 2017; Knuckles et al., 2017; Roundtree et al., 2017; Slobodin et al., 2020; Yang et al., 2019), all of which are either not parameterized by m6ADyn or not assumed to be different between m6A harboring and non-harboring transcripts (of note: we do expand m6ADyn, below, to test a subset of these). Yet, despite its simplicity, the predictions of this framework are in many cases non-intuitive, and modeling them provides an opportunity to explore which aspects of m6A dynamics can or cannot be explained using only the above-described simple and well-supported assumptions.

Testing m6ADyn predictions under steady state conditions

We sought to assess the extent to which relationships predicted by m6ADyn under steady state conditions are supported by experimental data. We began with two predictions on how mRNA metabolism impacts m6A levels. A first prediction of m6ADyn concerns heterogeneity in m6A levels across sub-cellular fractions: m6ADyn predicts that nuclear m6A levels per gene would be invariably higher than the corresponding cytoplasmic levels, given the selective removal of m6A harboring transcripts from the cytoplasm (γ^m6A = S*γ^A, with S > 1). In m6ADyn, nuclear and cytoplasmic m6A can be analytically calculated for each gene in each compartment, defined as the fraction of transcripts harboring m6A divided by the overall number of transcripts in that compartment (Material and Methods Eq. 2-3). To confirm this experimentally, we fractionated NIH3T3 cells into nuclear and cytoplasmic fractions in duplicates and applied m6A-seq2 to all samples, following which we calculated m6A levels at the sample level and gene levels (m6A-SI & m6A-GI) across both compartments. Consistent with m6ADyn predictions, we found that m6A levels were higher in the nuclear fraction both at the sample and at the gene level, with the vast majority of genes (∼ 82%) exhibiting higher methylation levels in the nucleus (Supplementary Figure 1a & Figure 2a). These trends were also recapitulated in an independent experiment performed in primary mouse embryonic fibroblasts (MEFs) (Supplementary Figure 1b).

Figure 2.

A second prediction of m6ADyn concerns heterogeneity in m6A levels between genes: m6ADyn predicts that m6A levels across the entire cell are higher in genes that are relatively more abundant in the nucleus, given that they are least susceptible to m6A-dependent decay (Figure 2b, left). To assess whether this could be experimentally confirmed, we performed RNA-seq on nuclear and cytoplasmic fractions derived from WT and METTL3KO mESCs, which allowed the transcriptome-wide estimation of relative nuclear:cytoplasmic localization. Correlating the WT fractionation data against a previously acquired dataset capturing m6A levels per gene in mESCs (Dierks et al., 2021) revealed that indeed m6A levels were higher in more nuclear localized genes, consistent with m6ADyn predictions (Figure 2b right). We made similar observations also in NIH3T3 and in HEK293T cells for which we acquired fractionation data, confirming the robustness of these results (Supplementary Figure 1c, f). To control for potential technical covariates in antibody-based mapping, we additionally used previously published dataset of 40109 m6A sites, captured at single-nucleotide resolution, identified via GLORI, an antibody independent method relying on resistance of m6A sites to chemical deamination (Liu et al., 2022). We aggregated the signal stemming from individual m6A sites within a gene into a GLORI-gene index (GLORI-GI) (Materials & Methods). A comparison of GLORI-GIs with our m6A-seq-derived m6A-GIs showed a high correlation of R = 0.71 (Supplementary Figure 1d), supporting the robustness of m6A-GIs. Moreover, GLORI-GIs were higher in more nuclear genes (Supplementary Figure 1e), consistent with our observations on the basis of m6A-GIs.

We next sought to assess m6ADyn predictions relating to how changes in m6A levels affect mRNA localization and stability. We first focused on mRNA subcellular localization. m6ADyn predicts that loss of methylation (as in Mettl3 KO conditions) in genes that are highly methylated in wild-type conditions would lead to increased cytoplasmic presence (and hence reduced nuclear:cytoplasmic levels), whereas such an effect would be less pronounced in poorly methylated genes. Consistently with this model prediction, a comparison of localization in WT and m6A-depleted mESCs revealed that loss of m6A results in an increased propensity for highly methylated genes to become cytoplasmic (Figure 2c). Notably, relationships between m6A levels and localization have previously served as potential indications of m6A impacting nuclear export (Edens et al., 2019; Hsu et al., 2019; Kim et al., 2021; Lesbirel et al., 2018; Roundtree et al., 2017). The results demonstrate that even when assuming identical export rates between m6A-containing and non-containing transcripts, substantial differences in m6A gene level can be explained solely by differences in subcellular localization and finally through cytoplasmic decay (Figure 2b-c).

A final set of predictions of m6ADyn concerns m6A-dependent decay. m6ADyn predicts that (1) cytoplasmic genes will be more susceptible to increased m6A mediated decay, independent of their m6A levels, and (2) more methylated genes will undergo increased decay, independently of their relative localization (Figure 2d left). To explore these predictions, we relied - in addition to the above-described fractionation and m6A-seq2 measurements - also on previously published Actinomycin D (ActD) based mRNA stability measurements in WT and Mettl3 KO mESCs (Ke et al., 2017). For both experimental data, and model simulations, mRNA half-lives (HL) were acquired via linear modeling of the logarithmic decrease in mRNA levels over the timecourse following transcriptional inhibition (Material & Methods). Strikingly, the experimental data supported the dual, independent impact of m6A levels and localization on mRNA stability (Figure 2d, right). This trend could also be recapitulated in HEK293T, when using available m6A-seq, subcellular fractionation, and mRNA stability measurements following treatment HEK293T cells with either STM2457 (a small molecule inhibitor of METTL3) or mock (DMSO) (Supplementary Figure 1g) (Yankova et al., 2021). In our previous studies, we were able to attribute roughly 25% of the variability in mRNA half-lives to m6A levels (Dierks et al., 2021; Uzonyi et al., 2023). Here, however, analyzing this relationship through the lens of subcellular localization revealed that among cytoplasmic genes - most susceptible to m6A-dependent decay - the association between localization and m6A is substantially stronger, with 44% of the variability in METTL3-dependent decay being attributable to m6A levels. Conversely, among nuclear genes this relationship is substantially weaker (R² = 0.19). These trends are qualitatively predicted by m6ADyn (Figure 2e) and were further recapitulated using HEK293T measurements (Supplementary Figure 1h).

We next sought to assess whether alternative models could readily predict the positive correlation between m6A and nuclear localization and the negative correlations between m6A and mRNA stability. First, we assessed how nuclear decay might impact these associations by introducing nuclear decay as an additional rate, δ. We found that both associations were robust to this additional rate (Supplementary Figure 2a-c). Next, we sought to assess how m6A-facilitated export (β^m6A = 10 * β^A) might impact the predictions by the model. Given the selective depletion of m6A from the nucleus, due to the facilitated export of the methylated transcript from the nucleus, the associations between m6A and nuclear localization were reversed (Supplementary Figure 2a-c). This model is, therefore, inconsistent with measurements. Finally, we also assessed an inverse model wherein m6A is inhibitive of export (β^m6A = 0.1 * β^A). This model recovers the positive association between m6A and nuclear localization but gives rise to a positive association between m6A and decay - again contrasting with measurements (Supplementary Figure 2a-c). Thus, among the tested models, only ones in which m6A negatively impacts stability allow qualitatively recovering both the association with localization and with stability.

Testing m6ADyn predictions upon perturbations of mRNA metabolism

Next, we sought to assess m6ADyn predictions for how perturbations of mRNA metabolism impact m6A levels. m6ADyn predicts that decreasing production values (i.e. decreasing ⍺ rates) results in a decrease in m6A levels. Moreover, m6ADyn predicts that the rate at which m6A levels decrease is inversely correlated with mRNA stability, i.e. that m6A levels will decrease faster in more short-lived genes (Figure 3a). To test this prediction, we employed m6A-seq2 to monitor m6A levels along a 6h-time course following the treatment of mESC cells with ActD. Consistent with the m6ADyn predictions (Figure 3b, left), gene-level measurements of m6A displayed a subtle yet consistent decrease in m6A levels (Figure 3b, right). This drop was also observed at the sample level (Supplementary Figure 3a), consistent also with previously reported mass-spec-based measurements that showed a decrease in m6A/A ratio following ActD treatment (Roundtree et al., 2017). Similarly, treatment of MCF7 cells with Camptothecin (CPT), a topoisomerase inhibitor that stalls RNAPII elongation, thereby reducing production rates (Listerman et al., 2006), led to reduced m6A levels, in this case with much more dramatic effect sizes on both the sample and gene-levels (Supplementary Figure 3b, Figure 3c). As predicted by m6ADyn (Figure 3d, left), the decrease in m6A was stronger in genes with lower half-lives (Figure 3d, right). Consistently, following both ActD and CPT administration, we observed that the decrease in m6A levels was inversely proportional to the stability of the genes (Supplementary Figure 3c, Figure 3d).

Figure 3.

To further test predictions of m6ADyn, we next sought to perturb m6A-dependent decay to assess whether predictions of m6ADyn following perturbations of γ^m6A could be experimentally confirmed. m6A-Dyn predicts that elimination of m6A-dependent decay (e.g. γ^m6A = γ^A ) will result (1) in increased overall methylation levels (Figure 3e) and (2) in decreased nuclear/cytoplasm localization due to stabilization of cytoplasmic transcripts. Moreover, this impact on stabilization would be predicted to be stronger for the more methylated genes. Given that m6A-mediated decay is primarily thought to be mediated via the three cytoplasmic YTHDF proteins (YTHDF1, YTHDF2, YTHDF3) (Du et al., 2016; Lasman et al., 2020; Zaccara and Jaffrey, 2020), we next conducted a triple knockdown of the YTHDF proteins in NIH3T3 cells followed by whole-cell m6A-seq2 measurements and subcellular fractionation coupled with RNA-sequencing. While the knockdown substantially reduced RNA levels of YTHDF1 and YTHDF2, it did not reduce YTHDF3 levels(Supplementary Figure 3d). Nonetheless, there was an overall increase, albeit subtle, in m6A levels on the gene and sample levels (Supplementary Figure 3e, Figure 3f). Moreover, as anticipated by m6ADyn, we noted a reduction in the nuclear/cytoplasmic levels of more highly methylated genes following this perturbation (Figure 3g). These analyses thus support the relevance of m6ADyn for interpreting dynamic changes in m6A levels following alterations in mRNA metabolism. Nonetheless, the relatively subtle observed effect sizes hint at limitations of the model and/or at pleiotropic effects induced by the metabolic perturbations (see Discussion) .

Changes in m6A levels in the biological response to heat stress are secondary to changes in mRNA metabolism

Finally, we tested whether our model can serve as a framework for exploring how altered mRNA metabolism during a biological response leads to alterations in m6A levels. The cellular response to heat shock is well-documented to be associated with significant alterations in mRNA metabolism, characterized by a pronounced decrease in global export and an increase in the production of heat-stress response (HSR) genes (Mahat et al., 2016; Shalgi et al., 2014; Trinklein et al., 2004a, 2004b). Additionally, the cellular heat shock response has been implicated in differential m6A methylation of heat shock response genes (HSR-gene) at hundreds of m6A sites (Knuckles et al., 2017; Liu et al., 2022; Meyer et al., 2015; Zhou et al., 2015). Both the decrease in mRNA production rates and the reduction in export are predicted by m6ADyn to result in increasing m6A levels, thereby rendering heat shock an attractive model for testing to what extent observed changes in m6A levels could be directly explained by changes in mRNA metabolism (Figure 4a).

Figure 4.

We first sought to confirm previous findings pertaining to the heat-shock-induced export block. We applied a single-molecule FISH protocol using a poly(dT) probe directed against poly(A) mRNA in mouse embryonic fibroblasts (MEFs). We assessed global changes in mRNA localization in cells grown under normal conditions (non-heat shock, NHS) at 37°C, cells subjected to 1.5 hours of heat shock at 43°C (heat shock, HS) and cells in recovery at 37°C for 1 hour and 4 hours, following the stress. We observed a dramatic increase in the nuclear:cytoplasmic ratio of mRNA upon heat stress treatment (Figure 4b-c), consistent with previously documented nuclear retention of mRNA upon heat stress (Shalgi et al., 2014).

We next applied m6A-seq2 to a heat stress and recovery time course. The experiment was performed in primary MEF cells, in biological triplicates, using the same time points as in the FISH measurements with the addition of an early time point 45 minutes following heat shock (Figure 4d). We obtained m6A-GIs for ∼8100 genes expressed at sufficient levels across all samples. In general, m6A levels across genes were highly correlated between all time points, consistent with previous studies demonstrating that m6A is, to a considerable extent, a constitutive feature of the mRNA. Nonetheless, subtle quantitative differences were apparent (Supplementary Figure 4a). To capture the most dominant trends in m6A levels along this time course, we performed a PCA analysis on the m6A-GIs of all samples in the time course. The first principle component (‘PC1’) captured a continuous increase (or decrease) of m6A methylation along the heat stress time course (Figure 4e). We hence use PC1 below as an unbiased metric for ranking genes based on whether they undergo induced m6A levels over time (‘high PC1’) or reduced (‘low PC1’). While this manuscript focuses primarily on m6A levels per gene, we confirmed the change of m6A methylation during the heat shock time course with SCARLET for three chosen m6A sites in genes that showed a dynamic increase in m6A upon heat stress (high PC1 genes) (Supplementary Figure 4b).

We next sought to understand whether – consistently with m6ADyn predictions – the changes in m6A levels could be linked with differences in mRNA subcellular localization during heat stress. We monitored relative subcellular localization over the course of heat stress via fractionation of cells into nuclear and cytosolic fractions, followed by RNA sequencing. Remarkably, we observed that the genes undergoing increased methylation upon heat stress (‘high PC1 genes’) became dramatically more nuclear during heat stress (Figure 4f). Conversely, the genes with decreased methylation (‘low PC1 genes’) exhibited a reversed trend (Figure 4f). Overall, PC1 loading of genes and nuclear:cytoplasmic fractionation display a correlation of R = 0.4 (Supplementary Figure 4c). These results thus tie alterations in m6A levels to alterations of subcellular localization induced during heat stress and are consistent with the notion that disrupted export during heat stress could underlie the systematic increase in m6A levels observed in heat stress.

In our analysis of the dataset and attempting to understand why only certain genes exhibit induced m6A levels than others, we made two observations: First, high PC1 genes (displaying induced m6A levels during heat stress) tend to have substantially higher turn-over rates (shorter half-lives) in comparison to low PC1 genes. Second, we observed that high PC1 genes become more nuclear upon heat stress (Figure 4f).To assess whether differences in nuclear localization upon export block, combined with differences in m6a-dependent half-lives between genes were sufficient to account for observed differences in m6A levels under heat stress conditions, we used m6ADyn to predict changes in methylation levels following a global reduction of export (down to 10% of steady-state export rates). We observed that METTL3-dependent half-life and Nuc:Cyt ratio of genes following export block were conditionally independent of each other, with each of them individually (and both of them additively) negatively correlating with extent of m6A induction (Figure 4g, left). Remarkably, highly similar trends were observed when analyzing the experimental dataset in a similar manner, binning genes based on their METTL3-dependent stability and based on their relative Nuc:Cyt localization after 90 minutes heat-stress treatment (Figure 4g, right). This analysis thus highlights how the complex interaction between diverse aspects of RNA metabolism and m6A levels can be captured by m6ADyn.

A closer inspection of high PC1 genes (in which m6A is induced during the cellular stress response) revealed that these were enriched in ‘cellular stress response’ gene ontology, comprising classical heat stress response (HSR-) genes such as Hspa1a, Hspb1, Hsp1d, Hsp1e, Hsph1 and Dnajb1 (Supplementary Figure 4e). An analysis of Pro-seq data, which evaluates RNA Pol II occupancy as a proxy for transcriptional rate in MEFss during heat stress, confirmed that genes that showed the highest transcriptional induction are genes that have a high PC1 loading, i.e. show increased m6A level upon heat stress (Figure 4h). This was of interest to us, as a prediction of m6ADyn is that increased production rates will result in increased m6A levels (Figure 4a). We, therefore, sought to validate that the increased methylation observed in these genes was indeed a consequence of their transcriptional induction. We reasoned that perturbing Heat shock factor 1 (HSF1), the central transcription factor governing the heat stress response, could allow assessing whether the increased m6A levels in HSF1-transcriptionally induced targets are truly attributed to increased production or potentially due to some other mode of regulation. To explore this, we applied m6A-seq2 to a heat stress time course in HSF1-KO primary MEFs. We first confirmed that also within this dataset, HSF1-independent genes displayed similar m6A methylation dynamics as observed in the WT cells (Figure 4i, bottom). In stark contrast, while HSF1 target genes remained expressed at basal levels also within HSF1-KO cells, their m6A dynamics were completely lost (Figure 4i, top). Consistently, fractionation experiments in HSF1 KO cells revealed that HSF1 target genes were no longer enriched in the nucleus, in stark contrast to the HSF-independent counterparts. Furthermore, fractionation experiments of heat-stress treated HSF1 KO cells followed by RNA-seq confirmed that the HSF1-target genes do not experience an increase in nuclear enrichment during heat stress, unlike in the wild-type cells. Conversely, HSF1-independent genes demonstrate a similar increase in nuclear enrichment in both HSF1-KO and wild-type cells (Figure 4j). These analyses thus demonstrate that the m6A increase observed within HSF1-target genes is a consequence of the transcriptional induction, in line with m6ADyn predictions that induced production should be sufficient to lead to increased m6A levels.

To summarize, our results indicate that during the heat stress response, m6A accumulates within a large set of genes, primarily within targets that undergo increased nuclear localization. This is likely primarily due to a block in export, but - in a subset of genes - due to increased transcription. In both cases, the impact on nuclear:cytoplasmic levels and on m6A levels are captured well by the m6ADyn framework. These results thus highlight the utility of the m6ADyn framework in capturing and predicting m6A dynamics in the context of a complex cellular response.

Discussion

Until recently, the rules governing the selectivity of m6A deposition were shrouded in mystery. In the absence of an understanding of how m6A deposition is governed, we also lacked a framework for considering rules governing m6A dynamics. Consequently, diverse models were proposed for both aspects, typically relying on some form of ‘active’ recruitment of the methyltransferase complex and/or of m6A erasers to target sites (Anders et al., 2018; Huang et al., 2019; Knuckles et al., 2018, 2017; Lichinchi et al., 2016; Liu et al., 2022; Patil et al., 2016; Su et al., 2018; Wen et al., 2018; Zhou et al., 2015). Recently, it was uncovered that m6A deposition could, to a considerable extent, be attributed to simple exclusion (“100-nt rule”), wherein m6A is installed by default at all eligible consensus motifs, unless they are within ∼100 nt of a splice junction (Uzonyi et al., 2023) in which case their deposition is inhibited by the exon-junction complex (He et al., 2023; Luo et al., 2023; Uzonyi et al., 2023; Yang et al., 2022). Here, we propose that changes in m6A levels, across compartments or upon stimuli, may be similarly predictable via a simple model, wherein they may be attributable - at least in part - to changes in mRNA metabolism. Under this model, no changes in active deposition or removal of methylation from specific sites are needed in order to explain changes in m6A levels, which instead merely reflect consequences of m6A-driven mRNA metabolism.

Passive, rather than active, control over m6A levels is somewhat reminiscent of the passive control over DNA methylation in mammalian systems. DNA methylation, like RNA methylation, is controlled by ‘writer’ enzymes, and can be reversed via ‘erasers’, hence rendering the ‘active’ model a viable and attractive one. Nonetheless, considering the dramatic dynamics to which DNA methylation is subject in the germline, these are mediated primarily passively, by dilution of DNA methylation over the course of cellular divisions (Rougier et al., 1998). Our findings also connect with recent findings in neural cells, revealing that differences in relative mRNA localization within cell bodies vs neurites could, to a considerable extent, be attributed to the stability of the mRNA. While mRNA localization in neurons had been primarily thought to be governed actively, via specific RNA binding proteins involved in directing mRNAs to their target locations, this study thus proposed a passive model, wherein the inherent stability of an RNA plays a role in determining its relative localization (Loedige et al., 2023).

The tight - yet indirect - interconnectedness of m6A with mRNA metabolic rates (e.g. production and degradation) highlights the need to exercise caution in the interpretation of experiments along these dimensions. In previous studies, it was often hypothesized that changes observed in m6A levels are likely to be predominantly attributed to active recruitment of methylases or demethylases (Knuckles et al., 2017; Su et al., 2018; Wei et al., 2018; Zhao et al., 2014; Zhou et al., 2015). Our study points at the possibility that such changes may be passive in nature, and attributed to changes in mRNA metabolism, including differences in production, export or decay. Similarly, in previous studies involving perturbation of the m6A pathway were often found to result in diverse changes on the composition (e.g. selective expression of certain splice isoforms) and in the subcellular whereabouts of the transcriptome, based on which the possibility of a direct and causal connection with m6A was hypothesized (Mao et al., 2019; Wang et al., 2015; Wen et al., 2018). Our findings raise the possibility that such changes could, at least in part, also be indirect, and be mediated by the redistribution of mRNAs secondary to loss of cytoplasmic m6A-dependent decay. Nonetheless, m6ADyn should only be considered as a framework for testing parsimonious assumptions, and cannot rule out previously proposed mechanisms regarding the impact of m6A on mRNA metabolism or the mechanisms governing m6A dynamics. For instance, we tested whether experimental measurements would be consistent with m6A-facilitated export, and found that under this scenario experimental measurements are not consistent with predictions (Supplementary Figure 2a-c). While this may suggest that m6A does not impact export, contrasting a number of studies on the topic, it could also in principle be consistent with a more complex model wherein m6A both impacts export and mRNA decay (Edens et al., 2019; Hsu et al., 2019; Kim et al., 2021; Lesbirel et al., 2018; Roundtree et al., 2017). Under such complex scenarios predictions of m6ADyn become a quantitative outcome of the relative impact of m6A on different processing steps, no longer lending themselves to qualitative interpretation.

While the power of m6ADyn lies in its conceptual simplicity, some of its key caveats also lie therein. First, our exploration of m6ADyn is simulation-based: Predictions of m6ADyn are explored using randomly sampled alpha, beta and gamma rates. Being able to experimentally estimate these parameters per gene would allow obtaining quantitative predictions for how m6A levels in each gene would be impacted by changes in production, export and degradation rates. While systematic acquisition of rates capturing different steps in RNA metabolism has been challenging, this is an area marked by recent advances on the basis of metabolic labeling combined with subcellular fractionation (Müller et al., 2023; Smalec et al., 2022). Second, m6ADyn operates on the assumption of a binary methylation state for genes (methylated vs. unmethylated). This simplification overlooks the likely reality of discrete methylation states for all transcripts of a gene with multiple m6A sites, wherein the susceptibility of a gene to m6A-mediated decay may be a function of its overall m6A levels. Third, m6ADyn may be over-simplistic in many of its assumptions on mRNA metabolism and outcomes of m6A. mRNA metabolism encompasses more than only production, methylation, nuclear export and decay. Nuclear decay, for example, could also play a major role in some cases, although several previous studies have found its role to generally be negligible (Müller et al., 2023; Smalec et al., 2022). Alternatively, localization with specialized non-membranous subcellular compartments such as P or TIS granules could also impact the susceptibility of transcripts to m6A-mediated decay (Cougot et al., 2004; Sheth and Parker, 2003). M6ADyn also does not account for nuclear envelope breakdown during cell division, which would redistribute transcripts between nucleus and cytoplasm. M6ADyn may also be over-simplistic with respect to its assumptions regarding the roles played by m6A. While m6ADyn only considers the well-established role of m6A in controlling cytoplasmic stability, m6A has also been implicated in a range of additional steps on the mRNA life cycle including splicing, nuclear degradation and export. Thus, while the ability of m6ADyn to predict diverse changes in m6A between and across compartments is indicative of its utility, it should not be interpreted as evidence ruling out the existence or relevance of regulatory levels that are not captured as part of m6ADyn.

In line with the above, we wish to highlight several results in which m6ADyn predictions were not fully, or not at all, supported. First, ActD treatment (perturbation of production) and knockdown of the three YTHDF proteins (perturbation of m6A-dependent decay) both resulted in only relatively subtle alterations in m6A levels. Second, there were two predictions pertaining to the relationship between nuclear/cytoplasmic localization and m6A levels for which we were unable to obtain experimental support in the context of production and YTHDF perturbations. (i) Inhibition of production is predicted by m6ADyn to lead to a stronger decrease in m6A levels for cytoplasmic genes than for nuclear counterparts, which we do not observe experimentally (Supplementary Figure 5a). (ii) Inhibition of m6A-dependent decay is predicted by m6ADyn to lead to a stronger increase in m6A levels in more cytoplasmic genes, in comparison to more nuclear counterparts, which we also do not observe experimentally (Supplementary Figure 5b). Third, the quantitative discrepancy between the relatively mild decrease in m6A levels following ActD treatment versus the much more dramatic decrease following CPT treatment is not accounted for by m6ADyn. Finally, also the correlation between degradation rates and the decrease of m6A level was more pronounced following CPT treatment in comparison to ActD treatment (Figure 3d, Supplementary Figure 2d). These inconsistencies may partially be accounted for by experimental limitations, such as insufficient knockdown of the YTH proteins (or additional factors promoting m6A decay), or complex pleiotropic effects on mRNA metabolism being invoked secondary to the treatment of cells via different lethal inhibitors. Alternatively, these inconsistencies may also reflect oversimplified assumptions of m6ADyn regarding m6A metabolism and the mRNA life cycle.

To conclude, it is important to highlight that changes in m6A levels, both across subcellular compartments and across stimuli, are relatively rare, and typically also subtle. Our findings suggest that many of these changes, when present, may be accounted for by changes in mRNA metabolism, rather than by active control over m6A at specific sites.

Materials & Methods

Cell line & cell culture

NIH3T3 (mouse), primary MEFs (mouse), MCF7 (human) and HEK293T (human) cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Biological industries, 01-052-1A) supplemented with 10% fetal bovine serum (FBS).

Cytosol-Nucleus Fractionation Protocol

Initially, two million cells were collected in cold PBS and centrifuged at 300 g for five minutes. After centrifugation, the supernatant was discarded, and the cells were resuspended in 150 µl of buffer A (15 mM Tris-Cl pH 8 , NaCl at 15 mM, 60 mM KCl, 1 mM EDTA pH 8 , 0.5 mM EGTA pH 8) with an addition of RNase inhibitor at 10U/ml added fresh. The cell suspension was then mixed with 150 µl of 2x lysis buffer and incubated for 8 minutes on ice. The suspension was centrifuged at 400 g for 5 min, and 200 µl of the supernatant was carefully transferred to a new tube, identified as the cytoplasmic fraction. Remove and discard the supernatant sup from nuclear pellet, and resuspend the pellet in one ml RLN buffer (Tris•Cl pH8.0 50mM, 140mM NaCl, 1.5mM MgCl2, 0.5 % Np40, 10 mM EDTA pH8). Then the supernatant was removed followed by RNA isolation was performed using TRIzol (Cat 12183-555) according to the manufacturers instructions.

Cell Treatment Assays

For the stability experiments, cells were plated in 6cm and grown to 80-90% confluency. The cells were treated with 10 μg/ml Actinomycin D (ActD) for the respective timepoints, and immediately followed by whole cell RNA extraction. For the Actinomycind D timecourse after inhibitions of METTL3 with STM2457 (Yankova et al., 2021), HEK cells were treated with 5μM inhibitor for 6h before the treatment with Actinomycin D and and timepoints were taken subsequently. For CPT treatment, cells were treated with 6 μg per ml medium on 80-90% confluency grown cells. 4h timepoints were taken with according control samples of cells grown in a similar volume of only DMSO. For the heat shock time course experiments, cells were grown on 6 cm plates up to 80-90 % confluency and then incubated at 43°C for the heat shock samples (45 min or 90 min) or 90 min 43°C followed by a recovery incubation at 37°C (1h or 4h). Control samples grown at 37°C permanently. For si-RNA mediated knock-down, cells were plated a day prior and grown to 80-90% confluency. Transfections for si-mediated knockdown were performed according to the Dharmafect protocol. For the triple knockdown, amounts of siRNA and siControl were adapted accordingly.

RNA isolation was performed using TRIzol (Cat 12183-555) according to the manufacturer’s instructions. For library preparation, polyA selection was performed according to the manufacturer’s protocol using Dynabeads™ mRNA Purification Kit (Cat 61006).

single-molecule Fluorescence In Situ Hybridization (smFISH)

One day prior to the experiment, cells were plated on coverslips in 6 well plates. Cells were fixed using cold 4% Paraformaldehyde (PFA) for 10 minutes at room temperature, followed by a rinse with PBSx1 and permeabilization with 70% cold Ethanol for 2 hours (until ON) at 4°C. On the following day, samples were washed with SSCx2 (1x 150 mM sodium chloride and 15 mM trisodium citrate, adjusted to pH 7.0 with HCl) for 5 minutes at room temperature, and formamide wash buffers concentrations (20% formamide) were prepared. Hybridization buffers (Dextran Sulfate 10% v/v, Formamide 15 %, 2x SSC buffer concentration, BSA 0.02 % w/v, 2mM Vanadylribonucleside, 1mg/ml E.Coli tRNA, Nuclease free water) at these concentrations were mixed with the probes (in a dilution 1:100) and incubated overnight at 30°C. Finally, post-hybridization washes were conducted with preheated wash buffer (WB 20% formamide) for 30 minutes at 30°C, and samples were mounted using a freshly prepared GLOX buffer, with DAPI used at a 1:8000 dilution for 2 minutes at room temperature for nuclear counterstaining, preparing them for microscopy analysis.

smFISH Image analysis

Images were taken by Nikon in .nd2 stacked format. From the 11 z-stacks the central five were stacked into a single image. DAPI channel was used to determine the ROI for nuclear areas. Cytoplasmic areas were determined as a fixed ellipse around the nucleus with a vertical and horizontal radius twice the size of the nucleus. Afterwards, GFP channel was subjected to a Gaussian blur (sigma = 2) and subtracted from the signal channel to eliminate background stemming from autofluorescence. For measuring the polyDt signal we measured the signal intensity and compared the mean signal of the respective ROIs (nuclear to cytoplasmic). For the smFISH stainings targeted against individual genes, the signal was converted to a binary mask (‘otsu’), and particles in the respective ROIs were measured, and the signal was compared.

Site-specific cleavage and radioactive-labeling followed by ligation-assisted extraction and thin-layer chromatography - SCARLET

SCARLET based on (Liu et al., 2013) was conducted as in (Garcia-Campos et al., 2019) for three selected m6A sites in genes that exhibited increased m6A methylation during the cellular heat shock response. SCARLET was conducted on polyA-selected mRNA samples from primary MEFs.

Transcriptomic Analysis & Metrics

Each dataset was first trimmed with cutadapt (Martin, 2011). The subsequent genome paired-end alignment for human and mouse mRNA samples were aligned to the mm9/hg19 reference genome using STAR v.2.5.3a with default parameters (‘--alignIntronMax 1000000’ parameter) (Dobin et al., 2013). Uniquely aligned reads were extracted for further use. For gene expression estimation (TPM) or to calculate the expected read counts, rsem-calculate-expression (rsem/1.3.3, (Li and Dewey, 2011)) was used.

m6A Gene Index & m6A Sample Index

m6A-GI and m6A-SI were computed according to Dierks et al. 2021 with bam2ReadEnds.R script (Garcia-Campos et al., 2019) to extract transcriptome-wide paired-end reads coverage to infer the local m6A IP enrichment and normalize it to the input coverage. In the case of the heat shock timecourse libraries, TMM normalization of Input and IP fraction was performed prior to generation of the m6A-GIs, in order to reduce variability between replicates.

m6A GLORI gene level

To calculate the m6A gene level, we downloaded the GLORI sequencing data (GSE210563, Liu et al. 2022) based on mouse embryonic fibroblasts. We processed the fastq files according to (Liu et al., 2022). Alignment was performed with HISAT-3N (Zhang et al., 2021) under default setting (with --base-change A,G). Bam files were processed into an .rds file (that can be processed in R Studio) with txtools (Garcia-Campos and Schwartz, 2024) using an mm9 reference genome and a canonical UCSC mm9 gene annotation. GLORI scores were calculated for annotated adenosines that lay centered in a DRACH motif environment. Using the ratio of adenosines (“A”, that were not converted to a “G”) to the coverage of the specific site, filtering for sites that have a coverage of at least 20 and a GLORI score higher than 0.1. To calculate the GLORI m6A gene index (GLORI-GI) we used the sum of GLORI scores derived from all eligible DRACH motifs and normalized it by the number of adenosines of the gene.

Nuclear:Cytoplasmic relative localization

We used RSEM (Li and Dewey, 2011) to determine expected counts of the nuclear and the cytoplasmic fractions. Cross-sample normalization for each of the two subcellular fractions was conducted using the TMM function (‘NOISeq’ package in R) (Tarazona et al., 2011), following which log fold-changes between normalized nuclear and cytoplasmic counts were calculated.

mRNA half-life estimation

mRNA decay-rate estimation was performed by linear modeling of the log-transformed normalized expression of the Actinomycin D time course samples against the respective time points. Decay rates that were positive or ones emerging from non-significant fits of the linear model (p>0.05) were removed.

‘m6ADyn’ model calculation

m6ADyn is a differential equation-based model for modeling mRNA metabolism In m6ADyn, transcripts are generated in the nucleus with the production rate α, exported from the nucleus to the cytoplasm with the export rate β, and subject to cytoplasmic decay at a rate of γ. One set of rates is defined for non-m6A harboring transcripts (⍺^A, β^A and γ^A), and another for m6A-harboring counterparts (⍺^m6A, β^m6A and γ^m6A). Unless specified otherwise, in our implementations and simulations γ^m6A >γ^A and β^A=β^m6A. The following two differential equations define this model:

Whereby Nuc and Cyt refer to the levels of a gene in the nucleus and cytoplasm, respectively.

Nuclear and cytoplasmic levels of a gene can be derived under steady-state conditions, where no change occurs over time .

Availability of nuclear and cytoplasmic abundances of m6A containing and lacking counterparts, allows to derive the following values under steady state conditions:

For deriving nuclear and cytoplasmic abundances outside of steady-state, we used the R-package ‘MOSAIC’ which simulates changes over time (Pruim et al., 2017). This package was also used to simulate stability rates by fitting a linear model to m6ADyn-derived log transformed mRNA levels along a timecourse following inhibition of production (⍺^A=⍺^m6A=0), mimicking an experiment relying on transcriptional inhibition.

For the simulation of genes, the necessary rates were sampled independently (⍺^A, ⍺^m6A, β, γ). We sampled rates from a gamma distribution, which, by default, has a rate parameter of 1 and a shape parameter of 1, unless specified otherwise. Unless specified, export rates (β) were equal for methylated and unmethylated transcripts (β^A=β^m6A ) and degradation rates for the methylated component was higher by a factor (S) of ten (γ^m6A = 10 * γ^A ).

m6ADyn with nuclear degradation

For this case we also include nuclear degradation (δ) as a rate. Resulting in minor changes of the equations.

Solving the equations for steady-state conditions

Data availability

All the m6A-seq2 datasets that were generated in the course of this paper are stored in GSE###

Code availability

The used in-house python script for demultiplexing as well as the script (bam2ReadEnds.R) used for paired-end transcript coverage are stored in https://github.com/SchwartzLab

Table 1.

Table 2.

Acknowledgements

SS is funded by the Israel Science Foundation (grant no. 913/21), the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant no. 101000970).

R.S.S. is supported by the Israel Science Foundation (grant number 395-21), and the European Research Council (ERC grant number 101043300). R.S.S. is the incumbent of the Robert and Yadelle Sklare Professorial Chair in Biochemistry