mRNA N6 -methyladenosine is critical for cold tolerance in Arabidopsis.

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, m6 A, m5 C 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 m6 A in plant growth and development has been appreciated, although its significance under stress conditions is still underexplored. To assess the role of m6 A 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 m6 A analysis revealed large-scale shifts in this modification in response to low temperature stress. Because m6 A is known to affect transcript stability/degradation and translation, we investigated these possibilities. Interestingly, we found that cold-enriched m6 A-containing transcripts demonstrated the largest increases in transcript abundance coupled with increased ribosome occupancy under cold stress. The significance of the m6 A epitranscriptome on plant cold tolerance was further assessed using the mta mutant in which the major m6 A 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 non‐coding RNAs, and microRNAs (Dominissini et al., 2012). Among the mRNA modifications, N 6‐methyladenosine (m6A) 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 m6A is added to an RNA by a large, multi‐component protein ‘writer complex’ containing two methyltransferases [METHYLTRANSFERASE 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 m6A 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 m6A 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, m6A 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 m6A have largely been determined from the analysis of mutants that are defective in writer, reader and eraser complexes. For example, the disruption of Arabidopsis m6A 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 (Arribas‐Hernandez et al., 2020). In rice, m6A 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). m6A has also been shown to play roles in fruit ripening in tomato (Zhou et al., 2019). Furthermore, m6A methylome studies conducted in various tissues of Arabidopsis and rice revealed tissue‐specific changes in m6A levels, which in turn affects gene expression (Li et al., 2014; Wan et al., 2015). These studies illustrate the importance of m6A in plant growth and development but their roles in plant adaptation to stresses are only emerging. Only recently has the potential involvement of m6A in salt stress responses been reported in Arabidopsis (Anderson et al., 2018; 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 m6A 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, 2018). 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 (Shi et al., 2018). 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., 2017, 2018; Thomashow, 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; Shi et al., 2018). 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 post‐transcriptional regulatory mechanism could also play important roles in plant adaptation to low temperature stress. To investigate the importance of m6A mRNA modifications on cold stress responses, we generated methylated RNA immunoprecipitation (MeRIP)‐sequencing (seq), mRNA‐seq and polysomal RNA‐seq. The significance of m6A 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 m6A epitranscriptome modification and cold tolerance in plants.

Results

MeRIP‐seq identifies transcripts that are m‐enriched or m‐depleted under cold stress

To understand the dynamics of m6A methylation in the overall transcriptome under low temperature stress, MeRIP‐seq was performed on polyA+ selected RNA in 3‐week‐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 peak‐calling algorithm macs2 was employed to identify m6A 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 high‐confidence m6A 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 m6A peaks were determined based on the peak abundances in the control and cold‐treated samples. This analysis identified approximately 1318 and 5829 m6A peaks that were designated as either cold‐enriched or control‐enriched (cold‐depleted) m6A peaks, respectively, and the remaining 11 943 m6A peaks shared between the two conditions, implying that their methylation profiles did not change upon cold stress (Figure 1a, Table S1). Cold‐enriched m6A 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 m6A peaks) and cold‐depleted (5829 m6A peaks) are represented by 1198 cold‐enriched and 4545 control‐enriched (cold‐depleted) m6A containing transcripts, respectively. Among these, 233 transcripts contain both a control‐enriched and a cold‐enriched m6A 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.

Figure 1

Arabidopsis m6A shows a 3′ UTR positional bias, although the bias is stronger for cold unresponsive peaks. (a) High‐confidence m6A peaks obtained by overlapping biological replicates of cold and control treated samples revealed 5829 control‐enriched (cold‐depleted) peaks, 1318 cold‐enriched m6A 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 m6A peak regions. (d) Percentage distribution of cold‐enriched, control‐enriched, and shared m6A peaks that overlap with the 5′ UTR, start codon, CDS, stop codon, and 3′ UTR. (e) Nucleotide level probability metric for methylation in the start and stop codon boundary regions for cold‐enriched, control‐enriched, and shared m6A peak containing transcripts. Each group is normalized to its maximim value in the region. [Colour figure can be viewed at wileyonlinelibrary.com]

To validate the MeRIP‐seq profiles, several cold‐enriched m6A containing transcripts (AT2G46690, AT4G19530, AT4G25990, AT5G17300, and AT5G62360) were analyzed using m6A‐immunoprecipitation quantitative polymerase chain reaction (qPCR) (Table S3). As expected, these transcripts were revealed to be more enriched in m6A pulldown fractions after cold treatment compared to untreated controls (Figure S1). Furthermore, we compared the identified m6A peaks in the present study with two previously published studies that reported m6A peaks from the adult leaves and cotyledons (Anderson et al., 2018; Shen et al., 2016). Unsurprisingly, the overlap between the m6A peaks in the present study and the previous studies was more than 88% (Figure S2).

To identify the localization of m6A peaks on the different positional segments of mRNA [5′ UTR, coding sequence (CDS) and 3′ untranslated region (UTR)], a relative coverage of the MeRIP‐seq library reads was plotted across transcript UTR and CDS positions. This revealed that the m6A pulldown reads were highly enriched in the 3′ UTR and near the stop codon of Arabidopsis transcripts, a typical feature of m6A 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 m6A localization identified in previous reports from multiple eukaryotes including plants (Anderson et al., 2018; Kramer et al., 2020; Luo et al., 2014; Luo et al., 2020, 2021; Zhou et al., 2015).

To identify motifs that are enriched for m6A 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 m6A consensus motif in plants, observed in several plant species including Arabidopsis (Hu et al., 2021; Wei et al., 2018), rice (Li et al., 2014), and maize (Luo et al., 2020). Other m6A motifs such as the AAGM and the GAD motifs were also found with high probability of a canonical ‘A’ as a site for m6A 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 peaks that are enriched or depleted under stress are less strictly localized in the 3′ UTR and have a relatively higher probability of occuring in CDS and 5′ UTR regions

Previous studies have suggested that the localization of m6A peaks on various mRNA segments are dynamic and can shift towards the 5′ 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 m6A peak regions identified in our cold‐enriched, control‐enriched, and shared m6A containing transcripts were mapped onto the annotated 5′ UTR, start codon, CDS, stop codon, and 3′ UTR (TAIR10 database). This analysis has revealed a small but noticeable increase in the percent of peaks that overlapped with the 5′ 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 m6A enrichment under cold stress, both the cold‐enriched and cold‐depleted m6A peak containing transcripts show a higher percentage of peaks in 5′ 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 m6A 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′ 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 m6A 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 control‐enriched (cold‐depleted) m6A 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 m6A 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 m6A peaks that are insensitive to cold stress tend to localize in 3′ UTRs, the more dynamic m6A peaks that are either enriched or depleted under cold have a relatively higher probability of being found in regions upstream of the 3′ UTR of the methylated transcript. Notably, the position of maximum probability of these groups, particularly for cold‐enriched m6A 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 cold‐enriched m6A 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 m6A in other locations, as well as 5′ UTR/start codon region m6A transcripts that were control‐enriched or shared (Figure S5b). These observations suggests a role for m6A in the 5′ 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, 2018). 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 stress‐altered transcriptome analysis in Arabiopsis (Fowler & Thomashow, 2002; Seki et al., 2001; Zhao et al., 2016).

Figure 2

Transcriptomic analysis after 24 h of cold stress shows increased mRNA abundance for major known cold responsive genes. (a) MA plot showing overall cold‐induced transcriptome changes in cold treated compared to control Col‐0 plants. Red points represent transcripts with statistically significant changes in RNA abundance (adjusted P < 0.05) and gray represents transcripts with no change. Blue points highlight transcripts known to be associated with the plant cold stress response. (b) Heatmap showing the GO terms associated with transcripts that are significantly more or less abundant during cold stress. (c) Fold change of mRNA abundances between cold treatment and control conditions compared between transcripts that are enriched for m6A peaks in cold and control, as well as shared between them. [Colour figure can be viewed at wileyonlinelibrary.com]

Transcripts containing cold‐enriched m peaks show increased abundance

Previous reports indicated a link between m6A and transcript stability in animal and plant cells (Anderson et al., 2018; Huang et al., 2018; Park et al., 2019; Kramer et al., 2020). To uncover any correlation between m6A enrichment and transcript abundance in the present study, the overall change in transcript abundances of the cold‐enriched, control‐enriched, and shared m6A 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 m6A containing transcripts showed a significantly higher increase in transcript abundance under stress compared to the transcripts that contain control‐enriched (cold‐depleted) and shared m6A peaks (P < 0.05 in all comparisons; Wilcoxon rank sum test) (Figure 2c). These observations suggest that transcripts that contain cold‐enriched m6A peaks have increased abundance. As expected, these cold‐enriched m6A 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 peaks show increased polysomal association

Previous studies have shown that m6A can act as a mark that promotes translation when recognized by members of the YTHDF2 reader protein family in mammalian systems (Mao et al., 2019; Wang et al., 2015). These studies have suggested an m6A reader‐dependent mechanism forming a link between increased ribosome loading of methylated transcripts and increased elongation. To assess whether m6A 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, 2018). 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 m6A peaks demonstrated a significantly more positive change under cold compared to transcripts containing cold‐depleted (control‐enriched) m6A 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 cold‐enriched m6A 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 m6A containing transcripts belonging to specific functional categories.

Figure 3

Cold‐enriched m6A peak containing transcripts show higher mRNA abundance and ribosome occupancy in response to cold treatment. (a) Ribosome occupancy differences for transcripts between control and cold treatments compared between transcripts enriched for m6A peaks in cold only and control only, as well as shared in both conditions. (b) Heatmap showing the breakdown of categories of Biological Function from a GO analysis of cold‐enriched m6A peak containing transcripts. (c) Genome browser views comparing meRIP‐seq (labeled m6A), mRNA‐seq (labeled mRNA), and polysome profiling (labeled as Ps‐seq) for major known cold regulatory transcripts. The arrows represent the direction in which the gene is transcribed. [Colour figure can be viewed at wileyonlinelibrary.com]

Major cold responsive ( ) 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 m6A 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 m6A detection, genome browser views were generated for several random transcripts. We found several examples where enrichment of m6A, 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.

Many cold‐enriched m containing transcripts undergo alternative polyadenylation events

mRNA undergoes polyadenylation which can affect its overall stability, localization, and translation (Di Giammartino et al., 2011). Recent studies have established a role for the larger isoform of cleavage and polyadenylation specificity factor30 (CPSF30‐L) in m6A‐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 m6A, 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 m6A 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 m6A reader CPSF30‐L was not altered under cold stress. This suggests that the association between m6A and APA under cold stress is not based on changing the transcript levels of CPSF30‐L. Additionally, control‐enriched and shared m6A containing transcripts were also well‐represented in this group of APA transcripts; thus, a m6A‐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 m6A 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 m6A 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 m6A levels in total RNA was reduced significantly in the mta mutant compared to Col‐0 plants (Figure S13). To address whether disrupting the m6A deposition pathway affects cold stress response, we performed cold tolerance assays on 3‐day‐old 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 m6A (Figure 4e). Consistant with our observation, increased ROS accumulation in m6A loss‐of‐function mutant (vir‐1) under salt stress has been recently reported (Hu et al., 2021).

Figure 4

Cold response in mta mutant compared to Col‐0 plants reveal differences in mRNA abundance and polysome loading. (a) Relative mRNA level of MTA transcripts in mta mutant compared to Col‐0 plants under control and cold conditions measured using qPCR. (b) Phenotypic assay and quantification of cold response of mta mutant compared to Col‐0 plants as measured for (c) root length (control seedlings were grown at 22°C for 10 days only, whereas, for cold treatment, the seedlings were grown for 52 days) and (d) fresh weight. (e) ROS accumulation in leaf tissues after long term chilling stress. [Colour figure can be viewed at wileyonlinelibrary.com]

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 m6A RNA methylation in freezing tolerance, the electroylate leakage was determined in detached leaves of non‐acclimated 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 m6A epitranscriptome mark in plant freezing tolerance.

Figure 5

The mta mutant is hypersensitive to freezing stress as assessed by an electrolyte leakage assay. (a) Electrolyte leakage assay performed on 5 week‐old non‐acclimated (NA) and (b) cold‐acclimated (CA, 4°C for 7 days) Col‐0 and mta mutant plants. Asterisks indicate statistically significant difference (P < 0.05, t test). [Colour figure can be viewed at wileyonlinelibrary.com]

The and 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 m6A 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 m6A 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 m6A 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).

Figure 6

qPCR analysis suggests potential correlations between cold‐enriched m6A 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. [Colour figure can be viewed at wileyonlinelibrary.com]

To gain insights into the effect of m6A 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., 2004, 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 m6A (i.e. m6A 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 m6A modifications on these transcripts does not significantly affect their ribosome loading. Taken together, compared to Col‐0, the mta mutant impaired in m6A 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 non‐acclimated 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 m6A, m5C, 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, m6A has been implicated in various RNA regulatory processes including transcript stability and promoting translation in both plants and animals (Anderson et al., 2018; Huang et al., 2018; Park et al., 2019). However, the influence of m6A on stress responsive transcriptomes is still poorly understood, particularly in plants.

Our analyses indicated that the majority of m6A peaks were localized near the stop codon and 3′ UTR regions of mRNA transcripts, in agreement with previous studies (Anderson et al., 2018; Kramer et al., 2020; Luo et al., 2014; Luo et al., 2020, 2021; Zhou et al., 2015). Although the overall distribution of m6A showed a similar 3′ UTR bias in both cold and control conditions, when separating the transcripts based on their enrichment of m6A under stress, transcripts with cold and control‐enriched m6A peaks showed a larger percentage of modification peak localization in their upstream regions (CDS and 5′ UTR) and lower probability of being localized in the 3′ UTR compared to the common m6A peaks that were identified in both conditions. This suggested that the stress sensitive dynamic peaks are more likely to deviate from the 3′ UTR location compared to those that remain unaltered. Among the dynamic peaks, cold responsive m6A peaks had a relatively higher probability of being present in the 5′ UTR region. Previous studies have suggested a similar trend under heat stress in human cell lines where transcripts show an increase in 5′ UTR m6A (Zhou et al., 2015). In Arabidopsis, a slight shift in m6A peak localization was observed under salt stress relative to the peaks in control samples (Anderson et al., 2018; 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 m6A site in mammalian systems suggesting this motif may be less relevant in plants, especially in the context of cold stress. Instead, our m6A data showed an enrichment of a UGUA motif that has previously been published as a major plant m6A motif able to bind to ECT2, a plant m6A reader protein (Wei et al., 2018).

Our data indicate that transcripts containing cold‐enriched m6A 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 m6A 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 m6A peak containing transcripts of Arabidopsis were more stable under salt stress (Anderson et al., 2018). There are other potential mechanisms leading to increased transcript abundance such as possible interactions of m6A 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 m6A containing transcripts via m6A reader proteins has been reported in mammalian system (Wang et al., 2014), suggesting that such consequences of m6A in RNA fate are context‐dependent.

It was previously demonstrated that m6A plays an important role in the translation of mRNAs through a pathway involving m6A reader proteins. In mammals, the small rRNA subunit has an m6A reader YTHDC2 binding site near the mRNA entry site, suggesting that m6A can help facilitate efficient translation of transcripts (Wang et al., 2015; Mao et al., 2019). In plants, less is known about how m6A affects the translation, particularly in response to stress conditions. Our ribosome profiling data suggest that, on average, transcripts that are specifically enriched for m6A 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 m6A peaks. Upon further parsing of the cold responsive m6A 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 m6A 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 m6A 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 m6A 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 m6A 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 m6A‐enrichment during cold stress response and the association of these transcripts with the polysomes.

Recent studies have identified a role for m6A 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 m6A 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 m6A methylation show changes in APA site usage. Interestingly, an m6A 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 cold‐induced m6A and alternative polyadenylation, this is a possibility that requires further attention.

and 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 cold‐inducible 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, m6A 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 m6A mark in WT plants could negatively influence these CBFs (Figure 6 and Figure S14). Such a negative role could be part of a fine‐tuning 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 m6A in the Col‐0 plants).

How m6A modification can bring about dual but opposite effects (i.e. promoting stability of certain cold‐regulated 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 m6A containing transcripts by YTHDF2 and YTHDF3 reader proteins promotes destabilization. Conversely, if the m6A peak containing transcripts are recognized by IGF2BP or Tudor SND1 (m6A 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 m6A in mta mutant plants were relatively less stable under salt stress (Anderson et al., 2018). Specifically, the endonuclease‐mediated cleavage of transcripts was far greater in the mta mutant in which m6A levels on mRNAs were significantly decreased. By contrast, m6A‐enriched salt‐responsive transcripts in Col‐0 with normal m6A levels are prevented from such endonucleolytic degradation, and, consequently, the transcripts were more stable during salt stress (Anderson et al., 2018; Kramer et al., 2020). Thus, m6A 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 m6A under cold stress leads to a general increase of those mRNA abundances. Furthermore, transcripts that show m6A enrichment under cold stress also demonstrates an increase in ribosome occupancy, which is a proxy for translation, suggesting that the m6A 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 m6A 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™ poly(A) mRNA Isolation Kit (Arraystar Inc., Rockville. MD, USA). The purified mRNA was chemically fragmented in 20 μl of fragmentation buffer (10 mm Zn2+, 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‐m6A rabbit polyclonal antibody (Synaptic Systems) in a 500‐μL immunopreciptation (IP) reaction at 4°C for 2 h. 20 μl Dynabeads™ 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 m6A‐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 phenol–chloroform and precipitated with ethanol and dissolved in nuclease‐free water. The m6A antibody enriched m6A 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.

Polysomal RNA sequencing (translatome profiling)

Polysomes were isolated from 3‐week‐old Arabidopsis seedlings grown under control or cold (4°C for 24 h) as described by Li et al. (2017b, 2018). Two biological replicates of control and cold‐treated samples were ground into a fine powder using liquid nitrogen. Approximately 1250 μl of polysome extraction buffer (200 mm Tris, pH 9.0, 200 mm KCl, 26 mm MgCl2, 25 mm EGTA, 100 μm 2‐mercaptoethanol, 50 μg ml–1 cycloheximide, 50 μg ml–1 chloramphenicol, 1% Triton X‐100, 1% Brij‐35, 1% Tween‐40, 1% IGEPAL CA‐630, 2% polyoxyethylene 10 tridecyl ether, 1% deoxycholic acid) was added to the powdered tissue (500 mg) and mixed thoroughly and incubated on ice for 10 min with occasional mixing by inverting tubes. The mixture was spun for 2 min at 16 900 (4°C). Approximately 600 μl of supernatant was layered onto the top of the sucrose density gradients (20–60% sucrose, w/v) and centrifuged for 120 min at 275 000  in a Optima LE‐80 centrifuge (Beckman, Brea, CA, USA). The optical density of the gradient fractions was measured by using a UA‐5 detector and a Gradient Fractionator (model 640; ISCO, Lincoln, NE, USA) at 254 nm. Sucrose fractions with more than two ribosomes were pooled and used for polysome RNA isolation using TRIzol reagent. Approximately 8 μg of polysome RNA was used for mRNA purification and subsequently the RNA‐seq library was constructed and sequenced as indicated above (MeRIP‐seq).

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‐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 peak calling

To identify m6A peaks, the peak calling software macs2 was used on the alignment bam files from star using the parameters as described previously (Anderson et al., 2018). 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.

m probability

Transcripts (genes) were selected from the GFF annotation only if they had at least one annotated 5′ UTR, exon, and 3′ UTR. If more than one 5′ UTR annotation was available for a transcript, the longest and the most upstream one was used, and, if more than one 3′ 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′ UTR and the most upstream end of the 3′ UTR were set as index (position) 0 for the 5′ UTR and 3′ UTR graphs, respectively. Finally, the probability of m6A for each category and graph was calculated as: where s(i) and T(i) represent the number of transcripts having an m6A site at index i and the number of transcripts having index i, respectively; the max term represents the maximum probability of m6A over the indices ranging from (inclusive) −n to n; and argmax represents the index at which the probability of m6A 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 Occupancycold – Ribo Occupancycontrol).

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 m6A 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 cold‐acclimated 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.

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.

Figure S1. m6A peak validations under cold stress.

Figure S2. m6A peaks identified in the present study show high overlap with previously published m6A peaks. In addition, thousands of new peaks identified in the present study.

Figure S3. Average relative read coverage of input mRNA seq compared to m6A IP reads across 100‐nucleotide binned transcript in (a) control and (b) cold conditions.

Figure S4. (a) Principal component analysis plot for the three replicates of cold and input control mRNA seq (labeled IN) and MeRIP libraries (labeled IP) after log transformation and normalization by deseq2. (b) Heatmap showing clustering of the three replicates of cold and control mRNA‐seq libraries between replicates rather than between conditions.

Figure S5. (a) Comparison of mRNA abundance for cold‐enriched m6A containing transcripts methylation in the 5′ UTR and start codon regions compared to those downstream. (b) Change in mRNA abundance of 5′ 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 m6A 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, m6A‐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 m6A in the mta mutant with and without cold stress.

Figure S14. Cold‐enriched m6A 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). m6A levels of CBF1 and CBF2 transcripts in Col‐0 and mta after 3 h of cold stress.

Click here for additional data file.

Table S1. List of genes with m6A peaks

Table S2. List of genes with APA between cold and control conditions

Click here for additional data file.

Table S3. Sequence of primers used in this study

Click here for additional data file.