Baoqi Li,Mengmeng Zhang,Weinan Sun,Dandan Yue,Yizan Ma,Boyang Zhang,Lingfeng Duan,Maojun Wang,Keith Lindsey,Xinhui Nie,Xianlong Zhang,Xiyan Yang
N6 ‐methyladenosine RNA modification regulates cotton drought response in a Ca2+ and ABA ‐dependent manner Baoqi Li et al.
Cotton,Ca2+ signalling pathway,Drought stress,GhALKBH10B,mRNA stability,N6‐methyladenosine
Summary N6‐methyladenosine (m6A) is the most prevalent internal modification present in mRNAs, and is considered to participate in a range of developmental and biological processes. Drought response is highly regulated at the genomic, transcriptional and post‐transcriptional levels. However, the biological function and regulatory mechanism of m6A modification in the drought stress response is still poorly understood. We generated a transcriptome‐wide m6A map using drought‐resistant and drought‐sensitive varieties of cotton under different water deficient conditions to uncover patterns of m6A methylation in cotton response to drought stress. The results reveal that m6A represents a common modification and exhibit dramatic changes in distribution during drought stress. More 5'UTR m6A was deposited in the drought‐resistant variety and was associated with a positive effect on drought resistance by regulating mRNA abundance. Interestingly, we observed that increased m6A abundance was associated with increased mRNA abundance under drought, contributing to drought resistance, and vice versa. The demethylase GhALKBH10B was found to decrease m6A levels, facilitating the mRNA decay of ABA signal‐related genes (GhZEP, GhNCED4 and GhPP2CA) and Ca2+ signal‐related genes (GhECA1, GhCNGC4, GhANN1 and GhCML13), and mutation of GhALKBH10B enhanced drought resistance at seedling stage in cotton. Virus‐induced gene silencing (VIGS) of two Ca2+‐related genes, GhECA1 and GhCNGC4, reduced drought resistance with the decreased m6A enrichment on silenced genes in cotton. Collectively, we reveal a novel mechanism of post‐transcriptional modification involved in affecting drought response in cotton, by mediating m6A methylation on targeted transcripts in the ABA and Ca2+ signalling transduction pathways.
N6 ‐methyladenosine RNA modification regulates cotton drought response in a Ca2+ and ABA ‐dependent manner Baoqi Li et al.

N6‐methyladenosine (m6A) is the most prevalent internal modification present in mRNAs, and is considered to participate in a range of developmental and biological processes. Drought response is highly regulated at the genomic, transcriptional and post‐transcriptional levels. However, the biological function and regulatory mechanism of m6A modification in the drought stress response is still poorly understood. We generated a transcriptome‐wide m6A map using drought‐resistant and drought‐sensitive varieties of cotton under different water deficient conditions to uncover patterns of m6A methylation in cotton response to drought stress. The results reveal that m6A represents a common modification and exhibit dramatic changes in distribution during drought stress. More 5'UTR m6A was deposited in the drought‐resistant variety and was associated with a positive effect on drought resistance by regulating mRNA abundance. Interestingly, we observed that increased m6A abundance was associated with increased mRNA abundance under drought, contributing to drought resistance, and vice versa. The demethylase GhALKBH10B was found to decrease m6A levels, facilitating the mRNA decay of ABA signal‐related genes (GhZEP, GhNCED4 and GhPP2CA) and Ca2+ signal‐related genes (GhECA1, GhCNGC4, GhANN1 and GhCML13), and mutation of GhALKBH10B enhanced drought resistance at seedling stage in cotton. Virus‐induced gene silencing (VIGS) of two Ca2+‐related genes, GhECA1 and GhCNGC4, reduced drought resistance with the decreased m6A enrichment on silenced genes in cotton. Collectively, we reveal a novel mechanism of post‐transcriptional modification involved in affecting drought response in cotton, by mediating m6A methylation on targeted transcripts in the ABA and Ca2+ signalling transduction pathways.

N6‐methyladenine (m6A) is a methylation modification occurring at the 6th N atom of adenine in RNA, and is widely distributed, in mammals, plants, fruit flies, yeasts, and viruses (Cantara et al., 2011; Wei et al., 2017). The levels of m6A mRNA methylation are promoted by methyltransferase complexes, named as “writers” (Liu et al., 2014b; Ping et al., 2014; Wang et al., 2016b), and removed by m6A demethylases, named as “erasers” (Jia et al., 2011; Zheng et al., 2013). The m6A‐modified transcripts are recognized by the “reader” proteins, which mediate the downstream effects of the m6A modification (Wang et al., 2015; Xu et al., 2014). In mammals, METTL3/METTL14 act as core methyltransferase to form a complex with other subunits (such as WTAP and VIRMA, etc.), which synergizes with m6A methylation modification in mature mRNA (Liu et al., 2014b; Ping et al., 2014; Yue et al., 2018). FTO and ALKBH5 are responsible for removing m6A modifications on adenine (Zheng et al., 2013). Plant m6A methyltransferase complex components such as MTA, MTB, FIP37 (homologous to mammalian METTL3, METTL14, and WTAP), VIR, HAKAI (Ping et al., 2014; Růžička et al., 2017; Zhong et al., 2008), demethylase ALKBH9B and ALKBH10B (Martinez‐Perez et al., 2017; Zheng et al., 2013), and reading proteins such as ECT2, ECT3, ECT4 and CPSF30‐L (homologous to the YTH family in mammals) (Arribas‐Hernandez et al., 2018; Scutenaire et al., 2018; Song et al., 2021; Wei et al., 2018) have also been identified and characterized. Through the interaction between methyltransferases and demethylases, m6A is dynamically balanced in organisms, resulting in docking sites for m6A binding proteins and appropriate assembly of RNA secondary structures (Liu et al., 2015). Most studies have shown that m6A is enriched near the 3'UTR in mammals and plants, and reading proteins participate in the regulation of mRNA stability or translation process by targeting this region (Luo et al., 2020; Miao et al., 2022; Zhou et al., 2022). Moreover, the m6A methylation occurs in conserved sequence context as RRACH and plant‐specific motif URUAY, with some specific sequence preferences in different plants or different tissues (Duan et al., 2017; Luo et al., 2020; Wei et al., 2018; Zhou et al., 2022). As a dynamic and reversible post‐transcriptional modification, the m6A modification affects almost all aspects of RNA processing, including RNA splicing, RNA stability, nuclear retention, and translation efficiency (Anderson et al., 2018; Dominissini et al., 2012; Fu et al., 2014; Liang et al., 2020; Meyer et al., 2012; Wang et al., 2014; Zhao et al., 2014, 2017; Zhou et al., 2015), thus participating in different developmental and biological processes, including cell fate determination, trichome branching, flowering transition, stress responses, spore development, and fruit maturation in plants (Arribas‐Hernandez et al., 2018; Hu et al., 2021; Scutenaire et al., 2018; Song et al., 2021; Zhang et al., 2019; Zhou et al., 2019, 2021). Drought is one of the most important environmental factors affecting crop distribution, plant development and yield production worldwide. During the long evolutionary process, plants have developed many response mechanisms to drought stress at physiological, biochemical, cellular and molecular levels, including the regulation of osmotic regulatory substances, reactive oxygen defence systems, stomatal conductance, as well as the activation of signal transduction pathways, such as Ca2+ signalling and ABA signalling pathways, transcription factor activation and genes responsible for controlling ion homeostasis (Dittrich et al., 2019; Gupta et al., 2020; Rodrigues et al., 2019; Saxena et al., 2016; Su et al., 2013; Tardieu et al., 2018). Consequently, m6A‐mediated abiotic stress responses were studied in maize (Miao et al., 2020), wheat (Sun et al., 2020), pak‐choi (Liu et al., 2020), sweet sorghum (Zheng et al., 2021), rice (Cheng et al., 2021), tomato (Yang et al., 2021), sea buckthorn (Zhang et al., 2021), cotton (Wang et al., 2022), apple (Hou et al., 2022) and barley (Su et al., 2022). These results reveal that m6A modification participates in controlling the responses of crops to abiotic stress, including drought stress. However, the characteristics and regulatory mechanisms of m6A in cotton response to drought stress have not been defined. Cotton (Gossypium hirsutum L.) is an important worldwide commercial crop grown as a source of natural fibre and edible oil. In a previous study, we characterized two upland cotton varieties which were identified as drought‐sensitive accession ZY7 and drought‐resistant accession ZY168, using an automatic phenotyping platform to acquire several i‐traits used as indicators for drought response in cotton (Li et al., 2020). Here, we describe transcriptome‐wide m6A methylation identification, and investigate the dynamic regulation of m6A in cotton under different water deficient conditions using these two varieties, and determine its biological significance by disturbing the expression of the m6A eraser ALKBH10. The results reveal that m6A represents a prevalent modification and exhibit dramatically dynamic distribution during drought stress. We show that decreased m6A modification by the demethylase GhALKBH10B reduces mRNA stability of genes in the Ca2+ and ABA signalling pathways, thereby facilitating the drought response of cotton. Moreover, virus‐induced gene silencing (VIGS)‐mediated knockdown of genes related to the Ca2+ signalling pathway reduced drought tolerance through effects on ion exchange. Our study uncovers the regulatory effects of m6A methylation on drought response and reveals a direct role for m6A methylation in the regulation of key elements in the Ca2+ signalling pathway.

To investigate whether m6A methylation participates in modulating cotton drought response, liquid chromatography tandem mass spectrometry (LC–MS/MS) was used to characterize the m6A levels during dynamic drought stress using two upland cotton varieties (Gossypium hirsutum), which were identified as drought‐sensitive (accession ZY7) and drought‐resistant (accession ZY168). The drought treatments were conducted as in our previous study (Li et al., 2020) (Figure 1a). RNA was isolated from leaf of the two varieties at three time points i.e. mild drought (MD), severe drought (SD), and re‐watering after drought (RD). These samples were named as MT (samples from mild drought treatment), ST (samples from severe drought treatment) and RT (samples from re‐watering after drought) and their parallel controls (MC, SC and RC) were used to conduct the experiments. LC–MS/MS analysis revealed a dramatic change of m6A distribution under different water deficient conditions, whereby the methylation levels of N6‐methyladenine (m6A/A) in mRNA of the two cotton varieties at the seedling stage ranged from 0.36% to 0.59% (Figure 1b). The overall m6A/A showed no significant difference between treatment and control under MD, but significant increasing was observed under SD in both cotton varieties, with the ratio of m6A/A increasing from 0.44% to 0.57% in ZY7 and from 0.51% to 0.59% in ZY168. Meanwhile, the ratios of m6A/A were higher in drought‐resistant ZY168 than drought‐sensitive ZY7 at MD under treatments and controls. However, the value of m6A/A decreased in droughted plants and showed no significant difference compared with the control after relieving drought stress. The results show that drought stress basically promotes the m6A content in cotton, however with significant differences in drought‐sensitive and drought‐resistant varieties.

m6A‐immunoprecipitation (IP) and matched input (non‐IP control) libraries of ZY7 and ZY168 were constructed and sequenced under MD, SD and RD conditions and subjected to massively parallel sequencing (Dominissini et al., 2013). The average data size of valid bases in MeRIP‐seq was 9.59G and 6.05G for input and IP samples, and the mean transcript coverage of input and IP was 41.25X and 26.08X in clean data (Table S1). The correlation analyses showed a good duplication between biological replicates that R2 ranged from 0.947 to 0.994 in input samples, and 0.965 to 0.990 in IP samples (Figures S1 and S2). These results suggest the dynamic change of mRNA transcription under different drought treatments as well as good consistency between the biological replicates. Peak numbers from MeRIP data reached saturation point when the input and m6A‐IP of mRNA reads reached a plateau of 40 million for each sample (Figure S3). The high‐confident consistent peaks (identified in both biological replicates) were used for the following analysis (Table S2). The results of peak calling showed that total number of consistent m6A peak could be identified ranging from 20 755 to 25 731 between the two cotton varieties under different watering regimes. The m6A profile of all peaks identified in two varieties under different drought condition in treatment with its parallel control panel was shown in Figure 1c. Comparing the conservative m6A peaks in the two cotton varieties, it was found that the number of peaks decreased with severe water stress (Figure S4). The differentially methylated peaks (with differentially m6A methylated enrichment under Treatment/Control, FDR <0.05) and the related differentially methylated transcripts (DMeTs) were obtained for further analysis (Table S2). 1266 and 1202 differentially methylated peaks were presented in ZY7 under MD and SD respectively, and 1841 and 1294 differentially methylated peaks under MD and SD respectively in ZY168 (Figure 2a). After RD, the DMeTs decreased to 732 and 526 in ZY7 and ZY168, respectively. Meanwhile, a total of 4420 non‐redundant DMeTs with average peaks (ratio of peak and transcript number) ranged from 1.04 to 1.12 across different conditions in these two cotton varieties (Figure 2a). Next, HOMER was applied to analyse the motif enrichment sequence within m6A peaks in cotton. Notably, the widely‐conserved motif RRACH (R = A/G; H = A/C/U) and plant‐specific motif UGUA found in Arabidopsis (Wei et al., 2018), tomato (Zhou et al., 2019), maize (Luo et al., 2020) and apple (Hou et al., 2022) was identified in cotton (Figure 2b). Moreover, we also newly identified DGCAG (D = A/G/U) and 5' deposited motif CAAUG motif in cotton (Figure 2b; Table S3). We then evaluated the distribution of m6A peaks in the whole transcriptome of cotton during dynamic drought treatments. m6A peaks were largely enriched in 3'UTRs, which is consistent with previous reports in other plants and mammals (Zhou et al., 2022). For drought‐sensitive variety ZY7, the proportion of m6A peaks located in different regions showed little changes (<1%) in the percentage ratio in treatments compared with controls from MD to RD. Surprisingly, an obvious apiciform was found in drought‐resistant variety ZY168 under both MD and SD with m6A peaks enriched near the 5'UTR (in the 5'UTR or the 1st exon) in drought treatments compared with the corresponding controls, with the proportions increased from 6.4% to 8.8% under MD, and from 6.4% to 9.5% under SD. However, the differences in the proportions were not obvious between the droughted and control plants under RD. We also identified more 5'UTR deposited m6A peaks and m6A modified transcripts in drought‐resistant variety ZY168 (762 peaks on 594 transcripts) than in drought‐sensitive variety ZY7 (310 peaks on 288 transcripts) (Figure 2c, d; Tables S4 and S5). These data not only suggest that drought has a critical influence in m6A peak distribution, but also demonstrates the dynamics and complexity of m6A modification in cotton.

To explore whether m6A methylation influences the abundance of DMeTs under drought conditions, the significant differentially m6A methylated peaks (FDR <0.05, m6A foldchange>1.5) within DMeTs were used for further analysis (Tables S4 and S5), following the reported criterion (Zhou et al., 2021). As a result, the changes of significant DMeTs abundance displayed a positive correlation with the fold change of m6A enrichment for the drought‐sensitive variety (ZY7) under MD and both varieties under SD, but not under RD, which suggested m6A methylation might affect the abundance of significant DMeTs under drought stress conditions (Figure S5). As shown in Figure 3a, 667, 680 and 202 downregulated m6A peaks and 292, 151 and 249 upregulated peaks were identified under MD, SD and RD in ZY7 (Table S4). Significant differences could be seen when comparing downregulated and upregulated group under MD (P = 1.02e‐09), SD (P = 1.13e‐05), and RD (P = 7.04e‐05). We did the same analysis in ZY168, and 228, 334 and 35 downregulated m6A peaks were identified under MD, SD and RD, with 957, 669 and 267 upregulated peaks (Table S5). Highly significant differences (P < 0.01) were only found between downregulated and upregulated group under SD (P = 2.2e‐16), but not under MD (P = 0.02) and RD (P = 0.11). We analysed the relationship between significant DMeTs abundance and m6A enrichment regions of the transcripts (Tables S4 and S5). The percentage of peaks enriched in the 5'UTR increased to 33.8% and 27.5% under MD and SD in ZY168, compared to 10.3% and 11.3% in ZY7. Significant differences (P < 0.05) were found between significant DMeTs levels of the group with m6A enrichment at the 3'UTR and 5'UTR regions under droughted conditions. While, the m6A‐deposited regions also affected the mRNA abundance of DMeTs, for that the significant DMeTs with exon‐deposited m6A methylation exhibited relatively lower down‐regulation than 5'UTR‐deposited m6A methylation in ZY168 under MD and in ZY7 under SD, and no significance between different regions under RD ‐ this was also seen in the cumulative curve results (Figure 3b; Figure S6). Interestingly, under SD, peaks enriched at the 5'UTR (27.5%) exhibited higher mRNA abundance of significant DMeTs compared to 3'UTR in ZY168, but lower abundance compared to 3'UTR enrichment in ZY7. We would propose that there is differential m6A deposition at the 5'UTR in drought‐resistance ZY168 and drought‐sensitive ZY7 under drought conditions, with a positive effect on drought resistance through the effects on abundance of significant DMeTs. In these significant DMeTs, no single transcript could be identified that was simultaneously modified with m6A at the 5'UTR, exon and 3'UTR. Most are modified with just one peak located in the 3'UTR, though a few were enriched with two peaks in the same region or two different regions (Figure 3c, d). The number of significant DMeTs with two peaks occurred under MD and SD rather than RD, and in ZY168 under RD that no two peaks could be found. It therefore suggested that drought affects the number of m6A peaks in cotton.

We further screened the unique and common significant DMeTs with differentially methylated peaks in the two cotton varieties under different drought conditions. Firstly, ZY7 (677) and ZY168 (651) shared larger numbers of significant DMeTs under MD, and smaller numbers under RD (Figure 4a). While, 99, 70 and 235 common significant DMeTs were identified in both varieties under MD, SD and RD, which were the basic response of cotton plants to drought stress, with 837, 739 and 405 unique target transcripts identified in ZY7, and 1029, 894 and 265 unique target transcripts identified ZY168. And for the same variety, 81 and 256 shared significant DMeTs was identified from MD to SD (under both MD and SD) and decreased to 31 and 62 after re‐watering (under both SD and RD) in ZY7 and ZY168, respectively (Figure S7). The differentially expressed significant DMeTs (DMeETs) were obtained with the fold change (FC) in both m6A enrichment (m6A FDR <0.05, FC >1.5) and mRNA expression level (mRNA FDR <0.05, FC >2) for further analysis. It was found that only 8 DMeETs were enriched for m6A modification and mRNA abundance in both ZY7 and ZY168, and 94 and 114 unique DMeETs were identified under MD for the representative varieties (Figure 4b). More DMeETs (20) were found in both varieties under SD, as well as 204 unique DMeETs identified in ZY7 and 299 in ZY168 (Figure 4c). As we described above, peaks decreased with re‐watering after drought, and only one DMeETs (Ghir_D08G013670.2, GhSAMS2) was identified in both (Figure 4d), and the number of unique DMeETs decreased to 36 in ZY7 and 0 in ZY168. This suggests that homeostasis of m6A modification on the transcript was not severely disrupted by drought stress. To further investigate the biological significance of DMeETs, gene ontology (GO) analysis was performed to uncover the biological functions of dynamic m6A regulation, except for ZY168 under RD due to insufficient DMeETs detected (Tables S6–S10), and we focused on biological processes that involved the methylated RNAs. DMeETs were annotated with GO terms in multiple biological processes, particularly in translation, RNA methylation, gene expression, and amide biosynthesis in both ZY7 (Figure 4e; Tables S6–S8) and ZY168 (Figure S8; Tables S9 and S10). These results suggest a role for m6A methylation in the biological response to drought stress by the regulation gene expression and translation through organonitrogen and ion transport pathways. Encyclopedia of Genes and Genomes (KEGG) analysis of these DMeETs was also performed. The most significant enriched pathways for these transcripts under drought were photosynthesis, carbon metabolism, porphyrin and chlorophyll metabolism, and pentose phosphate pathway, while cysteine and methionine metabolism, and amino acid biosynthesis were enriched after re‐watering in both varieties (Figure S9; Tables S11 and S12). This indicates m6A modification is widely involved in and related to photosynthesis and carbon metabolism pathways in cotton.

The plant ALKBH protein family are demethylases responsible for RNA demethylation, and ALKBH9 and ALKBH10 are homologous to the mammalian ALKB homologous 5 (ALKBH5). In the current study, we identified 22 ALKBH homologous genes in the allotetraploid upland cotton genome, including 4 GhALKBH9 and 4 GhALKBH10 gene members (Figure S10; Table S13). The homologous gene pair GhALKBH10B‐D (Ghir_D08G007610) and GhALKBH10B‐A (Ghir_A08G007520) was found to be predominantly expressed in vegetative organs, especially in root, stem and leaf. Another pair, GhALKBH10A‐D (Ghir_D11G012740) and GhALKBH10A‐A (Ghir_A11G012800), as well as homologous genes of ALKBH9, exhibit no obvious tissue specificity and are expressed at relatively low levels (Figure 5a). We further examined drought expression patterns of these ALKBH genes. GhALKBH10B‐D and GhALKBH10B‐A displayed significant expression changes during dynamic drought stress conditions, with up‐regulation under drought, especially under SD, and down‐regulation after re‐watering (Figure 5b). GhALKBH10A‐D and GhALKBH10A‐A however maintained low expression levels from the beginning of drought application to re‐watering in both accessions (Figure S11). CRISPR‐Cas9 gene editing technology was used to create the alkbh10a, alkbh10b, alkbh9a and alkbh9b mutant plants using two sgRNAs on each gene (Figure 5c; Figures S12 and S13). As a result, alkbh10b mutants with frame‐shifted protein sequence (Figure S14) showed higher m6A/A ratios than WT under normal water condition, and even increased to relatively high levels under drought stress condition in LC–MS/MS assays (Figure 5d). We also detected the levels of m6A modification using m6A dot blots, which confirmed the m6A demethylase activity that the m6A levels was significantly increased in the alkbh10b mutants compared with WT (Figure 5e). Accordingly, alkbh10b mutants showed drought resistance phenotypes compared with WT for the withering cotyledons and leaves on WT plants (Figure 5f). Both malondialdehyde (MDA) and hydrogen peroxide (H2O2), two important physiological indices of cell damage under drought stress, were significantly lower in alkbh10b than WT after drought which indicated alkbh10b plants suffered lower degree of stress injury (Figure 5g, h). While, alkbh10a mutants also displayed higher m6A/A ratios than WT under normal water condition, and also increased to high levels under drought stress condition, but without obvious phenotypic difference on drought resistance between mutant plants and WT plants (Figure S13). However, the m6A/A ratios in alkbh9a and alkbh9b did not exhibit significant differences compared with WT, and no significant difference phenotype on drought resistance and physiological indicators under drought stress (Figure S15). These results indicate that GhALKBH10B is an essential demethylase in cotton which might contribute to drought resistance. However, this does not prove that GhALKBH9A and GhALKBH9B has no RNA demethylation activity, since functional redundancy is common in the tetraploid cotton genome.

We examined the regulatory function of ALKBH10B in drought response pathways with focusing on the enriched stress‐responsive signal transduction pathways. Three transcripts, annotated as protein phosphatase 2C 37‐like protein (PP2CA), 9‐cis‐epoxycarotenoid dioxygenase 4 (NCED4) and zeaxanthin epoxidase (ZEP) (ABA1) isoform 1, were enriched in the ABA biosynthesis and signalling pathway. Both GhZEP and GhPP2CA had an enriched m6A peak located in the 1st exon in IP samples under drought stress, and GhNCED4 had an enriched peak located in both first exon and 3'UTR (Figure 6a–c). The expression levels of these three genes were upregulated under drought stress in both varieties or at least in one variety (Figure 6d), indicating that at least some ABA‐related genes are induced by drought stress. At the same time, the m6A enrichment of these three genes was also upregulated under drought condition compared to watered control, but with different fold change for different genes in either ZY168 or ZY7 (Figure 6e). For instance, GhZEP was significantly upregulated in mRNA expression and m6A enrichments in ZY168, while no significance been observed in ZY7. The expression of GhNCED4 was significantly upregulated in both varieties with little difference between ZY7 and ZY168. The expression of GhPP2CA was also significantly upregulated in both varieties but only showing significance of m6A enrichment in ZY7 under control and drought condition. (Figure 6d, e). We further confirmed that these three ABA‐related genes exhibited relatively higher levels of m6A enrichment in alkbh10b plants than control, both in alkbh10b‐#1 and alkbh10b‐#8 (Figure 6f), but only GhPP2CA presented higher levels in alkbh10a than WT (Figure S17). To better understand how m6A methylation affects transcript abundance, we investigated the mRNA stability of GhZEP, GhPP2CA and GhNCED4 by calculating the mRNA degradation rate using actinomycin D to inhibit transcription. It was seen that mRNA degraded after 3 h treatment, and the degradation rate decreased for mRNAs of GhZEP and GhPP2CA in alkbh10b mutant plants (Figure 6g, h), but only a little difference was seen for GhNCED4 in alkbh10b‐#8 (Figure 6i). Concomitant with m6A deposition in the mutant plants, m6A slowed down degradation and therefore stabilized mRNAs for these key ABA‐related genes (Figures 6j and S16). We also identified some DMeETs related to Ca2+ signal transduction. GhANN1 (cotton annexin D1), GhCML13 (Probable calcium‐binding protein CML13), GhECA1 (Calcium‐transporting ATPase 1, endoplasmic reticulum‐type‐like protein), and GhCNGC4 (Cyclic nucleotide‐gated ion channel 4‐like protein) were confirmed to be enriched in m6A modification in both Ghalkbh10b‐#1 and Ghalkbh10b‐#8 plants (Figure 7a), but only GhECA1 and GhCNGC4 showed higher m6A levels in alkbh10a than WT (Figure S17). The methylated regions were in the 3'UTR of these genes except for GhCML13, where it was located in the 1st exon near to the 5'UTR (Figure 7b–e). The stability of these genes' mRNA was also significantly increased in knockout plants under actinomycin D treatment, similar to the findings for ABA‐related genes (Figures 7f–i and S16). These results suggest that GhALKBH10B contributes to mRNA stability of genes in the Ca2+ and ABA signalling pathways.

Further functional validation of the m6A targeted Ca2+ signal related genes was carried out by virus‐induced gene silencing (VIGS) technology. GhECA1 and GhCNGC4 were two key genes involved in Ca2+ signal pathway in response to drought that GhECA1 was downregulated and GhCNGC4 was upregulated under drought stress in two varieties (Figure S18). By applying drought stress to the seedlings, plants silenced for TRV::GhECA1 and TRV::GhCNGC4 exhibited earlier wilting compared to controls (TRV::00; Figure 8a, b) with the decreased m6A enrichment on silenced genes (Figure 8c). PD6, a well‐illustrated image‐based trait, which could indicate the drought ability of cotton plants with high value presenting lower drought resistance ability in our previous paper (Li et al., 2020), was significantly higher than that of control plants (Figure 8d). Meanwhile, both MDA and H2O2 showed significant increases in the silenced TRV::GhECA1, TRV::GhCNGC4 and TRV::00 plants under drought compared to controls (Figure 8e, f). Moreover, the roles of these two genes in mediating calcium ion flow were further investigated by Non‐invasive Micro‐test Technology (NMT). Ca2+ flux direction was changed in the mesophyll cells of silenced plants, and specific Ca2+ exhibited efflux in silenced plant while controls retained influx after drought stress, despite a higher influx than in normal water conditions (Figure 8g), which also indicated the decreasing drought resistant in silenced plants. This demonstrates the positive regulatory functions of GhECA1 and GhCNGC4 in the drought response.

Drought response is highly regulated at the genomic and transcriptional levels (Gupta et al., 2020; Wang et al., 2016a). Due to the complexity of cotton genome (Chen et al., 2020), the drought response mechanism in cotton is poorly defined compared with food crops such as rice, maize and wheat. Previously identified changes in transcription and epigenetic modification during drought stress have greatly improved our understanding of the mechanism related to cotton drought response (Li et al., 2020; Liang et al., 2016; Liu et al., 2014a; Shang et al., 2020). In the present study, we show that a considerable number of mRNAs exhibit drought‐regulated methylation, with differential levels of m6A modification under drought stress in both drought‐resistant and drought‐sensitive varieties (Figure 1). These observations suggest that m6A modification might be another important contributor to the drought response mechanism. The overall m6A levels was slightly increased under mild drought condition (MD), but significantly increased under severe drought conditions (SD), up to the point of recovery following re‐watering (RD), both in drought‐resistant and drought‐sensitive varieties; these are different to what reported in maize and sea buckthorn (Miao et al., 2020; Zhang et al., 2021). Moreover, drought triggered different m6A modifications in drought‐resistant versus drought‐sensitive varieties under different drought conditions. m6A peaks were enriched mostly in the 3'UTR in cotton (with the proportion > 60%), consistent with previous reports in other plants (Zhang et al., 2019; Zheng et al., 2021; Zhou et al., 2022). The drought‐sensitive variety ZY7 exhibited a similar proportion of m6A modification under MD and SD. While, the distribution of m6A modification changed dramatically, with the enrichment of m6A near the 5'UTR region increasing significantly in drought‐resistant variety ZY168 under MD and SD (Figure 2c). Accordingly, we identified more transcripts within 5'UTR deposited m6A peaks in drought‐resistant variety ZY168 (762 peaks on 594 transcripts) than in drought‐sensitive variety ZY7 (310 peaks on 288 transcripts) (Tables S4 and S5). GO and KEGG analyses revealed that numerous genes regulating ion transport, photosynthesis, or carbon metabolism acquire m6A marks in a drought‐dependent manner, with the expression of 748 genes under tight regulation by their methylation status (Table S14). We found that some genes related to stress response and signalling, including ABA biosynthesis and signalling and Ca2+ transport and signalling genes, undergo m6A‐mediated post‐transcriptional regulation during drought response in cotton. In the seven genes we selected, four genes (GhZEP, GhPP2CA, GhNCED4 and GhCML13) harboured m6A site near the 5'UTR (in the 5'UTR or in the 1st exon), with three genes (GhZEP, GhNCED4 and GhCML13) harbouring cotton specific CAAUG motif. The m6A levels of these genes increased in alkbh10b plants compared to the control. Increased m6A abundance is known to improve mRNA stability for some key stress‐responsive genes in Arabidopsis (Anderson et al., 2018), sorghum (Zheng et al., 2021) and apple (Hou et al., 2022). Meanwhile, there also appears to be a complex relationship between the positions of m6A peaks in RNAs and their abundance, with positive relations in mRNA levels with m6A modification regulated by water conditions either in drought‐sensitive or drought‐resistant varieties. The m6A modification for GhNCED4 was both enriched in the 5'UTR and 3'UTR region, while the peak was in the 1st exon close to the 5'UTR for GhZEP and GhPP2CA. These genes showed upregulated 5'UTR m6A modification, which was accompanied by up‐regulation of mRNA levels under drought stress (Figure 6). While, downregulated 5'UTR m6A modification of GhCML13 resulting in the down‐regulation of mRNA levels (Figure 7). The m6A modification for GhCNGC4, GhECA1 and GhANN1 was mainly enriched in the 3'UTR region, with upregulated m6A modification of GhCNGC4 and GhANN1 linked to the up‐regulation of m6A levels after drought stress (Figure 7). These results indicate that the different distributions of m6A modification under different water conditions in different varieties, which might contribute to the different ability of drought resistance in different cotton varieties, and this was not reported in any other study. The dynamic RNA methylation levels during drought stress in cotton reveal its importance to maintain the balance of RNA methylation under drought stress, which is accomplished by both m6A methyltransferases and demethylases. Understanding the functional roles of reversible m6A methylation in mammals has been paramount since the discovery of the first human mRNA m6A demethylase FAT MASS AND OBESITY‐ASSOCIATED (FTO) (Jia et al., 2011). In Arabidopsis, there are 13 potential RNA demethylases, among which ALKBH10B functions in floral transition by regulating the methylation of key flowering‐related mRNAs for FT, SPL3 and SPL9 (Duan et al., 2017), while ALKBH9B modulates infection by alfalfa mosaic virus (Martinez‐Perez et al., 2017) and ALKBH6 functions in growth and stress responses (Huong et al., 2020). There are fewer RNA demethylases elucidated in other plants. Tomato DNA methylation‐dependent RNA demethylase SlALKBH2 functions in fruit ripening by regulating the abundance of SlDML2 mRNA (Zhou et al., 2019). ZmALKBH10 was induced by drought stress and acts to decrease m6A abundance in maize (Miao et al., 2020). A search of the upland cotton genome led to the identification of 22 homologous ALKBH genes (Table S13). Gene expression analysis indicated that the expression of GhALKBH10B increases dramatically after drought stress in cotton, which was consistent with the reported in maize and sea buckthorn (Miao et al., 2020; Zhang et al., 2021). However, there was unharmony that the overall methylation level also increased after drought stress in cotton, implying the complexity in studying m6A in response to drought, potentially due to be co‐regulated by methyltransferases. Consequently, RNA m6A methyltransferase is another important factor for the maintenance of RNA methylation at appropriate levels (Chen et al., 2018; Lu et al., 2020; Pendleton et al., 2017; Zhang et al., 2019; Zhou et al., 2021). Recently, MdMTA was reported to enhance drought tolerance by increasing m6A modification on drought‐responsive genes in apple by using overexpression and RNAi technologies (Hou et al., 2022). Here, 6 homologous METTL3 (AtMTA) genes and 2 homologous WTAP (AtFIP37) genes were identified in upland cotton genome (Table S15). We created mutant plants of the homologous methyltransferase genes (GhMTA and GhFIP37) in cotton, but the T0 generation of knockout transgenics (mta and fip37) grew weakly and exhibited abnormal reproductive growth, with poor seed set, and even unbolling (Figure S19), which was consisted with their functions in reproductive development in rice (Zhang et al., 2019). Dose‐dependent mutation by manipulating on promoters was needed for further study. Rapid advances have been made in understanding the functional diversity and regulatory mechanisms of m6A marks on abiotic stress response in plants. In the recent study, we generated a transcriptome‐wide m6A map in cotton and uncovered the dramatic patterns of m6A modifications associated with drought response using MeRIP‐seq data in cotton. We propose the regulatory role of m6A methylation in drought response of cotton by targeting some vital stress‐related genes such as genes in Ca2+ and ABA signalling pathways. Our study corroborated and complemented the conclusion that RNA methylation facilitates mRNA stability contribute to drought resistance as reported in apple and other plants, and even better to confirmed function of m6A on the stability of Ca2+ and ABA signalling pathway related genes mediated by GhALKBH10B. Other than the rough map in other studies, we used different water condition and varieties with different drought resistant ability to generate the global map. We highlight that increasing of m6A enrichment near the 5'UTR region under drought stress in drought‐resistant variety might contribute to the different ability of drought resistance of different cotton varieties. Moreover, this study not only expended our understanding of ALKBH10B functioned with affecting floral transition in Arabidopsis, but also provides new insights into the drought response network as well as new strategy for drought resistance breeding of cotton. Both writers and erasers had the global impact on whole transcriptome, and knock out them using CRISPR/Cas9 could lead to changes in m6A modification on many transcripts. Future research should focus on site‐specific m6A modifications of targeted mRNA, which could precisely decipher the regulated mechanism of m6A in diversiform biological process.

Drought experiments were carried out using two representative cotton varieties of drought‐resistant ZY168 and drought‐sensitive ZY7 as our previous study (Li et al., 2020). Germinated cotton seedlings were planted in pots filled with mixed nutrient soil in the greenhouse. Drought stress was imposed by withholding water supply continuously after the first leaf had fully spread, at three time‐points, that is, mild drought (MD), severe drought (SD), and re‐watering after drought (RD). Absolute soil water content (ASWC) was used to evaluate the different drought degrees at these time‐points. At the beginning of the experiment, the weight of dry mixed nutrient soil (4.5 kg) and pot (0.2 kg) was 4.7 kg for each. To record the ASWC during the growth of cotton, the whole weight (WW1) including pot and soil was weighed every day. The ASWC (%) was calculated as (WW1‐4.7) / 4.5 × 100. And the ASWC at ca. 10% for MD, at ca. 5% for SD, ASWC and at ca. 15% for RD. At the same time, a parallel control panel with normal water supply was maintained at an ASWC of ca. 15% at each point. Newly emerged second true leaves were sampled and mixed at MD, SD and 6 days after RD. Leaf samples were obtained at the same time (11:00 am), immediately frozen in liquid nitrogen, and stored at −80 °C until sequencing.

Total RNA was extracted from the leaves of drought‐treated and parallel control cotton plants under MD, SD and RD conditions using the TG‐DP441 RNA Kit (TIANGEN, Beijing, China) following the manufacturer's procedure. The quality and quantity of total RNA were analysed using a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number >7.0. The total m6A/A ratio in cotton was determined by LC–MS/MS as previously described (Zhou et al., 2019) with minor modifications. mRNA was isolated from total RNA using Dynabeads Oligo(dT)25 (Invitrogen, 61 006). 200 ng was digested to 5'‐mononucleotides with 1 unit of Nuclease P1 (NEB, M0660S) in 50 μL reaction buffer (10 mM ammonium acetate, pH 5.3, 25 mM NaCl, and 2.5 mM ZnCl2) for incubation at 37 °C for 6 h. Then the reaction product was mixed with 1 unit of alkaline phosphatase (Sigma‐Aldrich, P6774‐1) and 5 μL 1 M fresh NH4HCO3, and incubated at 37 °C for another 6 h. The supernatant of digested sample was collected and diluted to a suitable concentration for LC–MS/MS assay following a published method (Zhou et al., 2018) with minor modifications. In brief, nucleoside standards including N6‐methyladenosine and adenosine nucleoside were used to determine the peak position of each nucleoside. Nucleosides were separated by reverse‐phase ultra‐performance liquid chromatography on a C18‐ODS column (Shimadzu). The AB SCIEX Triple Quad 5500 System liquid chromatography mass spectrometer was used for online mass spectrometry detection, set to multiple reaction monitoring in positive electrospray ionization mode. The mobile phase consisted of buffer A (5 mM ammonium acetate) and buffer B (100% acetonitrile). Nucleosides were quantified using the nucleoside‐to‐base ion mass transitions of 282.0–150.1 (m6A) and 268.0–136.0 (A). Standard curves were generated by running a concentration series of pure commercial A (TargetMol, T0853) and m6A (TargetMol, T6599). Quantification was determined by comparison with the standard curves obtained from nucleoside standards running at the same time. Three independent biological replicates were performed for each sample.

The MeRIP‐seq was based on a previously described method with minor modifications (Dominissini et al., 2013). Briefly, mRNA was purified by poly‐T oligo, then fragmented and constrated input and IP libraries for sequencing through Illumina Novaseq™ 6000 platform. The detailed method for MeRIP‐seq and related data analysis was supplied in Method S1.

Blastp (identity >40%, e‐value <1e‐30, score > 150) was used to identify homologous m6A methylases in upland cotton TM‐1 (Wang et al., 2019). A phylogenetic tree was construct using protein sequences with the software MEGA‐X (Kumar et al., 2013) with the Tets‐Neighbour‐Joing method, and the calculated result was exported to iTOL (http://itol.embl.de/) for drawing. CRISPR‐Cas9 gene editing technology was used to construct vector for GhALKBH9A (Ghir_A04G014040), GhALKBH9B (Ghir_D11G018740), GhALKBH10A (Ghir_D11G012740) and GhALKBH10B (Ghir_D08G007610), using published procedures (Wang et al., 2018). In brief, targets sequences were searched from the CRISPR‐P website (http://cbi.hzau.edu.cn/crispr/), and each pair of sgRNAs, named sgRNA1 and sgRNA2, were designed to be integrated in a demethylase vector. Fragments containing tRNA‐sgRNA1 fusion and gRNA‐tRNA‐sgRNA2 fusion amplified from pGTR plasmid were fused together using an overlapping extension PCR. The tRNA‐sgRNA1‐gRNA‐tRNA‐sgRNA2‐gRNA assembled fragments were finally ligated to the linearized pRGEB32‐NPT II‐U6.7 vector for creating demethylase knockout vectors. The primers used in vector construction are listed in Table S16. All the constructed plasmids were separately introduced into Agrobacterium tumefaciens strain GV3101 by electroporation. The hypocotyls of wild‐type cotton (Jin668) cultured in a 28 °C incubator for approximately 7 days were cut into stem segments of 0.5–0.8 cm, and were used transferred into the activated agrobacterium tumefaciens. After 3–5 min of infection, the hypocotyls were dried slightly, placed on the co‐culture medium with filter paper for 36 h at 20 °C. And then cultured on the selective medium containing appropriate antibiotics until embryogenic callus grows, which were then selected into the differentiation medium to make it grow into embryoids, and then develop into transgenic cotton seedlings. The cotton genetic transformation method is as described previously (Jin et al., 2006). Mutants were detected using the published Hi‐TOM program (Liu et al., 2019).

Total RNA was extracted from leaves in WT (Jin668) and mutant plants (alkbh9a, alkbh9b, alkbh10a and alkbh10b) at 28 days after seed germination. mRNA was isolated from total RNA using Dynabeads Oligo(dT)25 (Invitrogen, 61 006) according to the kit instructions. 2 μg poly(A)+ mRNA was randomly fragmented into ~300 nt long fragments by an incubation at 94 °C for 30 s in fragmentation butter, and the reaction terminated by addition of EDTA, followed by phenol‐chloroform extraction and ethanol precipitation to purify the fragmented mRNAs according to Zhou et al. (2019). Fragmented mRNA was precipitated and resuspended in 150 μL RNase‐free water. NanoDrop spectrophotometry was used to measure concentrations of fragmented mRNA, and fragmented mRNA size was validated in a 1.5% (wt/vol) agarose gel (Figure S20). 10 μL fragmented mRNA was retained and used as the input sample (control). 100 μL fragmented mRNA was incubated with 5 μg m6A‐specific antibody at 4 °C for 2 h, which was later used to perform MeRIP experiments according to the instructions for the Magna MeRIP m6A Kit (Merck, A‐17‐10 499). Fold enrichment of the specific region of transcripts was determined by qRT‐PCR using 2−ΔCT methodology, the value of ΔCT (cycle threshold) was normalized by ΔCTIP, ΔCTinput and input dilution factor. The designed primers for m6A‐IP‐qPCR were shown in Table S16. Statistical significance was determined by a two‐sided t‐test on three biological replicates: *P < 0.05, **P < 0.01.

Dot blot assay was performed as described previously with minor modifications (Ren et al., 2022). Total RNA from the WT and the alkbh10b mutant were extracted and heated at 95 °C for 5 min to remove the secondary structure and then was chilled on ice. RNA was spotted onto Hybond™‐N+ membrane. The membrane was crosslinked using an HL‐2000 HybriLinker for 5 min and then washed with TBST (TBS with Tween‐20) for 5 min to remove the unbound RNA and blocked in 5% non‐fat milk for 1 h at room temperature. The membrane was then incubated with an m6A antibody (1:2000; ABclonal) overnight at 4 °C. The membrane was washed using TBST and incubated with the secondary antibody for 1 h at room temperature. Finally, the super ECL detection reagent (Yeasen) was used for visualization.

The procedure is based on the previously described method (Duan et al., 2017) with minor modifications. Briefly, 28‐d‐old cotton leaves of wild type (Jin668) and mutant plants (alkbh10b‐#1 and alkbh10b‐#8) were cut into leaf discs and then added to actinomycin D (AbMole, M4881) solution with a final concentration of 20 μg mL−1. After infiltration for 1 h, five leaf disks were collected and considered as time 0 controls, and subsequent samples were harvested every 3 h in triplicate, and qRT‐PCR conducted with ABI 7500 system (Applied Biosystems, Foster city, CA) was used to determine mRNA expression levels, with the internal control GhUBQ7. The 2−ΔCT method was used to present relative changes in gene expression levels (Schmittgen and Livak, 2008). The mRNA degradation was calculated by the percentage of gene expression at sampled time relative to control. All primers used in transcription inhibition assay are listed in Table S16.

Drought stress was imposed to the seedlings of generation T2 of mutant lines (alkbh10a‐#2, alkbh10a‐#3, alkbh10b‐#1, alkbh10b‐#8, alkbh9a‐#4 and alkbh9b‐#1) by withholding water supply continuously after 10 days post germination in an environments controlled chamber (CONVIRON, MTR30‐190316) at temperature of 28–30 °C with 16 h light and 8 h dark (30% humidity, 100 μMol multi‐chip LED lights) for growing and natural drought treatment. The phenotypes and physiological indices were taken at time of the ASWC (%) at about 12% in treatment panels. At the beginning of experiment, all cotton seeds were sown in the same pot (55 g nutrient soil and 15.4 g pot) with saturated water in a same big container. The whole weight (WW2) including pot and soil was weighed every day. The ASWC (%) was calculated as (WW2‐70.4) / 55 × 100. Physiological indices of MDA and H2O2 were measured by kits provided by COMIN Technology CO., Ltd. (Suzhou, China), 4–5 seedlings were pooled for each biological replicate and at least three biological replicates were measured in experiments. Statistical significance was determined by a two‐sided t‐test on three biological replicates: *P < 0.05, **P < 0.01. Virus‐induced gene silencing (VIGS) technology was used to verify the function of m6A modified transcripts under drought treatment in cotton. The procedure of VIGS was based on a previous method (Gao et al., 2013). 300–500 bp fragment of the target gene (GhECA1 or GhCNGC4) was cloned into the pTRV2 vector, and pTRV2 vector (empty vector) without target gene was used as a control. Both vectors were transformed into A.tumefaciens GV3101 by electroporation. Agrobacterium solution containing pTRV1 or pTRV2 vectors was cultured overnight in liquid LB medium with kanamycin 50 mg/L and rifampicin 25 mg/L, and collected by buffer solution (10 mmol/L MES, 10 mmol/L MgCl2, 40 μg/mL AS, pH 5.6) and resuspended and adjusted to OD600 of 0.6–0.8, then pTRV1 and pTRV2 (TRV:GhECA1 or TRV:CNGC4) were mixed in equal quantities and placed in a 28 °C shaker at 150 r/min for activation for 2 h, before agroinfiltration into the cotyledons of Jin668 by syringe inoculation when the two cotyledons had fully spread (12 days after sowing). TRV:CLA1 (GhCLA1, a key enzyme for chloroplast synthesis) plants were used as a positive control. After covering with black plastic film (removed after 24 h), these inoculated plants were placed to the controlled chamber, at temperature of 25 °C, light and humidity was the same to drought treatment described above. Then, watering to saturation for the initiation of natural drought treatment. The phenotypes and physiological indices of MDA and H2O2 were analysed at the time of the ASWC at about 12.7%. The image‐based traits were extracted for the images followed our previous methods (Li et al., 2020).

Net Ca2+ flux of mesophyll cells was measured using NMT (NMT‐YG‐100; Younger USA, Amherst, MA, USA) as previously described with minor modifications (Li et al., 2021). Leaves sampled from cotton following VIGS were immobilized in the measuring buffer (0.1 mM KCl, 0.1 mM CaCl2, 0.1 mM MgCl2, 0.5 mM NaCl, 0.3 mM MES, 0.2 mM Na2SO4, pH 6.5) for 10 min to recover from mechanical damage, and the steady‐state ion fluxes in leaf mesophyll cells were then continuously recorded for the indicated times at room temperature. During measurements, the ion microsensor was placed about 3 μm above the surface of a cell and moved repeatedly between two points (10 μm in distance) perpendicular to the cell (Figure S21). The results of net Ca2+ ion fluxes were calculated by imFluxes V2.0 software (YoungerUSA LLC, Amherst, MA).

XY, XZ and XN designed and guided the project. KL and MW gave comments on the project. BL and MZ performed experiments and wrote the manuscript. BL and BZ carried out the drought treatment experiments. DY help creating the mutants plants. WS help performing the NMT experiments. LD help extracting imaging traits. BL and YM analysed the data. XY, KL and XZ revised the manuscript.

This work was supported by grants from the National Natural Science Foundation of China (32201725 and 32171942), the Bintuan Science and Technology Program (2022DB012), the Natural Science Foundation of Hubei Province (2022CFB756), the China Postdoctoral Science Foundation (2022M711279) and Fundamental Research Funds for the Central Universities (2662020ZKPY011).

The authors declare that they have no competing interests.

Click here for additional data file. Click here for additional data file. Click here for additional data file.