Skip to main content

m6A RNA methylation impairs gene expression variability and reproductive thermotolerance in Arabidopsis


Heat-imposed crop failure is often attributed to reduced thermotolerance of floral tissues; however, the underlying mechanism remains unknown. Here, we demonstrate that m6A RNA methylation increases in Arabidopsis flowers and negatively regulates gene expression variability. Stochastic gene expression provides flexibility to cope with environmental stresses. We find that reduced transcriptional fluctuation is associated with compromised activation of heat-responsive genes. Moreover, disruption of an RNA demethylase AtALKBH10B leads to lower gene expression variability, suppression of heat-activated genes, and strong reduction of plant fertility. Our work proposes a novel role for RNA methylation in the bet-hedging strategy of heat stress response.


Reproductive tissues of plants are more susceptible to heat stress than in other developmental stages [1, 2]. Previous studies suggested that unfolded protein response factors play an important role in the reproductive thermotolerance [3,4,5,6]. For example, basic leucine zipper transcription factors, bZIP17, bZIP28, and bZIP60, are essential for maintaining fertility under heat stress [3,4,5]. In addition, SQUAMOSA promoter-binding protein-like (SPL) transcription factors, SPL1 and SPL12, are required for thermotolerance at the reproductive stage of Arabidopsis [7]. Moreover, various small RNAs were identified in the heat-stressed flowers of flax, soybean, and maize [8,9,10]; however, their mode of action is still elusive.

Sporadic gene expression is an important cellular feature that increases the adaptability of an organism to changing environments [11,12,13]. In unicellular organisms, noisy transcription enhances fitness in adverse growth conditions [11, 14,15,16]. Stochastic gene expression is also observed in multicellular organisms; for instance, clonal heterogeneity of gene expression is important for lineage choice of mouse hematopoietic progenitor cells [17], the intestinal differentiation of nematode is also influenced by gene expression variability [18], and a robust transcription of NF-κB is attributed to both intrinsic and extrinsic variability [19]. In Arabidopsis, it was previously suggested that genes with variable expression were enriched with factors involved in stress responses [20,21,22,23]. Therefore, noisy gene expression is an evolutionary conserved and biologically relevant mechanism that enhances flexibility of cellular responses.

In this study, we show that N6-methyladenosine (m6A), a widespread RNA modification in mRNAs [24,25,26], is strongly increased in Arabidopsis flowers and inversely correlated with gene expression variability, which then compromises reproductive thermotolerance. Our work suggests a novel role for RNA methylation and gene expression variability in heat tolerance, which can help mitigate the heat-imposed reproductive failure of crop plants.

Results and discussion

Plants’ resilience to environmental challenges varies during development and floral tissues are particularly more sensitive to heat stress [1, 2]. Such heat susceptibility during the reproductive development can cause the reduction of yield and quality of fruit and cereal crops. To understand the molecular mechanisms underlying the heat-imposed reproductive failure in Arabidopsis, we carried out transcriptomic analyses of the heat-stressed flower and leaf samples and found that flowers exhibited a distinct pattern of gene expression under heat stress (Additional file 1: Fig. S1). Genes in clusters 1 and 10, for instance, were strongly activated upon heat treatment in leaves, while in flowers, their expression levels were not altered dramatically (Additional file 1: Fig. S1). Unfortunately, the fundamental mechanisms driving such divergence of heat responsiveness in different tissues are not well understood.

Previous studies suggested that RNA modifications (often referred to as epitranscriptomic marks) are implicated in cellular response to environmental challenges [27,28,29]; however, the functional roles of RNA modification in the heat stress response of Arabidopsis remain largely unknown. In order to test if the developmentally divergent heat response can be attributed to m6A, the most widespread RNA modification type in mRNA, we carried out m6A-RNA immunoprecipitation (RIP)-seq experiments using the same samples tested in Fig. S1 (Fig. 1a and b; Additional file 1: Fig. S2 and S3). Unlike the vast changes of gene expression levels observed in Fig. S1, the m6A levels did not show any dramatic difference in the heat stress condition and between leaves and flowers (Additional file 1: Fig. S4). We then measured the m6A enrichment (determined by fold change of m6A-immunoprecipitated to total RNA) and found that the gene expression levels were poorly correlated with the m6A enrichment (Additional file 1: Fig. S5). In addition, Pearson correlation coefficient analysis revealed that heat stress-induced gene expression and m6A enrichment changes are only marginally correlated (Additional file 1: Fig. S6). This is in agreement with previous studies that suggested a negligible effect of RNA methylation in transcriptome-wide mRNA steady-state levels in Arabidopsis [30]. Taken together, heat stress causes distinct transcriptomic changes in the Arabidopsis leaves and flowers and triggers little changes of m6A levels.

Fig. 1
figure 1

RNA methylation and expression variability is inversely correlated. a Density distribution of m6A in the control and heat-stressed leaves and flowers of Arabidopsis. UTR and CDS are marked as dark and light grey boxes, respectively. CS, control sample; HS, heat-stressed sample. b Consensus sequence motif of m6A-modified sites in the non-stressed leaf sample. The top 1000 most m6A-enriched regions were used for the analysis. c Number of genes with m6A enrichment score greater than 0. Enrichment of m6A was determined by log2-converted fold change of m6A-immunoprecipitated to input RNA levels. d, e Number of HVGs (d) and LVGs (e) in leaves and flowers subjected to the control and heat stress treatments. HVGs and LVGs are genes with log2CV2 greater than 5 and lower than − 10, respectively. CV2 refers to coefficient of variation of a gene which was determined by FPKMs of RNA-seq performed for three biological replications. f Fraction of m6A-modified transcripts in quartiles of gene expression variability. Methylated transcript was as defined in c

Our m6A RNA methylation data revealed that the flower samples display distinct m6A pattern from leaves exhibiting higher levels along the CDS and lower levels around the 3′ UTR (Fig. 1a). More importantly, a remarkably stronger enrichment of m6A was detected in flowers (Fig. 1c; Additional file 1: Fig. S7), which is consistent with previous studies [31, 32]. We therefore hypothesized that strong m6A enrichment in flowers might play a role in the reproductive thermosensitivity by modulating certain RNA processes. Intriguingly, we found a strong negative correlation between m6A enrichment and variability of gene expression (Fig. 1d–f). Gene expression variability refers to stochastic and noisy expression of genetically identical cells and can be determined by the squared coefficient of variation (CV2) that is a measure of dispersion of data points. In this study, we generated RNA-seq datasets for three independent biological replicates and interestingly, genes with invariable expression in biological replications significantly overlapped with those invariable genes identified from single-cell (sc) and single-plant RNA-seq datasets (Additional file 1: Fig. S8; Additional file 2: Table S1). This suggests that biological replications can be used to infer gene expression variability.

As shown in Fig. 1d and e, while highly variable genes (HVGs) were found similar in leaves and flowers, lowly variable genes (LVGs) were in greater numbers in flowers where stronger m6A enrichment is observed. Figure 1f also shows that transcripts with low gene expression variability contain more methylated mRNAs. The negative correlation of gene expression variability and m6A levels was also observed for transcripts with higher number of m6A peaks which display low expression variability (Additional file 1: Fig. S9; Additional file 3: Table S2). We also examined the m6A levels of HVGs and LVGs determined from other datasets. In the study of Cortijo et al., HVGs and LVGs were identified from inter-individual expression variability of Arabidopsis [23], and LVGs contained more m6A-methylated transcripts (Additional file 1: Fig. S10). Similar result was also found for HVGs and LVGs identified from a scRNA-seq dataset. We analyzed a previously published scRNA-seq data generated from Arabidopsis root tip [33] and determined both HVGs and LVGs. Similarly, LVGs contained more m6A-modified transcripts than HVGs (Additional file 1: Fig. S10). It is also worth noting that the LVGs found in flowers were strongly enriched with genes involved in abiotic stress response (Additional file 1: Fig. S11), which partly indicates a functional association between expression variability and stress resistance. In conclusion, the increased RNA methylation in flowers is associated with lower gene expression variability.

In search of possible regulators of m6A in Arabidopsis flower, we focused on AtALKBH10B (hereinafter 10B) because its expression level is the highest in the floral tissue among the five RNA demethylases encoded in the Arabidopsis genome (Additional file 1: Fig. S12). To test if 10B plays a role in the m6A-associated gene expression variability, we generated the m6A-RIP-seq and RNA-seq data using the 10b-1 mutant (Additional file 1: Fig. S13). In the m6A enrichment analyses, we found only limited number of transcripts with altered m6A levels in the leaf sample of 10b-1 (Additional file 4: Table S3) and were unable to observe any significant differences in m6A distribution between the wild type (wt) and 10b-1 leaves (Additional file 1: Fig. S14). However, the flower sample of 10b-1 exhibited distinct m6A distribution as compared with wt, showing stronger peak around the stop codon and weaker enrichment along the CDS (Fig. 2a). Most importantly, the m6A enrichment level was greatly increased in the flowers of 10b-1, while in leaves the difference of RNA methylation was only marginal (Fig. 2b; Additional file 4: Table S3), which collectively suggests that 10B mediates RNA demethylation mostly in flowers of Arabidopsis. We thus assessed the gene expression variability of the 10B-regulated genes, which were defined as those with increased RNA methylation in the 10b-1 mutant. Interestingly, we found that the 10B-regulated transcripts exhibit reduced level of gene expression variability as compared with randomly selected mRNAs (Fig. 2c; Additional file 1: Fig. S15). These data further support the notion that increased m6A RNA methylation leads to reduction in gene expression variability.

Fig. 2
figure 2

AtALKBH10b is required for reproductive thermotolerance. a Density distribution of m6A in the wt and 10b-1 mutant flowers. UTR and CDS are marked as dark and light grey boxes, respectively. CS, non-stressed control sample. b Cumulative distribution of the m6A enrichment fold change in 10b-1 normalized against Col-0. P value was obtained by the one-tailed Wilcoxon rank sum test. c Cumulative distribution of expression variability difference between wt and 10b-1 determined for the 10B-regulated and randomly selected transcripts. The 10B-regulated genes are those with at least two-fold increase of m6A enrichment in 10b-1 compared to wt. P value was obtained by the one-tailed Wilcoxon rank sum test. d Heatmap of FC expression upon heat stress in different samples. Genes upregulated by the heat treatment in wt leaves by at least two-fold were selected. e Fraction of normal and short siliques of plants subjected to the control and heat stress treatments. CS, control sample; HS, heat-stressed sample. Siliques that are longer than 10 mm were considered normal and those shorter than 10 mm were counted as short siliques. f A representative image for siliques of the non-stressed and heat-stressed wt and 10b-1 plants. Bar is 10 mm. g Schematic illustration of m6A-associated gene expression variability. Bars represent the level of gene expression of individual samples in heat. Floral tissues exhibit stronger m6A RNA methylation and lower expression variability. Dashed line indicates a critical threshold required for heat tolerance. Reduced expression variability in flowers diminishes the gene activation to a level below the critical point

Previously, several studies proposed that stochastic gene transcription can be dictated by promoter structure, DNA methylation, and chromatin status [16, 23, 34]. Further to the control at DNA or chromatin levels, we showed that an epitranscriptomic mark adds to the complex regulation of gene expression variability. Given that m6A-methylated mRNAs are preferably located to cytoplasmic stress granules (SGs) [27, 35,36,37], we speculated that the physical confinement of m6A-modified transcripts to SGs might account for the invariability of gene expression. To test this possibility, we compared LVGs and random genes for their SG enrichment score determined in our previous study [38], and LVGs indeed showed stronger SG enrichment (Additional file 1: Fig. S16). However, the precise mechanisms as to how low gene expression variability is attributed to RNA sequestration to SGs require further investigation.

Lastly, in order to understand the functional relevance of gene expression variability in the heat stress response of Arabidopsis, we examined the heat responsiveness of genes in leaves and flowers of wt and 10b-1. For this, the heat-activated genes in the wt leaves were selected and their fold change expression in heat was compared. Noticeably, the fold change of heat-activated genes was drastically compromised in the wt flowers, and it was further reduced in the 10b-1 flowers (Fig. 2d). Importantly, the 10b-1 mutant plants exhibited more severe fertility defects in the heat stress condition (Fig. 2e and f), indicating that 10B is required for heat resistance of Arabidopsis flowers. In conclusion, 10B is an m6A RNA demethylase acting in the floral tissue and contributes to both variable gene expression and heat responsiveness, all of which is required for reproductive success under environmentally challenging conditions.


It has been previously suggested that stochastic gene expression is advantageous to unicellular organisms with respect to survival under threatening environmental conditions [11, 15, 20, 23]. Such behavior of variable gene expression was also observed in higher eukaryotes including mammals and was found to be associated with disease susceptibility [13, 39, 40]. Gene expression variability can be determined at different levels from cellular to population scales. Recently, Cortijo et al. investigated inter-individual gene expression variability of Arabidopsis during the diurnal circadian oscillation and discovered many HVGs involved in environmental responses [23]. In this study, we showed that Arabidopsis flowers exhibit reduced variability of gene expression, which limits the activation of heat-responsive genes to a level below certain critical point and consequently compromises fertility under heat stress (Fig. 2g). Therefore, our work suggests a novel framework that gene expression variability is an important factor impacting heat resilience of Arabidopsis and m6A RNA methylation is key to gene expression variability.

Materials and methods

Plant materials and growth condition

Arabidopsis Col-0 and atalkbh10b-1 (SALK_004215C) mutant plants were grown at 22 °C under 16 h light/8 h dark day/night cycle. Heat stress treatment was performed to 5-week-old plants by elevating the growth temperature to 37 °C for 3 h. Rosette leaves (n = 6) and stage 1 to 12 floral buds (n = 12) were collected for both RNA-seq (3 biological replicates) and m6A-RIP-seq (2 biological replicates).

NGS library construction

Total RNA was isolated using the Trizol reagent (Invitrogen). For m6A-RIP-seq, 50 μg of total RNA was used. Briefly, poly(A)-RNA was selected using the oligo-d(T)25 magnetic beads (Thermo Fisher) and was fragmented using Magnesium RNA Fragmentation Module (NEB) at 86 °C for 7 min. Cleaved RNA fragments were then incubated for 2 h at 4 °C with m6A-specific antibody (No. 202003, Synaptic Systems) in the IP buffer (50 mM Tris-HCl, 750 mM NaCl and 0.5% Igepal CA-630). Library preparation was performed using the NEBNext Ultra Directional RNA Library Prep Kit (NEB) following the manufacturer’s instructions. Finally, 150-bp paired-end sequencing (PE150) was carried out on an Illumina Novaseq 6000. For RNA-seq, 3 μg of total RNA was used, and the libraries were constructed by following the same method as described above.

NGS data analyses

Raw sequencing reads were cleaned using Trimmomatic (version 0.39) to remove reads containing adapter and low-quality sequences [41]. Trimmed reads were then aligned to the Arabidopsis reference genome (TAIR10) with default settings using Hisat2 (version 2.2.1) [42]. FPKM values were calculated by StringTie (version 2.1.7) [43]. Visualization of the sequencing data was performed using the Integrative Genomics Viewer (IGV) [44]. For gene expression variability analysis, we used the squared coefficient of variation obtained from the FPKM values of three biological replicates. Coefficient of variation is the ratio of standard deviation to mean. For m6A peak calling, MACS2 (version was run with the following parameters: --nomodel,--extsize 50, -p 5e-2, and -g 65084214 [45]. The -g option accounts for the size of the Arabidopsis transcriptome. m6A peaks identified in both two replicates were annotated by CHIPseeker [46]. Distribution of m6A peaks along mRNAs were analyzed by adopting the MeRIP-PF and Guitar scripts [47, 48]. Consensus sequence motif was analyzed using the top 1000 most significant m6A-enriched regions. RRACH sequence motif was searched in the selected regions, and its frequency was calculated and visualized using the weblogo tool [49]. The NGS datasets generated in this study is summarized in Additional file 5: Table S4.


m6A-RIP-qPCR experiment was performed as described previously with minor modifications [50] and the Magna RIP™ RNA-Binding Protein Immunoprecipitation Kit (Merck) was used following the manufacturer’s instruction. Briefly, 300 μg of total RNA was fragmented using RNA fragmentation buffer [for 1 mL 10X reagents: 800 μL 1 M Tris-HCl (pH 7.0), 100 μL 1 M ZnCl2, 100 μL RNase-free H2O] to a size of about 250 nucleotides. Fragmented RNA was then precipitated using 2.5 volume of ethanol, 1/10 volume of 3 M NaOH, and 100 μg/mL glycogen at − 80 °C overnight. Precipitated RNA was resuspended in 55 μL RNase-free H2O, 5 μL of which was kept as input sample. The remaining RNA was incubated with 5 μg of m6A-specific antibody (No. 202003, Synaptic Systems) overnight at 4 °C. Immunoprecipitated RNA was isolated using magnetic beads and the beads were washed five times with 500 μL of cold RIP wash buffer. After resuspending the beads in 150 μL of proteinase K buffer (117 μL of RIP wash buffer, 15 μL of 10% SDS, 18 μL of 10 mg/mL proteinase K) and incubation at 55 °C for 30 min, both input and immunoprecipitated RNA was isolated by the phenol:chloroform method and resuspended in 20 μL of RNase-free H2O. Reverse transcription reaction was performed using the ReverTra Ace qPCR RT Master Mix (Toyobo). The resulting first-strand cDNA was diluted four-fold with DEPC-treated water and 1.5 μL was used for a 20-μL qPCR reaction. Real-time qPCR was carried out using ChamQ Universal SYBR qPCR Master Mix (Vazyme) in the CFX96 Connect Real-time PCR Detection system (BioRad). Actin2 was used as an internal control and the sequences of the primers are provided in Additional file 6: Table S5.

Heat resistance phenotyping

Heat tolerance test at reproductive stage was carried out as described previously [4]. Briefly, plants were grown under normal growth condition set at 22 °C and 16 h light/8 h dark day/night cycle. Heat stress of 37 °C was treated for 6 h when the first flower opened. After the heat treatment, plants were moved to normal growth condition and continued to grow for 10 days. Silique length was measured from eight individual plants per genotype.

Availability of data and materials

The datasets supporting the conclusions of this article are available in the SRA repository, PRJNA793364 ( [51] and summarized in Additional file 5: Table S4. The analyses were performed using the standard codes instructed by the tools described in the “Methods” section, and the custom codes used in this study are available under MIT license at GitHub ( and Zenodo ( [52].


  1. Jacott CN, Boden SA. Feeling the heat: developmental and molecular responses of wheat and barley to high ambient temperatures. J Exp Bot. 2020;71:5740–51.

    CAS  PubMed  PubMed Central  Google Scholar 

  2. Chaturvedi P, Wiese AJ, Ghatak A, Záveská Drábková L, Weckwerth W, Honys D. Heat stress response mechanisms in pollen development. New Phytol. 2021;231:571–85.

    PubMed  PubMed Central  Google Scholar 

  3. Gao J, Wang M-J, Wang J-J, Lu H-P, Liu J-X. bZIP17 regulates heat stress tolerance at reproductive stage in Arabidopsis. aBIOTECH. 2022;3:1–11.

    CAS  PubMed  Google Scholar 

  4. Zhang S-S, Yang H, Ding L, Song Z-T, Ma H, Chang F, et al. Tissue-specific transcriptomics reveals an important role of the unfolded protein response in maintaining fertility upon heat stress in Arabidopsis. Plant Cell. 2017;29:1007–23.

    CAS  PubMed  PubMed Central  Google Scholar 

  5. Deng Y, Humbert S, Liu J-X, Srivastava R, Rothstein SJ, Howell SH. Heat induces the splicing by IRE1 of a mRNA encoding a transcription factor involved in the unfolded protein response in Arabidopsis. Proc Natl Acad Sci. 2011;108:7247–52.

    CAS  PubMed  PubMed Central  Google Scholar 

  6. Iwata Y, Koizumi N. An Arabidopsis transcription factor, AtbZIP60, regulates the endoplasmic reticulum stress response in a manner unique to plants. Proc Natl Acad Sci. 2005;102:5280–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  7. Chao L-M, Liu Y-Q, Chen D-Y, Xue X-Y, Mao Y-B, Chen X-Y. Arabidopsis transcription factors SPL1 and SPL12 confer plant thermotolerance at reproductive stage. Mol Plant. 2017;10:735–48.

    CAS  PubMed  Google Scholar 

  8. Pokhrel S, Meyers BC. Heat-responsive microRNAs and phased small interfering RNAs in reproductive development of flax. Plant Direct. 2022;6:e385.

    CAS  PubMed  PubMed Central  Google Scholar 

  9. He J, Jiang Z, Gao L, You C, Ma X, Wang X, et al. Genome-wide transcript and small RNA profiling reveals transcriptomic responses to heat stress. Plant Physiol. 2019;181:609–29.

    CAS  PubMed  PubMed Central  Google Scholar 

  10. Ding X, Guo J, Zhang Q, Yu L, Zhao T, Yang S. Heat-responsive miRNAs participate in the regulation of male fertility stability in soybean CMS-based F1 under high temperature stress. Int J Mol Sci. 2021;22:2446.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. Grimbergen AJ, Siebring J, Solopova A, Kuipers OP. Microbial bet-hedging: the power of being different. Curr Opin Microbiol. 2015;25:67–72.

    PubMed  Google Scholar 

  12. Eldar A, Elowitz MB. Functional roles for noise in genetic circuits. Nature. 2010;467:167–73.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. de Jong TV, Moshkin YM, Guryev V. Gene expression variability: the other dimension in transcriptome analysis. Physiol Genomics. 2019;51:145–58.

    PubMed  Google Scholar 

  14. Liu J, Martin-Yken H, Bigey F, Dequin S, François J-M, Capp J-P. Natural yeast promoter variants reveal epistasis in the generation of transcriptional-mediated noise and its potential benefit in stressful conditions. Genome Biol Evol. 2015;7:969–84.

    CAS  PubMed  PubMed Central  Google Scholar 

  15. Elowitz MB, Levine AJ, Siggia ED, Swain PS. Stochastic gene expression in a single cell. Science (80- ). 2002;297:1183–6.

    CAS  Google Scholar 

  16. Jones DL, Brewster RC, Phillips R. Promoter architecture dictates cell-to-cell variability in gene expression. Science (80- ). 2014;346:1533–6.

    CAS  Google Scholar 

  17. Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S. Transcriptome-wide noise controls lineage choice in mammalian progenitor cells. Nature. 2008;453:544–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. Raj A, Rifkin SA, Andersen E, van Oudenaarden A. Variability in gene expression underlies incomplete penetrance. Nature. 2010;463:913–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Kellogg RA, Tay S. Noise facilitates transcriptional control under dynamic inputs. Cell. 2015;160:381–92.

    CAS  PubMed  Google Scholar 

  20. Hirao K, Nagano AJ, Awazu A. Noise–plasticity correlations of gene expression in the multicellular organism Arabidopsis thaliana. J Theor Biol. 2015;387:13–22.

    PubMed  Google Scholar 

  21. Araújo IS, Pietsch JM, Keizer EM, Greese B, Balkunde R, Fleck C, et al. Stochastic gene expression in Arabidopsis thaliana. Nat Commun. 2017;8:2132.

    PubMed  PubMed Central  Google Scholar 

  22. Bhosale R, Jewell JB, Hollunder J, Koo AJK, Vuylsteke M, Michoel T, et al. Predicting gene function from uncontrolled expression variation among individual wild-type Arabidopsis plants. Plant Cell. 2013;25:2865–77.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Cortijo S, Aydin Z, Ahnert S, Locke JCW. Widespread inter-individual gene expression variability in Arabidopsis thaliana. Mol Syst Biol. 2019;15:e8591.

    PubMed  PubMed Central  Google Scholar 

  24. Meyer KD, Jaffrey SR. Rethinking m 6 a readers, writers, and erasers. Annu Rev Cell Dev Biol. 2017;33:319–42.

    CAS  PubMed  PubMed Central  Google Scholar 

  25. Yue H, Nie X, Yan Z, Weining S. N6-methyladenosine regulatory machinery in plants: composition, function and evolution. Plant Biotechnol J. 2019;17:1194–208.

    PubMed  PubMed Central  Google Scholar 

  26. Jia G, Fu Y, He C. Reversible RNA adenosine methylation in biological regulation. Trends Genet. 2013;29:108–15.

    CAS  PubMed  Google Scholar 

  27. Scutenaire J, Deragon JM, Jean V, Benhamed M, Raynaud C, Favory JJ, et al. The YTH domain protein ECT2 is an m6A reader required for normal trichome branching in arabidopsis. Plant Cell. 2018;30:986–1005.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. Meyer KD, Patil DP, Zhou J, Zinoviev A, Skabkin MA, Elemento O, et al. 5′ UTR m6A promotes cap-independent translation. Cell. 2015;163:999–1010.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian SB. Dynamic m6 a mRNA methylation directs translational control of heat shock response. Nature. 2015;526:591–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. Anderson SJ, Kramer MC, Gosai SJ, Yu X, Vandivier LE, Nelson ADL, et al. N6-methyladenosine inhibits local ribonucleolytic cleavage to stabilize mRNAs in Arabidopsis. Cell Rep. 2018;25:1146–1157.e3.

    CAS  PubMed  Google Scholar 

  31. Wan Y, Tang K, Zhang D, Xie S, Zhu X, Wang Z, et al. Transcriptome-wide high-throughput deep m6A-seq reveals unique differential m6A methylation patterns between three organs in Arabidopsis thaliana. Genome Biol. 2015;16:272.

    PubMed  PubMed Central  Google Scholar 

  32. Wang C, Yang J, Song P, Zhang W, Lu Q, Yu Q, et al. FIONA1 is an RNA N6-methyladenosine methyltransferase affecting Arabidopsis photomorphogenesis and flowering. Genome Biol. 2022;23:40.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. Zhang T-Q, Xu Z-G, Shang G-D, Wang J-W. A single-cell RNA sequencing profiles the developmental landscape of Arabidopsis root. Mol Plant. 2019;12:648–60.

    CAS  PubMed  Google Scholar 

  34. Bashkeel N, Perkins TJ, Kærn M, Lee JM. Human gene expression variability and its dependence on methylation and aging. BMC Genomics. 2019;20:941.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. Alvarado-Marchena L, Marquez-Molins J, Martinez-Perez M, Aparicio F, Pallás V. Mapping of functional subdomains in the atALKBH9B m6A-demethylase required for its binding to the viral RNA and to the coat protein of alfalfa mosaic virus. Front Plant Sci. 2021;12:701683.

  36. Fu Y, Zhuang X. m6A-binding YTHDF proteins promote stress granule formation. Nat Chem Biol. 2020;16:955–63.

    CAS  PubMed  PubMed Central  Google Scholar 

  37. Ries RJ, Zaccara S, Klein P, Olarerin-George A, Namkoong S, Pickering BF, et al. m6A enhances the phase separation potential of mRNA. Nature. 2019;571:424–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Kim EY, Wang L, Lei Z, Li H, Fan W, Cho J. Ribosome stalling and SGS3 phase separation prime the epigenetic silencing of transposons. Nat Plants. 2021;7:303–9.

    CAS  PubMed  Google Scholar 

  39. Green DJ, Sallah SR, Ellingford JM, Lovell SC, Sergouniotis PI. Variability in gene expression is associated with incomplete penetrance in inherited eye disorders. Genes (Basel). 2020;11:179.

    CAS  PubMed  Google Scholar 

  40. Hagai T, Chen X, Miragaia RJ, Rostom R, Gomes T, Kunowska N, et al. Gene expression variability across cells and species shapes innate immunity. Nature. 2018;563:197–202.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    CAS  PubMed  PubMed Central  Google Scholar 

  42. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    CAS  PubMed  PubMed Central  Google Scholar 

  43. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33:290–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.

    PubMed  PubMed Central  Google Scholar 

  46. Yu G, Wang L-G, He Q-Y. ChIPseeker: an R/bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015;31:2382–3.

    CAS  PubMed  Google Scholar 

  47. Li Y, Song S, Li C, Yu J. MeRIP-PF: an easy-to-use pipeline for high-resolution peak-finding in MeRIP-Seq data. Genom Proteom Bioinform. 2013;11:72–5.

    Google Scholar 

  48. Cui X, Wei Z, Zhang L, Liu H, Sun L, Zhang S-W, et al. Guitar: an R/bioconductor package for gene annotation guided transcriptomic analysis of RNA-related genomic features. Biomed Res Int. 2016;2016:1–8.

    Google Scholar 

  49. Crooks GE, Hon G, Chandonia J-M, Brenner SE. WebLogo: a sequence logo generator. Genome Res. 2004;14:1188–90.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. Dominissini D, Moshitch-Moshkovitz S, Salmon-Divon M, Amariglio N, Rechavi G. Transcriptome-wide mapping of N6-methyladenosine by m6A-seq based on immunocapturing and massively parallel sequencing. Nat Protoc. 2013;8:176–89.

    CAS  PubMed  Google Scholar 

  51. Wang L, Zhuang H, Fan W, Zhang X, Dong H, Yang H, et al. Arabidopsis thaliana heat stress epitranscriptome. Datasets Gene Expression Omnibus. 2022;

  52. Wang L, Zhuang H, Fan W, Zhang X, Dong H, Yang H, et al. Custom codes for “m6A RNA methylation impairs gene expression variability and reproductive thermotolerance in Arabidopsis”. Zenodo. 2022.

Download references


Not applicable.

Peer review information

Wenjing She was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Review history

The review history is available as Additional file 7.


This work was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (XDB27030209), National Natural Science Foundation of China (31970518, 32150610473 and 32111540256), and General Program of Natural Science Foundation of Shanghai (22ZR1469100) granted to JC. HY is supported by the Shanghai Landscaping Administrative Bureau Program (G172404) and Shanghai Pujiang Program (16PJ1403000).

Author information

Authors and Affiliations



LW, HZ, WF, XZ, and HD conducted the experiments. LW, HY, and JC analyzed the data. JC was a major contributor in writing the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Hongxing Yang or Jungnam Cho.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information


Additional file 1: Fig. S1. Divergent transcriptomic changes of leaves and flowers in heat. Fig. S2. Comparison of m6A enrichment in two biological replications. Fig. S3. Detected sequence motifs of m6A-modified sites. Fig. S4. m6A peaks detected in different clusters. Fig. S5. Correlation between m6A and expression levels. Fig. S6. Correlation between FC m6A and FC expression. Fig. S7. Increased m6A RNA modification in Arabidopsis flowers. Fig. S8. Overlap of lowly variable genes identified from different datasets. Fig. S9. Inverse correlation of m6A levels and gene expression variability. Fig. S10. Lowly variable genes show stronger m6A RNA modification. Fig. S11. Enrichment of genes associated with low expression variability. Fig. S12. Expression profiling of genes encoding RNA demethylases. Fig. S13. An exemplary locus for increased m6A in 10b-1. Fig. S14. Density distribution of m6A enrichment in Arabidopsis leaves. Fig. S15. Exemplary loci with decreased expression variability in 10b-1. Fig. S16. Stress-granule association of lowly variable transcripts.

Additional file 2: Table S1. LVGs identified by different datasets.

Additional file 3: Table S2. Expression variability and m6A peak numbers.

Additional file 4: Table S3. Hypermethylated genes in leaves and flowers of 10b-1.

Additional file 5: Table S4. Summary of sequencing data generated in this study.

Additional file 6: Table S5. Sequences of oligonucleotides.

Additional file 7. Review history.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, L., Zhuang, H., Fan, W. et al. m6A RNA methylation impairs gene expression variability and reproductive thermotolerance in Arabidopsis. Genome Biol 23, 244 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: