m6ADyn simulates subcellular transcript localization and m6A methylation level sensitivity to parameters

a) Schematic representation of the circular relationship between m6A methylation level, nuclear:cytoplasmic localization, and mRNA degradation. (top) m6A enhances cytoplasmic degradation, while cytoplasmic m6A facilitated degradation decreases m6A level. (bottom) Cytoplasmic degradation enhances relative nuclear abundance of a gene, while increased nuclear localization reduces susceptibility to m6A-mediated cytoplasmic degradation, thereby increasing stability.b) Scheme of the m6ADyn parameter (left). Equations for nuclear (Nuc) and cytoplasmic (Cyt) gene abundance based on three rates, production ⍺, export β and cytoplasmic degradation γ. Solutions, under steady-state conditions, for nuclear and cytoplasmic abundances of methylated and unmethylated transcripts, for overall methylation levels, and for relative nuclear:cytoplasmic levels are show. c) Heatmap indicating the simulated steady-state m6A level (left) and steady-state Nuc:Cyt localization (right) for a specific set of parameters. In both cases a gene was simulated using 1 for all the rates other than the rates indicated in columns or rows, for which we used increasing values (0.08 - 3.5).

Predictive modeling unravels interactions of steady-state m6A level, mRNA localization, and stability

a) m6A-GIs of a nuclear fraction and a cytoplasmic fraction, based on m6A-seq2 measurements in NIH3T3 cells in two replicates. b) Impact of nuclear:cytoplasmic localization on m6A levels. (left) m6ADyn predicted m6A levels, plotted as a function of m6ADyn predicted nuclear:cytoplasmic ratios, for simulated genes (n=2000). Parameters for alpha, beta, and gamma were sampled independently from a gamma distribution (rate = 1, shape = 1). (right) Experimentally measured m6A gene levels (m6A-GIs, y-axis) as a function of relative nuclear:cytoplasmic steady-state localization measured in mESC WT, divided into four equally sized bins (right). c) Distribution of difference in nuclear:cytoplasmic localization as a function of m6A levels for modeled data and experimental data, binned into 4 equally sized bins. (left) Simulations were conducted for 1000 genes. m6A-depleted conditions were modeled using γm6A = γA ,i.e. no m6A facilitated cytoplasmatic degradation, while in the control conditions, we used γm6A = 10 * γA. (right) The difference in nuclear:cytoplasmic localization measured in WT vs METTL3 KO mESCs, plotted as a function of m6A-GIs measured in mESCs WT. d) Impact of m6A levels and nuclear:cytoplasmic localization on mRNA half-life. (left) Heatmap showing simulated m6A-dependent half-lives, binned by simulated steady-state Nuc:Cyt ratios and m6A-gene level based on 1000 simulated genes using m6ADyn (rate sampling performed as in a). Genes were binned into 25 bins based on their m6ADyn m6A levels (X axis) and their m6ADyn nucleus:cytoplasmic levels (Y axis). The median simulated m6A-dependent half-life per gene is color-coded. (right) Corresponding values for experimentally-derived data, on the basis of half-lives, nuclear localization, and m6A data generated in WT and METTL3 KO mESCs. e) Associations between METTL3-dependent half-lives and m6A levels as a function of nuclear:cytoplasmic localization, binned into four equally sized groups. (top) Scatterplots of m6ADyn-derived half-lives as a function of m6ADyn-derived m6A levels (n=2000 gene simulations, rates sampled as in a). (bottom) METTL3-dependent half-lives (Mettl3 KO / WT HL) as a function of m6A-GIs based on measurements in mESCs (x-Axis).

Perturbations of RNA production or m6A methylation lead to predicted outcomes in localization and methylation level

a) Simulation of normalized m6A gene level after abolishing the production rate (⍺ = 0) for three simulated genes with different cytoplsamatic degradation rates. The three genes were simulated with identical rates (all rates equal 1, except γA and γm6A). The ony difference were three increasing degradation rates ( γA , with γm6A = 10 * γA ) and are indicated. b) Comparison of model-derived to experimentally-measured changes in m6A levels following transcriptional inhibition. (Left) m6A levels as predicted by m6ADyn on the basis of 2000 sampled genes for time points following perturbations in which production rate (α) is set to zero (rates were sampled as in Figure 2a). (Right) Experimentally (m6A-seq2) measured m6A-GIs in an A9 mouse fibroblast cell line along the indicated time points following ActD-mediated transcriptional inhibition. In both cases, the heatmaps are scaled across rows. c) Heatmap of row-scaled m6A-GIs of MCF7 cells treated with CPT (4h) and the control group DMSO (4h) d) Changes in m6A levels upon reducing production as a function of half-life. (left) m6ADyn simulated values (n=1000). WT rates were sampled as in Figure 2a. Reduced production was simulated by setting the production rate (α) to 10% of its initial level. (right) Distribution of m6A-GI changes following CPT treatment. Genes were binned into five equally sized bins based on their half-lives (right), human mRNA half-lives were obtained from (Agarwal and Kelley, 2022) e) Simulated changes over time of m6A levels (blue) and nuclear:cytoplasmic (red) localization upon elimination of m6A dependent degradation, by setting γm6A to equal γA. f) Row-scaled m6A-GIs based on m6A-seq2 measurements of NIH3T3 cells treated with siCtrl or with siRNAs targeting all three Ythdf1, -2, -3 (siDF1-3) g) Distribution of the fold change of Nuc:Cyt ratio of the siCtrl compared to siDF1-3, with the genes binned by the steady-state m6A level (m6A-GI in the siCtrl sample)

Changes in m6A gene modification under heat shock conditions are reflected by changes in mRNA metabolism

a) m6ADyn prediction of m6A gene level (top) and the genes relative localization Nuc:Cyt ratio (bottom) modeled by m6ADyn over time after perturbation of the export or production rate. The perturbations are the increase of the production rate by a factor of 10 or the decrease of the export rate to 10% b) Representative FISH poly-A staining signal and DAPI signal merged for timepoints of the cellular heat shock timecourse experiment (Scale bar = 20 μm) . c) Quantification of the poly-dt staining signal comparing mean intensity for nucleus / mean cytoplasmic signal. d) Scheme of a heat shock time course experiment (timepoints: ctrl/NHS, heat stress at 43 °C 0.75 h or 1.5 h and two timepoints after the recovery at 37 °C at 1 h (r1h) and at 4 h (r4h) e) PC1 of a PCA of the m6A-GIs based on the m6A-seq2 experiment. f) Heatmap depicting m6A-GIs, nuclear:cytoplasmic localization data, and mRNA stability for 1000 genes with the highest and lowest PC1 loading(plotted from high to low). Left, depiction of row-scaled m6A-GIs over the indicated time points and replicates. Middle, row-scaled Nuc:Cyt ratios. Right, available half-life measurements of NIH3T3 cells (column scaled) (Agarwal and Kelley, 2022). g) Fold-change in m6A level following an export block, as a function of difference in nuclear:cytoplasmic localization and METTL3-dependent half-lives. Values are plotted in the form of a heatmap, binned into 6 bins along each dimension (yielding a total of 36 bins). (left) m6ADyn simulations (n=1000, rates sampled as in Figure 2a) comparing steady-state m6A level and m6A level after reduction of export rates, with beta reduced to 10% of its WT levels. (right) Experimentally measured fold-changes of m6A level of heat shock 1.5h in comparison to steady-state control level, depicted as a function of changes in nuclear:localization and METTL3-dependent half-lives (half lives in WT compared to METTL3 KO cells). h) Boxplot showing the distribution of PC1 genes of genes binned by PRO-seq data (Mahat et al., 2016). PRO-seq data was divided into two bins: genes displaying a 10-fold increase in PRO-seq coverage upon heat shock (“transcriptionally induced”) and all remaining genes (“rest”). i) Heatmap depicting row-scaled normalized expression level (left) and m6A-GIs (right) based on three time points (triplicates) of Hsf1-KO cells and the five time points of the WT cells. (Top) Nine HSF1-target HSR-genes and (Bottom) 500 highest and lowest PC1 (based only on the WT-sampled) genes. j) Distribution of relative Nuc:Cyt localization of HSF1-targets (left) and for the other HSF1-independent genes (right). Based on Nuc:Cyt fractionation experiment based on WT and Hsf1 KO MEFs One-sided Mann-Whitney U test P values are indicated, testing for a higher Nuc:Cyt values in the wild-type fractionation sample.

smFISH library probes

List of the all the NGS datasets used for this manuscript

a) m6A-seq2 based m6A-SI (m6A sample index) for 2 replicates of cytoplasmic fraction and 2 replicates of the nuclear fraction in NIH3T3 cells. b) Distribution of m6A-GIs (scaled per gene across fractions and replicates) as a function of subcellular fractions, based on m6A-seq2 data acquired in MEF cells. c) Distributions of m6A level (m6A-GIs, y-axis) of genes as a function of relative nuclear:cytoplasmic steady-state localization, divided into four equally sized bins (x axis), on the basis of measurements in 3T3 cells d) Scatterplot of GLORI-derived gene indexes (Y axis) as a function of m6A-seq derived counterparts (X axis), both in MEFs. e) Distributions of GLORI based m6A-level (y-axis) of genes binned (equal-sized) according to their relative nuclear:cytoplasmic ratio f) Distributions of m6A level (m6A-GIs, y-axis) of genes binned (equal-sized) according to their relative nuclear:cytoplasmic steady-state localization measured in HEK293T g) Heatmap showing measured median fold-change half-life (HEK293T STM2457 / WT) in tiles. Genes are binned (equal-sized) along m6A-level (columns) and steady-state Nuc:Cyt localization (rows). h) Scatterplots of m6A-GI (x axis) and the fold-change half-life of STM2457/WT. Binned by the steady-state Nuc:Cyt measurements. All measurements are based on HEK293T mRNA.

Comparison of m6ADyn predictions based on different assumptions.1) m6A-dependent cytoplasmic degradation is 10-fold higher than non-methylated counterparts. 2) Model as in (1) but with nuclear degradation at a rate of δ. 3) Model as in (1) but assuming m6A facilitated export, in which the export rate of the methylated gene (βm6A) is by a factor (S) higher than the export rate of the the unmethylated part (βm6A = S * βA, with S = 10) 4) Model as in (1) but assuming decreased m6A dependent export rate (βm6A) compared to the export rate on the unmethylated transcript (S = 0.1). Following analyses are based on 2000 independently sampled rates from a gamma distribution (shape = 1, scale = 1) a) Fold change of steady-state nuclear to cytoplasmic m6A level predicted by the respective model (as in Figure 2a ) b) Distributions of m6A level of modeled genes binned according to their relative nuclear:cytoplasmic steady-state localization simulated with the respective model (as in Figure 2 b ) c) Scatterplot of modeled m6A level (x-axis) to the simulated half-life (y-axis), spearman correlation indicated

a) m6A-seq2 based m6A-SIs of the ActD timepoints b) m6A-SIs based on the CPT experiment m6A-seq2 measurements c) Boxplot of fold changes of m6A-GIs 6-hours following ActD treatment vs at 0h as a function of gene half-lifes (half-lifes taken from (Agarwal and Kelley, 2022) d) Knock-down efficiency comparing normalized expression level (TPM) of the siRNA targets YTHDF1,-2,-3 (DF1, DF2, DF3) e) m6A-SI based on m6A-seq2 measurements of NIH3T3 cells upon 48h treatment with siCtrl or with siRNA targeting Ythdf1, -2, -3 (siDF1-3)

a) Heatmap of m6A-GI pairwise Pearson’s correlation of the heat shock time course m6A-seq2 experiment b) SCARLET measurements of m6A at three sites after heat shock and recovery. The m6A-sites (mm9 genome) are Hsph1 (chr5.150420086.-), KLF6 (chr13.5867824.+), and Syf2 (chr4.134493325.+). c) Scatterplot of the m6A-GI heat shock timecourse PC1 (x-axis) against nuclear to cytoplasm relative localization of heat shock 1.5 h (hs1.5). Pearson’s R indicated d) Heatmap of the 1000 highest and lowest PC1 genes (high to low). Depiction of scaled m6A-GI (left) and expression (right) over the samples e) Gene Ontology analysis based on the top 1000 PC1 genes. The barplot shows the 10 most significant GO clusters and indicates the p-value.

a) Distributions of simulated m6A-level fold-change based on 1000 genes simulated, with independent sampling of rates with a gamma distribution under control conditions γm6A = S*γA with S =10 and test conditions in which γm6A = γA. The genes are binned in equal-sized bins (n=250) according to their relative nuclear:cytoplasmic ratio (left). Fold-change of m6A-level of triple knockdown siDF1-3 compared to siCtrl, with the genes binned by their relative nuclear:cytoplasmic ratio. b) Distribution of m6A-GI changes following CPT treatment for 4h, genes binned by measured steady-state Nuc:Cyt ratio (right) and fold change m6A gene level based on 1000 sampled genes after reducing production to 10% compared to steady-state m6A gene level and binned by steady-state localization of MCF7 cells (left)