{'fraction_non_alphanumeric': 0.033988106535079395, 'fraction_numerical': 0.024957297399886127, 'mean_word_length': 4.848580149847377, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
mRNA N6 ‐methyladenosine is critical for cold tolerance in Arabidopsis

SUMMARY Plants respond to low temperatures by altering the mRNA abundance of thousands of genes contributing to numerous physiological and metabolic processes that allow them to adapt. At the post‐transcriptional level, these cold stress‐responsive transcripts undergo alternative splicing, microRNA‐mediated regulation and alternative polyadenylation, amongst others. Recently, m6A, m5C and other mRNA modifications that can affect the regulation and stability of RNA were discovered, thus revealing another layer of post‐transcriptional regulation that plays an important role in modulating gene expression. The importance of m6A in plant growth and development has been appreciated, although its significance under stress conditions is still underexplored. To assess the role of m6A modifications during cold stress responses, methylated RNA immunoprecipitation sequencing was performed in Arabidopsis seedlings esposed to low temperature stress (4°C) for 24 h. This transcriptome‐wide m6A analysis revealed large‐scale shifts in this modification in response to low temperature stress. Because m6A is known to affect transcript stability/degradation and translation, we investigated these possibilities. Interestingly, we found that cold‐enriched m6A‐containing transcripts demonstrated the largest increases in transcript abundance coupled with increased ribosome occupancy under cold stress. The significance of the m6A epitranscriptome on plant cold tolerance was further assessed using the mta mutant in which the major m6A methyltransferase gene was mutated. Compared to the wild‐type, along with the differences in CBFs and COR gene expression levels, the mta mutant exhibited hypersensitivity to cold treatment as determined by primary root growth, biomass, and reactive oxygen species accumulation. Furthermore, and most importantly, both non‐acclimated and cold‐acclimated mta mutant demonstrated hypersensitivity to freezing tolerance. Taken together, these findings suggest a critical role for the epitranscriptome in cold tolerance of Arabidopsis.


INTRODUCTION
The first RNA modification was identified in yeast in 1957 and, subsequently, over 160 distinct types of RNA modifications occurring predominantly on non-coding RNAs, such as transfer RNA, ribosomal RNA (rRNA), and small nuclear RNA, have been reported (Arribas-Hernandez & Brodersen, 2020;Boccaletto et al., 2018;Fray & Simpson, 2015;Hu et al., 2019;Liu & Pan, 2016;Vandivier & Gregory, 2018). Recent advances in detection methods have resulted in the identification of these modifications on very low abundant RNAs such as mRNAs, long noncoding RNAs, and microRNAs (Dominissini et al., 2012).
Among the mRNA modifications, N 6 -methyladenosine (m 6 A) is currently the most prevalent internal modification in eukaryotic mRNAs (Arribas-Hernandez & Brodersen, 2020;Hu et al., 2019;Kramer et al., 2018;Liu & Pan, 2016;Vandivier & Gregory, 2018). This modification on mRNAs plays important roles in affecting stability, degradation, and translation, amongst others. The m 6 A is added to an RNA by a large, multi-component protein 'writer complex' containing two methyltransferases [METHYL-TRANSFERASE A (MTA) and METHYLTRANSFERASE B (MTB)] and is removed by 'erasers' known as demethylases (ALKBH9 and ALKBH10) in plants ( Arribas-Hernandez & Brodersen, 2020;Duan et al., 2017). The m 6 A modification is highly dynamic as well as reversible and is known to be enriched in the sequence context of 'DRACH' (where D = A/G/U, R = A/G, H = A/C/U) of the mRNA) (Luo et al., 2014;Meyer et al., 2012;Wan et al., 2015). Although this canonical motif appears to be conserved in mammals, recent studies suggest the presence of other plant specific m 6 A motifs (Wei et al., 2018). The current consensus suggests that the presence of the motif alone is not sufficient to understand the dynamic transcript methylation patterns during different developmental stages and stress conditions, and thus is likely context-dependent.
In mammals, m 6 A has been implicated in various developmental and biological processes such as embryonic development and cell fate determination (Batista et al., 2014), circadian rhythms (Zhong et al., 2018), and cancer cell proliferation (Cai et al., 2019). In plants, the functions of m 6 A have largely been determined from the analysis of mutants that are defective in writer, reader and eraser complexes. For example, the disruption of Arabidopsis m 6 A writer components (MTA, MTB, FIP37, and VIR) results in developmental defects including embryo lethality (Zhong et al., 2008), reduced apical dominance (Bodi et al., 2012), over proliferation of shoot meristems (Shen et al., 2016), aberrant vascular formation in the root (Ruzicka et al., 2017), and leaf morphogenesis, as well as rate of leaf formation . In rice, m 6 A has been shown to be essential for microspore development based on the observation that Osfip (a component of writer complex) mutant plants demonstrate early degeneration of microspores (Zhang et al., 2019). m 6 A has also been shown to play roles in fruit ripening in tomato . Furthermore, m 6 A methylome studies conducted in various tissues of Arabidopsis and rice revealed tissue-specific changes in m 6 A levels, which in turn affects gene expression Wan et al., 2015). These studies illustrate the importance of m 6 A in plant growth and development but their roles in plant adaptation to stresses are only emerging. Only recently has the potential involvement of m 6 A in salt stress responses been reported in Arabidopsis Hu et al., 2021). Likewise, vir mutant plants demonstrated salt hypersensitivity (Hu et al., 2021). On the other hand, MTA overproducing Poplar plants in which the overall m 6 A levels were elevated showed enhanced drought tolerance (Lu et al., 2020).
Low temperatures suppresses leaf growth and expansion, overall plant growth and development, and flowering (Chinnusamy et al., 2010;Ritonga & Chen, 2020). Physiologically, cold stress inhibits the Calvin-Benson cycle, generates reactive oxygen species (ROS) and decreases membrane fluidity, thus affecting membrane functions. The low temperature-responsive gene expression profiles in Arabidopsis revealed extensive reprogramming of gene expression leading to the activation of COR genes (COR15A, COR78/RD29a, COR15B, KIN1, KIN2, and COR413-PMI) (Ding et al., 2019;Kidokoro et al., 2017;Ritonga & Chen, 2020;Thomashow, 2010;Yamaguchi-Shinozaki & Shinozaki, 2006). The accumulation of COR gene products is critical for the acquisition of cold acclimation and subsequent freezing tolerance in plants. Both CBF-dependent and CBF-independent pathways contribute to activating COR genes (Shi et al., 2017. The low temperature-responsive CBF-independent pathway includes various transcription factors, such as HSFC1, ZAT12, ZF, ZAT10, RAV1, CZF1, and HY5, all of which have been shown to activate COR gene expression . The CBF-dependent pathway genes (CBF1, CBF2, and CBF3) are transiently and rapidly induced during cold treatment to promote the subsequent activation of specific COR genes (Shi et al., 2017Thomashow, 2010). Interestingly, the regulation of CBFs themselves is quite complex, with positive regulation mediated by upstream factors such as ICE1, CAMTA3, BZR1/BES1, CESTA (CES), and circadian clock-associated 1/late elongated hypocotyl (CCA1/LHY), and negative regulation being mediated by factors such as MYB15, PIFs, EIN3, and SOC1 factors (Dong et al., 2011;Li et al., 2017a;. Thus, a very complex regulatory network influences the expression of CBFs during the cold stress response in plants. In addition to the transcriptional changes, post-transcriptional regulation such as microRNA-mediated regulation, alternative splicing, and alternate polyadenylation have been shown to play critical roles during cold stress responses (Calixto et al., 2018;Song et al., 2016;Sunkar, 2010). With the discovery of RNA modifications that could influence mRNA splicing, stability, degradation, and translation, this new posttranscriptional regulatory mechanism could also play important roles in plant adaptation to low temperature stress. To investigate the importance of m 6 A mRNA modifications on cold stress responses, we generated methylated RNA immunoprecipitation (MeRIP)-sequencing (seq), mRNA-seq and polysomal RNA-seq. The significance of m 6 A RNA methylation on cold tolerance was further assessed using the mta mutant, which demonstrated hypersensitivity as assesed by ROS accumulation and electrolyte leakage in response to freezing stress. Thus, our findings have revealed an important link between the m 6 A epitranscriptome modification and cold tolerance in plants.

RESULTS
MeRIP-seq identifies transcripts that are m 6 A-enriched or m 6 A-depleted under cold stress To understand the dynamics of m 6 A methylation in the overall transcriptome under low temperature stress, MeRIP-seq was performed on polyA + selected RNA in 3week-old Col-0 seedlings that were either cold treated or untreated (see Experimental Procedures) for 24 h. Using input polyA + selected RNA-seq as background, the peakcalling algorithm MACS2 was employed to identify m 6 A peaks across three biological replicates of both control and cold treated samples. The identified peaks showed between approximately 76% and 86% overlap among the replicates (Figure 1a) from which a list of unique highconfidence m 6 A peak for each condition was generated. In total, we identified 17 772 and 13 261 high-confidence peak regions from the control and cold-treated samples, respectively. The cold stress-responsive m 6 A peaks were determined based on the peak abundances in the control and cold-treated samples. This analysis identified approximately 1318 and 5829 m 6 A peaks that were designated as either cold-enriched or control-enriched (cold-depleted) m 6 A peaks, respectively, and the remaining 11 943 m 6 A peaks shared between the two conditions, implying that their methylation profiles did not change upon cold stress ( Figure 1a, Table S1). Cold-enriched m 6 A containing transcripts refer to transcripts that contain peaks that were enriched in all three replicates of cold-treated samples but present in two or less replicate samples of the control (not in all three replicates). The cold-enriched (1318 m 6 A peaks) and cold-depleted (5829 m 6 A peaks) are represented by 1198 cold-enriched and 4545 control-enriched (colddepleted) m 6 A containing transcripts, respectively. Among these, 233 transcripts contain both a control-enriched and a cold-enriched m 6 A peak, although in different positions of the transcript, implying that, under cold stress, they are depleted for one peak but are enriched for another peak.
To validate the MeRIP-seq profiles, several cold-enriched m 6 A containing transcripts (AT2G46690, AT4G19530, AT4G25990, AT5G17300, and AT5G62360) were analyzed using m 6 A-immunoprecipitation quantitative polymerase Figure 1. Arabidopsis m 6 A shows a 3 0 UTR positional bias, although the bias is stronger for cold unresponsive peaks. (a) High-confidence m 6 A peaks obtained by overlapping biological replicates of cold and control treated samples revealed 5829 control-enriched (cold-depleted) peaks, 1318 cold-enriched m 6 A peaks, and 11 943 common (unaltered) peaks. (b) Distribution of the meRIP-seq reads along binned transcripts. (c) RNA motifs identified by HOMER to be enriched in the indicated m 6 A peak regions. (d) Percentage distribution of cold-enriched, control-enriched, and shared m 6 A peaks that overlap with the 5 0 UTR, start codon, CDS, stop codon, and 3 0 UTR. (e) Nucleotide level probability metric for methylation in the start and stop codon boundary regions for cold-enriched, controlenriched, and shared m 6 A peak containing transcripts. Each group is normalized to its maximim value in the region. chain reaction (qPCR) ( Table S3). As expected, these transcripts were revealed to be more enriched in m 6 A pulldown fractions after cold treatment compared to untreated controls ( Figure S1). Furthermore, we compared the identified m 6 A peaks in the present study with two previously published studies that reported m 6 A peaks from the adult leaves and cotyledons Shen et al., 2016). Unsurprisingly, the overlap between the m 6 A peaks in the present study and the previous studies was more than 88% ( Figure S2).
To identify the localization of m 6 A peaks on the different positional segments of mRNA [5 0 UTR, coding sequence (CDS) and 3 0 untranslated region (UTR)], a relative coverage of the MeRIP-seq library reads was plotted across transcript UTR and CDS positions. This revealed that the m 6 A pulldown reads were highly enriched in the 3 0 UTR and near the stop codon of Arabidopsis transcripts, a typical feature of m 6 A containing transcripts in plants and animals ( Figure 1b). In comparison, the metaplot for background mRNA library shows that many of these reads tend to localize in the CDS region ( Figure S3). Overall, these results are consistent with m 6 A localization identified in previous reports from multiple eukaryotes including plants Kramer et al., 2020;Luo et al., 2014;Luo et al., 2020Luo et al., , 2021Zhou et al., 2015).
To identify motifs that are enriched for m 6 A peaks, we used the motif calling algorithm HOMER. We observed an enrichment of UGUA as the top motif in both cold and control samples (Fig. 1C). The UGUA motif has previously been identified as an m 6 A consensus motif in plants, observed in several plant species including Arabidopsis (Hu et al., 2021;Wei et al., 2018), rice , and maize (Luo et al., 2020). Other m 6 A motifs such as the AAGM and the GAD motifs were also found with high probability of a canonical 'A' as a site for m 6 A modification. We performed principal component analysis and generated heatmaps showing the clustering of meRIP-seq libraries and mRNA-seq libraries, which show distinct clusters of sequenced libraries ( Figure S4), verifying the quality and reproducibility of our various RNA-seq datasets.
m 6 A peaks that are enriched or depleted under stress are less strictly localized in the 3 0 UTR and have a relatively higher probability of occuring in CDS and 5 0 UTR regions Previous studies have suggested that the localization of m 6 A peaks on various mRNA segments are dynamic and can shift towards the 5 0 UTR region under stress conditions such as during heat stress in human cell lines (Zhou et al., 2015). To identify any such shifts in cold-stressed samples, the unique m 6 A peak regions identified in our cold-enriched, control-enriched, and shared m 6 A containing transcripts were mapped onto the annotated 5 0 UTR, start codon, CDS, stop codon, and 3 0 UTR (TAIR10 database). This analysis has revealed a small but noticeable increase in the percent of peaks that overlapped with the 5 0 UTR and start codon regions for the cold-enriched peaks compared to shared or control-enriched peaks (Figure 1d). Overall, compared to the 'shared' transcripts that do not show a change in m 6 A enrichment under cold stress, both the cold-enriched and cold-depleted m 6 A peak containing transcripts show a higher percentage of peaks in 5 0 UTR, start codon, and CDS ( Figure 1d). These observations imply that a larger percentage of the stress responsive peaks is localized in these upstream regions. However, MeRIP-seq cannot map the m 6 A peak at nucleotide level resolution, thus giving a range of nucleotide positions as a peak region. As such, the peaks that span two or three regions may be represented in all of those categories (5 0 UTR, start codon, and CDS).
Despite these limitations of the MeRIP-seq technique, we attempted to obtain a closer look at the segmental information of m 6 A on the transcripts by calculating the individual nucleotide probability of methylation within the three groups of transcripts without binning the transcripts (see Experimental Procedures). Probability maps around the stop codon regions reveal that the controlenriched (cold-depleted) m 6 A peaks have a higher relative probability of being localized in CDS regions upstream of stop codon compared to other two groups. Near the start codon, cold-enriched transcripts show a higher probability of having an m 6 A in the nucleotide positions surrounding the start codon compared to control-enriched or shared transcripts (Figure 1e). This observation along with the percentage localization analysis suggests that, although the shared m 6 A peaks that are insensitive to cold stress tend to localize in 3 0 UTRs, the more dynamic m 6 A peaks that are either enriched or depleted under cold have a relatively higher probability of being found in regions upstream of the 3 0 UTR of the methylated transcript. Notably, the position of maximum probability of these groups, particularly for cold-enriched m 6 A containing transcripts which show probability 'humps' immediately up and downstream of the start codon. Relatedly, a total of 101 transcripts were found to have a coldenriched m 6 A peak within these 'humps' located in the 100-nucleotide region spanning the start codon (Figure 1e). Although no functional group was enriched as assessed by Gene Ontology (GO) (http://geneontology. org) analysis, the overall mRNA fold change (abundances) for these 101 transcripts was slightly higher compared to the remaining transcripts containing m 6 A in other locations, as well as 5 0 UTR/start codon region m 6 A transcripts that were control-enriched or shared (Figure S5b). These observations suggests a role for m 6 A in the 5 0 UTR region in transcript abundance/stability, although additional studies are needed to validate this hypothesis.

Cold stress induces large-scale transcriptome changes
Cold stress induces large-scale physiological and biochemical changes in plants driven by altering expression of a large number of genes. To assess the change in mRNA levels after 24 h of cold stress, we performed mRNA-seq on Col-0 plants. Using DESEQ2, we observed major changes in mRNA abundance patterns resulting in significant increases (n = 2654) and decreases (n = 2440) in thousands of transcripts across the Arabidopsis transcriptome in response to low temperature treatment (Figure 2a). To assess the functionality of the proteins encoded by these transcripts, GO enrichment analysis was performed using the tool DAVID (Huang et al., 2009). As expected, the transcripts that are significantly (adjusted P < 0.05) upregulated during the cold stress response were enriched for biological terms including 'cold stress response' and 'cold acclimation', as well as 'rRNA processing and ribosome biogenesis' (Figure 2b). Importantly, among the upregulated transcripts, we found strong upregulation (a more than five-fold increase) for the mRNAs represented by well-known COR (COR4, COR15A, COR15B, COR27, and COR47), KIN (KIN1 and KIN2/COR6), and LTI mRNAs (LTI30 and LTI78/COR76) ( Figure S6a) that are direct downstream targets of CBF transcription factors (Shi et al., 2017. Several of these genes were validated using qPCR assays ( Figure S6b). On the other hand, general response related terms and single organism processes were significantly reduced at the mRNA level in response to cold treatment. Overall, our RNA-seq analysis was in agreement with the previously reported cold stressaltered transcriptome analysis in Arabiopsis (Fowler & Thomashow, 2002;Seki et al., 2001;Zhao et al., 2016).
Transcripts containing cold-enriched m 6 A peaks show increased abundance  Kramer et al., 2020). To uncover any correlation between m 6 A enrichment and transcript abundance in the present study, the overall change in transcript abundances of the coldenriched, control-enriched, and shared m 6 A peak containing transcripts were measured using mRNA-seq data. We plotted the overall log-transformed fold change for the abundance of each of these transcripts. This assessment revealed that the cold-enriched m 6 A containing transcripts showed a significantly higher increase in transcript abundance under stress compared to the transcripts that contain control-enriched (cold-depleted) and shared m 6 A peaks (P < 0.05 in all comparisons; Wilcoxon rank sum test) (Figure 2c). These observations suggest that transcripts that contain cold-enriched m 6 A peaks have increased abundance. As expected, these cold-enriched m 6 A peak cotaning transcripts are also belong to some of the well-known cold stress related transcripts with the GO terms such as 'response to stress, response to abiotic stimulus, response to cold and response to temperature stimulus' ( Figure S7).

Transcripts containing cold-enriched m 6 A peaks show increased polysomal association
Previous studies have shown that m 6 A can act as a mark that promotes translation when recognized by members of the YTHDF2 reader protein family in mammalian systems Wang et al., 2015). These studies have suggested an m 6 A reader-dependent mechanism forming a link between increased ribosome loading of methylated transcripts and increased elongation. To assess whether m 6 A affects the translational properties (promoting or inhibiting/decreasing their translation) of methylated transcripts during the cold stress response in Arabidopsis, we profiled mRNAs associated with the polysomes from control and cold-stressed Arabidopsis seedlings as described previosuly (Li et al., 2017b. The ribosome occupancy metric was calculated for each transcript where polysome profiling read count was normalized against its mRNA-seq read count to overcome potential bias as a result of high mRNA abundance affecting polysome loading (see Experimental Procedures). This allowed us to compute differences in polysome binding between cold and control treatments, as well as to evaluate whether transcripts are more or less likely to be translated under cold stress. As a control, we compared the global change in ribosome occupancy values between cold and control conditions and found no statistically significant changes ( Figure S8). However, the average ribosome occupancy value of transcripts containing cold-enriched m 6 A peaks demonstrated a significantly more positive change under cold compared to transcripts containing cold-depleted (control-enriched) m 6 A peaks (P = 8.1 × 10−7 ; Wilcoxon rank sum test) or those that remain unchanged in methylation status (P = 0.038; Wilcoxon rank sum test) ( Figure 3a). To further narrow down and identify which transcripts within the coldenriched m 6 A containing transcripts were driving this average positive occupancy change, we divided this group into those that show an actual positive or negative ribosome occupancy change and performed a GO analysis. The results showed separate clusters of significant GO terms where transcripts associated with rRNA processing, karrikin response, and response to cadmium show increased ribosomal occupancy, whereas transcripts associated with photosynthesis, leaf development, and chloroplast organization show decreased ribosome occupancy ( Figure 3b). Transcripts directly related to cold stress response, although seen in both groups (positive or negative ribosome occupancy), showed a higher significant enrichment in the category that has increased ribosome occupancy values. These findings suggest an underlying mechanism governing the polysome association of cold-induced m 6 A containing transcripts belonging to specific functional categories.
Major cold responsive (COR) genes show increased profiles in MeRIP-seq, mRNA-seq, and polysome profiling To consolidate the data from the various sequencing experiments, IGV browser views (Thorvaldsdóttir et al., 2013) were generated for the major COR genes (CBF2, COR15A, COR15B, COR27, COR78, and COR413) using their normalized read counts from our mRNA-seq, MeRIP-seq, and polysome profiling (Figure 3c). The browser views revealed that when exposed to 24 h of cold treatment, the overall normalized read counts for these transcripts are significantly increased in the mRNA-seq, as well as even more significantly increased in the MeRIP-seq and polysome profiling experiments, relative to the untreated controls. This demonstrates that the increase in m 6 A methylation goes hand in hand with the increased mRNA abundance and ribosome occupancy for the COR transcripts under cold stress. To inspect whether this trend is biologically significant or a result of bias introduced by higher transcript level resulting in higher polysome association, as well as m 6 A detection, genome browser views were generated for several random transcripts. We found several examples where enrichment of m 6 A, mRNA, and polysome-associated RNA did not go in parallel upon cold treatment ( Figure S9). Thus, our findings for the COR genes are likely of bona fide biological significance, as expected. specificity factor30 (CPSF30-L) in m 6 A-dependent polyadenylation site selection, which in turn affect transcripts stability (Song et al., 2021;Hou et al., 2021). Using the alternative polyadenylation site detection tool TAPAS on our mRNA-seq data, we found 711 transcripts that show differential APA site usage between cold and control conditions (adjusted P < 0.05) out of 24 527 total transcripts. Interestingly, the GO analysis of these 711 transcripts show a strong enrichment for those encoding proteins involved in abiotic stress responses, including cold stress (Figure S10, Table S2). In total, 69 of these transcripts, including the one encoding the known cold stress response protein TCF1, are enrchied for m 6 A, as well as undergo alternative polyadenylation during cold stress (Ji et al., 2015). Other cold responsive transcripts such as GRP7, COR15A, COR15B, LTI78, KIN2, and TCF1 were also enriched for m 6 A as seen in the genome browser views ( Figure 3C, Figure S11). Notably, there was a distal to proximal polyA site usage shift in these transcripts under cold stress ( Figure S12). Upon closer inspection at the transcript level, the m 6 A reader CPSF30-L was not altered under cold stress. This suggests that the association between m 6 A and APA under cold stress is not based on changing the transcript levels of CPSF30-L. Additionally, controlenriched and shared m 6 A containing transcripts were also well-represented in this group of APA transcripts; thus, a m 6 A-dependent mechanism affecting polyA site usage could not be established, but a closer analysis of these findings is worth considering in the future.

The mta mutant plants show hypersensitivity to cold stress
To further assess the role of m 6 A methylation in cold stress responses, a post-embryonic mutant for MTA (mta ABI3: MTA; hereafter referred to as mta) (Bodi et al., 2012), the gene coding for a core protein in the m 6 A writer complex, was evaluated for cold tolerance. Because the complete loss of function for MTA gene is embryonic lethal, this mutant allows for embryonic expression and shows highly reduced levels of MTA transcripts in adult plants, as also shown by our qPCR analysis ( Figure 4a). Interestingly, MTA transcript levels were found to be increased in both Col-0 and mta in response to cold treatment, although the levels remained comparatively very low in mta mutant plants under both conditions ( Figure 4a). Thus, there was a significant loss of this writer protein in the mutant system. Furthermore, as expected, the overall m 6 A levels in total RNA was reduced significantly in the mta mutant compared to Col-0 plants ( Figure S13). To address whether disrupting the m 6 A deposition pathway affects cold stress response, we performed cold tolerance assays on 3-dayold Col-0 and mta seedlings that were transferred to vertical plates containing MS agar and placed at 4°C for 55 days (cold stress) and subsequently the root growth was measured. For untreated control plants, the seedlings were grown at 22°C for 10 days only and then the root growth was measured and photographs were taken. The primary root growth in the mta mutant was significantly retarded compared to the Col-0 plants in response to cold treatment (Figures 4b,c). Additionally, the biomass of mta mutant seedlings was also significantly reduced compared to Col-0 plants (Figure 4d). Cold stress is often closely correlated with greater ROS accumulation resulting in oxidative stress in the affected plants. Thus, ROS accumulation was assessed in the leaves collected from 25-day-old pot grown plants exposed to 4°C for 2 days and stained with NBT (Figure 4e). Both Col-0 and mta mutant leaves showed substantial staining under cold treatment, although the intensity was stronger in the mta mutant leaves suggesting a greater accumulation of ROS in the absence of m 6 A (Figure 4e). Consistant with our observation, increased ROS accumulation in m 6 A loss-of-function mutant (vir-1) under salt stress has been recently reported (Hu et al., 2021). Cold-acclimation is an adaptive process, which increases the freezing tolerance of plants. Without cold acclimation, plants exposed to freezing temperatures could be severely injured as a result of ice crystal formation in the apoplast causing membrane damage (Armstrong et al., 2020;Zhao et al., 2016). Electrolyte leakage assays have been previously used to assess the differences in freezing tolerance of Arabidopsis plants. To assess the importance of m 6 A RNA methylation in freezing tolerance, the electroylate leakage was determined in detached leaves of nonacclimated and cold-acclimated Col-0 and mta mutant plants. Accordingly, 4-week-old plants were exposed to freezing temperatures of -6 and -8°C directly without previous cold exposure (non-acclimated plants). This analysis revealed that Col-0 and mta leaves demonstrated significant differences (i.e. electrolyte leakage was far greater in the mutant compared to Col-0), suggesting that the mta mutant plants are hypersensitive to cold stress (Fig. 5). However, at higher freezing temepraures (-10 and -12°C), the electrolyte leakage did not differ between Col-0 and mta plants. This, could be because these temperatures already caused more than 90% electrolyte leakage indicating severe membrane damage to both Col-0 and mta plants. On the other hand, when cold-acclmated plants were exposed to -8, -10, and -12°C, the electrolyte leakage was consistantly and significantly greater in the mta mutant compared to Col-0 plants ( Figure 5). These results suggest that mta mutant plants were undergoing greater membrane leakage during freezing stress compared to Col-0. Thus, both under non-acclimated and cold-acclimated conditions, the cold tolerance competence of the mta mutant was severely impaired compared to Col-0. These results provided strong evidence for the importance of the m 6 A epitranscriptome mark in plant freezing tolerance.
The CBF and COR gene expression levels in the mta mutant under cold stress The CBF regulon (i.e. CBF1, CBF2, and CBF3) genes are induced during cold treatment that in turn activates the expression of COR genes in Arabidopsis (Thomashow, 2010). As measured by qPCR assays, in both total RNA and polysomal RNA fractions, the transcripts levels for these three CBFs are increased in response to 24 h of cold stress, although the magnitude of increase differed greatly between the three CBFs. Specifically, we found that CBF2 is more strongly upregulated than CBF1 or CBF3 at 24 h of cold in Col-0 plants ( Figure 6). Interestingly, in Col-0, the m 6 A levels on CBF2 were greatly elevated but the methylation profiles for CBF1 or CBF3 were mostly unaffected in response to cold treatment (Figure 3c). Similar to CBF2, the m 6 A profiles of five different COR transcripts (COR78, COR15A, COR15B, COR27, and COR413PM1) were also highly increased in response to cold treatment in Col-0 plants (Figure 3c). Furthermore, the increased m 6 A levels on these transcripts also positively correlated with their transcript abundances both in total RNA and polysome RNA fractions in cold-treated Col-0 compared to untreated controls (Figure 3c and Figure 6).
To gain insights into the effect of m 6 A on CBF transcritps, we compared their expression profiles in both mRNA and polysome RNA fractions of wild-type (WT) and mta mutant exposed to cold stress for 24 h. Interestingly, CBF2 levels were greatly induced in the mta mutant (i.e. up to 25-fold in total RNA and approximately 40-fold in polysome RNA fractions), whereas their levels were only increased by < 10-fold in Col-0 under cold treatment (Figure 6). Likewise, CBF1 levels were also strongly elevated in the mta mutant compared to WT under cold stress (Figure 6). The CBF1 and CBF2 transcripts are known to be induced highly after 3 h of exposure to low temperatures (Gilmour et al., 1998;Park et al., 2015;Novillo et al., 2004Novillo et al., , 2007. Therefore, we also measured CBF levels in Col-0 and mta mutant plants after 3 h of cold treatment. As expected, both CBF1 and CBF2 transcript levels were more strongly induced at 3 h compared to the 24 h of cold stress in Col-0 ( Figure S14). Compared with Col-0 plants, even after 3 h of cold stress, the CBF1 and CBF2 levels were more strongly increased in the mta mutant ( Figure S14). These observations from the mta mutant suggests a potential negative role for m 6 A (i.e. m 6 A negatively affects the levels of CBF2 and CBF1 both in mRNA and polysome fractions in Col-0 plants subjected to cold stress). On the other hand, four of the six COR transcripts (COR15A, COR15B, KIN1, and COR47 transcript levels in the mRNA fractions) were significantly lower under cold stress in the mta mutant compared to WT Col-0 ( Figure 6). Overall, the COR mRNAs in the ploysome RNA were reduced or had no change in the mta mutant, suggesting that a lack of m 6 A modifications on these transcripts does not significantly affect their ribosome loading. Taken together, compared to Col-0, the mta mutant impaired in m 6 A methylation demonstrated cold sensitivity (lower biomass and suppressed root growth) and greater oxidative stress coupled with the differences in CBF and COR gene expression profiles during cold stress. Most importantly, freezing tolerance of both cold-acclimated and nonacclimated mta plants was severely diminished compared to the Col-0 plants. Thus, this epitranscriptome modification appears to be critical for proper cold stress response in plants.

DISCUSSION
At the molecular level, low temperature alters the expression of thousands of genes in plants, which in turn contribute to adaptation to stress. Recently discovered m 6 A, m 5 C, and many other mRNA modifications are known to either stabilize transcripts or promote their degradation, thus revealing an additional layer of post-transcriptional regulations. Indeed, m 6 A has been implicated in various RNA regulatory processes including transcript stability and promoting translation in both plants and animals Huang et al., 2018;Park et al., 2019). However, the influence of m 6 A on stress responsive transcriptomes is still poorly understood, particularly in plants.
Our analyses indicated that the majority of m 6 A peaks were localized near the stop codon and 3 0 UTR regions of mRNA transcripts, in agreement with previous studies Kramer et al., 2020;Luo et al., 2014;Luo Figure 6. qPCR analysis suggests potential correlations between cold-enriched m 6 A containing transcripts and RNA abundance, as well as ribosome occupancy. Changes in relative mRNA abundance and polysome loading for candidate CBF and COR genes measured as fold change between control and cold conditions in Col-0 and mta.  Zhou et al., 2015). Although the overall distribution of m 6 A showed a similar 3 0 UTR bias in both cold and control conditions, when separating the transcripts based on their enrichment of m 6 A under stress, transcripts with cold and control-enriched m 6 A peaks showed a larger percentage of modification peak localization in their upstream regions (CDS and 5 0 UTR) and lower probability of being localized in the 3 0 UTR compared to the common m 6 A peaks that were identified in both conditions. This suggested that the stress sensitive dynamic peaks are more likely to deviate from the 3 0 UTR location compared to those that remain unaltered. Among the dynamic peaks, cold responsive m 6 A peaks had a relatively higher probability of being present in the 5 0 UTR region. Previous studies have suggested a similar trend under heat stress in human cell lines where transcripts show an increase in 5 0 UTR m 6 A (Zhou et al., 2015). In Arabidopsis, a slight shift in m 6 A peak localization was observed under salt stress relative to the peaks in control samples Kramer et al., 2020). The mechanism and the importance behind this shift is unknown although the stress association appears to be evident. Similar to several other plant studies, we did not find the DRACH motif that is described as a canonical m 6 A site in mammalian systems suggesting this motif may be less relevant in plants, especially in the context of cold stress. Instead, our m 6 A data showed an enrichment of a UGUA motif that has previously been published as a major plant m 6 A motif able to bind to ECT2, a plant m 6 A reader protein (Wei et al., 2018).
Our data indicate that transcripts containing coldenriched m 6 A peaks show a higher overall mRNA abundance compared to transcripts that lose this enrichment under cold stress or those that remain unchanged in their methylation status under stress. This suggests that m 6 A gain may result in an increase in transcript abundance of this subset of transcripts under cold stress. Indeed, previously, it was shown that salinity-responsive m 6 A peak containing transcripts of Arabidopsis were more stable under salt stress . There are other potential mechanisms leading to increased transcript abundance such as possible interactions of m 6 A writer or reader complexes with histone modifications to affect RNA transcription (Huang et al., 2009;Shim et al., 2020). Based on our data alone, it is hard to predict whether this change in mRNA abundance of the cold-enriched m6A containing transcripts is the result of an increased stability or enhanced transcription, or a combination of both, and this needs further investigation. Indeed, stabilization and destabilization of m 6 A containing transcripts via m 6 A reader proteins has been reported in mammalian system , suggesting that such consequences of m 6 A in RNA fate are context-dependent.
It was previously demonstrated that m 6 A plays an important role in the translation of mRNAs through a pathway involving m 6 A reader proteins. In mammals, the small rRNA subunit has an m 6 A reader YTHDC2 binding site near the mRNA entry site, suggesting that m 6 A can help facilitate efficient translation of transcripts Mao et al., 2019). In plants, less is known about how m 6 A affects the translation, particularly in response to stress conditions. Our ribosome profiling data suggest that, on average, transcripts that are specifically enriched for m 6 A under cold show a significantly more positive change in their ribosome occupancy, implying that these transcripts are more likely to be translated compared to transcripts containing cold-depleted m 6 A peaks. Upon further parsing of the cold responsive m 6 A containing transcripts, those that are truly driving this positive change in ribosome occupancy are associated with GO terms involving rRNA processing and stress response. By contrast, those that show decreased occupancy are largely involved in photosynthesis and development. Intriguingly, transcripts directly involved in cold stress showed both increases and decreases in ribosome occupancy, although with a higher enrichment for the group showing increased ribosome occupancy. These results suggests there is a potential role for m 6 A in promoting translation for specific transcripts encoding important cold stress response regulators, including CBF1 and CBF2. However, it is important to note that not all transcripts showing cold-enriched m 6 A demonstrate this trend. The present study suggests a complex and probably context dependent regulation that dictates which of the methylated transcripts shows a higher ribosome occupancy and possibly a higher rate of translation.
The positive correlations between the cold-enriched m 6 A containing transcripts and increased mRNA abundance and polysome binding are noteworthy observations. These positive correlations could simply be a result of the overall greater abundance of the transcripts, leading to increased m 6 A detection and polysome association. However, our observation of numerous transcript examples that did not support this being a general trend in our data argues against the notion that high levels of transcript abundance are a cause of the observed positive correlations (Figure S9). This could suggest a potential biological importance for the overall positive correlations between m 6 Aenrichment during cold stress response and the association of these transcripts with the polysomes.
Recent studies have identified a role for m 6 A reader protein CPSF30-L in mRNA polyadenylation in Arabidopsis (Hou et al., 2021;Song et al., 2021). The larger isoform of CPSF30-L is a plant specific protein. It contains the m 6 A reader YTH domain that can bind methylated transcripts and alter the polyadenylation site selection giving rise to alternative polyadenylation, in turn generating multiple mRNA isoforms that could alter RNA stability, as well as translation (Hou et al., 2021, Song et al., 2021. The mRNA transcript level of CPSF30L did not chnage under cold stress, suggesting that the regulation is not dependent on changing CPSF30L abundances. However, alternative polyA site usage analysis suggested that some of the major cold responsive transcripts with m 6 A methylation show changes in APA site usage. Interestingly, an m 6 A reader protein ECT4 shows a significant increase (2.3-fold) in transcript abundance under cold stress conditions. Although we could not demonstrate a link between coldinduced m 6 A and alternative polyadenylation, this is a possibility that requires further attention.
CBF and COR transcript levels in WT and mta plants under cold stress CBFs/DREB1s are the best characterized transcription factors known to be important for cold tolerance in plants and they alone regulate approximately 20% of the coldinducible genes in Arabidopsis (Zhao et al., 2016). Based on cbf1, cbf2 and cbf3 individual mutant analysis, CBF2 was suggested to play a more important role under cold stress than CBF1 or CBF3 (Zhao et al., 2016). Notably, of the three CBFs, m 6 A profiles on CBF2 were strongly elevated under cold stress ( Figure 3c). Intriguingly, the observation that CBF2 and CBF1 levels were significantly elevated in both total RNA and polysome fractions of the mta mutant compared to Col-0 plants both at 3 h and 24 h of cold treatment suggests that the m 6 A mark in WT plants could negatively influence these CBFs (Figure 6 and Figure S14). Such a negative role could be part of a finetuning mechanism of CBF2 levels under stress conidtions (i.e. after a threshold level of transcript abundance/translation, both the CBF2 transcript stability and its loading onto the polysomes were negatively regulated by the presence of m 6 A in the Col-0 plants). How m 6 A modification can bring about dual but opposite effects (i.e. promoting stability of certain coldregulated transcripts at the same time as destabilizing other transcripts) (CBF2 and CBF1) is intriguing. Such contrasting outcomes are likely dependent on the nature of the reader protein, which recognizes the specific methylated sites. For example, in mammalian systems, recognition of m 6 A containing transcripts by YTHDF2 and YTHDF3 reader proteins promotes destabilization. Conversely, if the m 6 A peak containing transcripts are recognized by IGF2BP or Tudor SND1 (m 6 A reader proteins), then the transcripts are stabilized (Baquero-Perez et al., 2019;Huang et al., 2018). These and other possibilities are worthy of further exploration in the context of plant cold stress response.
Despite the highly elevated CBF2 and CBF1 levels in the mta mutant, the downstream target COR gene (COR47, COR15B, COR15A, and KIN1) expression levels were low or unaffected in the mRNA fractions and the COR78/RD29a and COR15A levels were significantly lower in the polysome RNA fractions compared to Col-0. It was recently reported that the transcripts lacking m 6 A in mta mutant plants were relatively less stable under salt stress . Specifically, the endonuclease-mediated cleavage of transcripts was far greater in the mta mutant in which m 6 A levels on mRNAs were significantly decreased. By contrast, m 6 A-enriched saltresponsive transcripts in Col-0 with normal m 6 A levels are prevented from such endonucleolytic degradation, and, consequently, the transcripts were more stable during salt stress Kramer et al., 2020). Thus, m 6 A is likely to have important roles in regulating the plant stress responsive transcriptome across the gambit of abiotic stressors that plants encounter.
Plants subjected to low temperature stress undergo massive physiological adjustments largely as a result of the large-scale change in gene expression profiles regulated at various levels including before and after transcription of the RNA. Using high-throughput sequencing strategies, the present study demonstrates that enrichment of m 6 A under cold stress leads to a general increase of those mRNA abundances. Furthermore, transcripts that show m 6 A enrichment under cold stress also demonstrates an increase in ribosome occupancy, which is a proxy for translation, suggesting that the m 6 A mark plays an important role during plant cold stress responses. Consistent with these observations, the mta mutant was found to be hypersensitive to cold stress. This was accompanied by a greater ROS accumulation coupled with the reduced expression of some of the COR genes, as well as diminished cold acclimation capacity of the mta mutant. All of these observations corroborate the suggestion that m 6 A RNA methylation plays a critical role in mounting the appropriate molecular and physiological responses that make the plant adapt to cold stress.

EXPERIMENTAL PROCEDURES
Stress treatment for sequencing MeRIP and polysomal RNA and mRNA profiles and gene expression analysis using a quantitative real-time PCR (qRT-PCR) Arabidopsis thaliana (Col-0) seeds were surface sterilized and plated on full strength MS agar medium with 1% sucrose. The plates were kept at 4°C for 2 days for vernalization and transferred to a growth chamber maintained under a 12/12 h light/dark photocycle at 22°C. For cold treatment, 3-week-old seedlings were transferred to 4°C for 24 h with a low-level light intensity (35 μmol photons m −2 sec −1 ) to avoid photoinhibition (Fowler & Thomashow, 2002), whereas control seedlings were maintained at 22°C. At the end of the treatment, both control and cold-treated seedlings were harvested and stored at −80°C.

MeRIP-seq
Total RNA was extracted from control, and cold stressed seedlings using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The quantity and quality total RNA were measured using a Nanodrop ND-1000 (Thermo Fisher Scientific) and 1.2% agarose gel, respectively. Approximately 120 μg of total RNA per replicate was used for poly(A) RNA purification using a Seq-Star TM poly(A) mRNA Isolation Kit (Arraystar Inc., Rockville. MD, USA). The purified mRNA was chemically fragmented in 20 μl of fragmentation buffer (10 mM Zn 2+ , 10 mM Tris-HCl, pH 7.0) at 94°C for 5 min. The size of the fragments was analyzed using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). An aliquot of the fragmented mRNA was kept as the input control, and the remaining fragmented mRNA was immunoprecipitated with 2 μg of anti-m 6 A rabbit polyclonal antibody (Synaptic Systems) in a 500-μL immunopreciptation (IP) reaction at 4°C for 2 h. 20 μl Dynabeads TM M-280 Sheep Anti-Rabbit IgG (Thermo Fisher Scientific) was prepared by blocking with 0.5 mg ml -1 BSA at 4°C for 2 h. The mRNA-antibody IP reaction mixture was incubated with the blocked Dynabeads for an additional 2 h at 4°C. Three washes with 1 × IP buffer (10 mM Tris-HCl at pH 7.4, 150 mM NaCl, 0.1% NP-40) and two washes with 1 × wash buffer (10 mM Tris-HCl at pH 7.4, 50 mM NaCl, 0.1% NP-40) were performed. The m 6 A-antibody immunoprecipitated mRNA fragments were eluted from the Dynabeads in 200 μl of elution buffer (10 mM Tris-HCl at pH 7.4, 1 mM EDTA, 0.05% SDS, 40 U of proteinase K) at 50°C for 30 min. The immunoprecipitated mRNA fragments were extracted by phenolchloroform and precipitated with ethanol and dissolved in nuclease-free water. The m 6 A antibody enriched m 6 A containing mRNA fragments and the input samples were used to construct the RNA-seq library using a KAPA Stranded mRNA-seq Kit (Illumina ® platform; Illumina Inc., San Diego, CA, USA). Each treatment had three biological replicates. The quality of constructed libraries was analyzed using an Agilent 2100 Bioanalyzer and subsequently sequenced on the Illumina HiSeq 4000 platform at the Arraystar in accordance with the HiSeq 3000/4000 SBS Kit protocol.
Validation of differentially expressed genes at transcriptional and translational levels using qRT-PCR qRT-PCR was performed to validate the differentially expressed genes identified in RNA-seq and polysome RNA-seq, as well as to determine transcript levels of CBF and COR genes in Col-0 and mta mutant. The first-strand cDNA was synthesized with 2 μg of total RNA using Superscript reverse transcriptase II (Invitrogen, Waltham, MA, USA). PCR reactions were performed in the total volume of 10 μl, with 0.5 μM of each forward and reverse primers and Power 2 × SYBR Green PCR mix (Applied Biosystems, Waltham, MA, USA) on a LightCycler 96 system (Roche, Basel, Switzerland). The PCR conditions included an initial denaturation at 96°C for 1 min, followed by 40 cycles of 5 sec at 94°C and 1 min at 60°C. The PCR was performed on two independent biological samples with two technical replicates. The relative expression levels of the target genes were calculated using 2 -ΔΔct (Livak & Schmittgen, 2001). Protein phosphatase 2A subunit A3 (PP2AA3) was used as reference gene for normalization.

Bioinformatics analysis of m 6 A-seq and polysomal-seq datasets
Read processing. Raw fastq files obtained from high throughput RNA sequencing were first checked for read quality using the java application FASTQC (Andrews, 2010). After validating that the reads were of satisfactory quality, the Illumina sequencing adapters were trimmed from raw reads using the tool TRIMMOMATIC-0.36 (Bolger et al., 2014). The trimmed fastq files were then mapped to Arabidopsis genome assembly TAIR10 using STAR ALIGNER (Dobin et al., 2013) with default parameters.
Differential RNA abundance. The bioinformatic package HTSEQ (Anders et al., 2015). was used to count raw reads for each transcript using the TAIR10 reference annotation file obtained from the ENSEMBL repository (https://www.ensembl.org/index.html). These raw counts were then used as input for the R statistical package DESEQ2 (Love et al., 2014) with the default settings, which uses normalization factors that incorporate library depth and gene wide dispersions for normalization of the reads and determines differentially expressed genes along with respective P-values. Transcripts that showed an adjusted P < 0.05 were considered statistically significant. m 6 A peak calling. To identify m 6 A peaks, the peak calling software MACS2 was used on the alignment bam files from STAR using the parameters as described previously . Because MACS2 (https://pypi.python.org/pypi/MACS2) is blind to strand information, we first split bam files using SAMTOOLS (Li et al., 2009) based on strand information and called peaks separately for each strand. The background file used for calling the peaks was the input mRNA sequencing data, therefore yielding peaks that were significantly enriched in the pulldown compared to all expressed mRNAs. The peak regions had lengths ranging from 50 to 3000 bp, with a median length of 249 bp. Peaks were filtered using P < 0.05. Validation of several peaks were performed using MeRIP-qPCR. To obtain transcript level information, peak coordinates were overlapped with annotated transcript loci in a strand-specific manner. The longest transcript was taken as the reference when a gene had multiple annotated isoforms to prevent repeated counting of the same peak.
Motif searches. The de novo motif search algorithm HOMER (http://homer.ucsd.edu/homer/motif/) was used to identify the enriched motifs in the identified peak regions. In simple words, HOMER analyzes for enriched motifs in given sequence regions against randomized background of nucleotides with a matched GC% content along with the possible false positive ratio. The input files for HOMER were the peak bed files containing positional information for the MACS2 (https://pypi.python.org/pypi/MACS2) output peaks and the genome fasta file for Arabidopsis. The Plant Journal, (2022), 111, 1052-1068 exon, and 3 0 UTR. If more than one 5 0 UTR annotation was available for a transcript, the longest and the most upstream one was used, and, if more than one 3 0 UTR annotation was available, the longest and the most downstream one was used. The bed peaks for each category were then aligned to the trimmed GFF. After aligning, the most downstream end of the 5 0 UTR and the most upstream end of the 3 0 UTR were set as index (position) 0 for the 5 0 UTR and 3 0 UTR graphs, respectively. Finally, the probability of m 6 A for each category and graph was calculated as: where s(i) and T(i) represent the number of transcripts having an m 6 A site at index i and the number of transcripts having index i, respectively; the max term represents the maximum probability of m 6 A over the indices ranging from (inclusive) −n to n; and argmax represents the index at which the probability of m 6 A is maximized within the given indices.

Ribosome occupancy.
To determine ribosomal occupancy of each transcript, the read per million (RPM) value was first calculated for the transcripts using raw transcript read count to normalize for library depth. As previously described, ribosome occupancy for each transcript was then calculated as: ribosome occupancy = log2 (RPM from polysome sequencing/ RPM in the mRNA-seq) To determine how the ribosome occupancy of a transcript changes under cold stress, we first determined the ribosome occupancy for the transcripts in control and cold samples separately. Then, we measured the change in ribosome occupancy between conditions by taking the difference of the log transformed values where a positive value would indicate an increased ribosome occupancy under cold treatment and a negative value would indicate a reduced ribosome occupancy compared to the untreated control sample (Ribo occupancy change = Ribo Occupancy cold -Ribo Occupancy control ).
APA analysis. APA analysis was carried out using the bioinformatic tool TAPAS on mRNA-seq data as described previously (Arefeen et al., 2018). Briefly, de novo APA sites were detected from the mRNA-seq data across all conditions and combined with previously annotated polyA + site in the TAIR10 database using ENSEMBL annotation. Differential usage of each polyA + site was measured across each condition and replicates to give a list of transcripts with statistically significant differential APA usage. Furthermore, normalized read counts for each polyA + site were given as an output by the tool, which was then used to calculate the proportion of distal to proximal polyA + site usage for candidate transcripts.
Dot blot assay. Approximately 200 ng of total RNA obtained from Col-0 and mta mutant plants was dispensed on an N + Hybond membrane (GE Healthcare, Chicago, IL, USA) and crosslinked twice using a Stratalinker 2400 (Stratagene, San Diego, CA, USA) at 1200 μJ for 30 sec. The membrane was washed with PBS-Tween 20, blocked in 5% milk and incubated with 1:1000 primary antibody (SYSY m 6 A antibody) overnight at 4°C. After three washes, the membrane was incubated in 1:10000 horseradish peroxidase cojugated secondary antibody for 1 h, washed three times and visualized using ECL reagent (Pierce, Rockford, IL, USA).
Cold tolerance assay. Arabidopsis thaliana (Col-0) and mta mutant seeds were surface sterilized and plated on full strength MS agar plates containing 1% sucrose. The plates were kept at 4°C for 2 days for vernalization and then transferred to a growth chamber maintained under a 12:12 h light/dark photocycle at 22°C. For root growth assay, 3-day-old Col-0 and mta mutant seedlings were transferred to a vertical MS agar plates, which were maintained at 4°C with a low-level light intensity (35 μmol m −2 sec −1 ), whereas seedlings that were maintained at 22°C served as controls. The primary root length and fresh weight of seedlings were recorded after 52 days of cold treatment. The untreated seedlings are only 13 days old (3 day-old seedlings were transferred to a vertical plate and maintained at 22°C for 10 days and primary root length and fresh weight were recorded and also photographed. Detection of ROS. For superoxide detection, leaves were detached from 25-day-old pot grown plants that were grown under control and cold conditions (4°C for 2 days) and immersed in 50 mM sodium phosphate buffer (pH 7.5) containing 1 mg ml -1 nitroblue tetrazolium (NBT) (N6876; Sigma-Aldrich, St Louis, MO, USA). After 12 h of incubation in the dark at room temperature, samples were then boiled in acetic acid:glycerol:ethanol (1:1:3) for 10 min and stored in 95% ethanol until photographs were taken. The average signal intensity across the leaf area was measured using ImageJ (NIH, Bethesda, MD, USA).
Freezing tolerance in cold-acclimated and non-acclimated WT and mta mutant plants. Electrolyte leakage assay was used to assess the differences in cold tolerance of Col-0 and mta plants as described previously (Zhao et al., 2016) with some minor modifications. Briefly, 4-week-old plants grown in soil under a 12:12 h light/dark photocycle at 22/20°C were divided into two groups; one group was transferred to 4°C for 7 days for cold acclimation, whereas the other group continued to grow under a 12:12 h light/dark photocycle at 22/20°C and served as the control group (non-acclimated). For electrolyte leakage analysis of coldacclimated and non-acclimated Col-0 and mta plants, seven replicates of two fully developed rosette leaves were excised and placed in 20 × 150 mm tubes containing 100 μl of deionized water. The tubes were placed in a freezing growth chamber (model PGR15; Conviron. Winnipeg, MB, Canada) at 0°C for 30 min to equilibrate. Then, the temperature was decreased to -12°C at a rate of 1°C per 30 min. When the temperature reached -1°C, tiny pieces of equal volume of ice were added to each tube for ice nucleation. The tubes were removed from the chamber at -6, -8, -10, and -12°C and immediately placed on ice for gradual thawing. Then, to each tube, 12 ml of deionized water was added and the tubes were gently shaken overnight at room temperature. The conductivity of solutions was measured using a conductivity meter (Orion Star A212; Thermo Fisher Scientific). Next, these tubes along with the samples were autoclaved at 121°C for 20 min, and then agitated for 2 h before the conductivity was measured again. The percentage of electrolyte leakage was calculated as the ratio of the conductivity before compared to that after autoclaving.

ACCESSION NUMBER
All raw and processed data have been uploaded to NCBI GEO under accession number GSE184056. Science Foundation, USA grant IOS-2023310 to BDG, IOS-1849708 to BDG and RS, and NSF EPSCoR Track II (1826836) as well as Stevens Endowed Chair to RS. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

AUTHOR CONTRIBUTIONS
RS conceived and designed the study. GG, Y-FL, and PM performed the experiments. BS, CA, JSR, and BDG analyzed the datasets. RS, BS, JSR, and BDG interpreted the data and wrote the manuscript with inputs from all of the coauthors.

CONFLICT OF INTEREST
The authors declare no conflict of interest.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article. Figure S1. m 6 A peak validations under cold stress. Figure S2. m 6 A peaks identified in the present study show high overlap with previously published m 6 A peaks. In addition, thousands of new peaks identified in the present study. Figure S3. Average relative read coverage of input mRNA seq compared to m 6 A IP reads across 100-nucleotide binned transcript in (a) control and (b) cold conditions. Figure S4.  Figure S5. (a) Comparison of mRNA abundance for cold-enriched m 6 A containing transcripts methylation in the 5 0 UTR and start codon regions compared to those downstream. (b) Change in mRNA abundance of 5 0 UTR and start codon m6A containing transcripts enriched in cold, control, or shared between conditions. Figure S6. (a) Change in transcript levels of highly upregulated genes under cold stress as identified by DESeq analysis. Protein IDs reveal important cold regulated proteins are upregulated. FC and SEM were calculated using DESEQ2. (b) Validation of candidate genes that show upregulation in transcript level in RNA-seq. Figure S7. Gene Ontology functional categories and -log(P-value) for transcripts that are enriched for m 6 A in control and cold conditions separated by a positive or negative change in mRNA abundance. Figure S8. Average ribosome occupancy across all transcripts calculated for cold treated and control samples. Figure S9. Example browser views of transcripts that do not show positive correlation between mRNA-seq, m 6 A-seq, and polysomal-RNA-seq seen in Figure 3c. Figure S10. Gene Ontology functional enrichment for 711 transcripts identified by TAPAS as differentially expressed for alternative polyadenylation site usage. Figure S11. Browser views of (a) TCF, (b) CPSF30, and (c) ECT4 transcripts comparing changes across all three libraries under cold stress. Figure S12. (a) Proportion of reads belonging to either proximal (P) or distal (D) polyadenylation sites for cold responsive transcripts with two identified polyA sites. (b) Proportion of reads belonging to either proximal (P), medial (M), or distal (D) polyadenylation sites for cold responsive transcripts with three identified poly A sites. Transcript data can be found in Table S2. Figure S13. Dot-blot showing the reduced levels of m 6 A in the mta mutant with and without cold stress. Figure S14. Cold-enriched m 6 A containing CBF1 and CBF2 transcript abundances in Col-0 and mta mutant after 3 h of cold treatment. (a). RT-qPCR analysis of CBF1 and CBF2 transcript levels after 3 h of cold stress. (b). m 6 A levels of CBF1 and CBF2 transcripts in Col-0 and mta after 3 h of cold stress. Table S1. List of genes with m 6 A peaks