Next Article in Journal
Book Review: Teucrium Species: Biology and Applications; Stanković, M., Ed.; Springer Nature: Cham, Switzerland, 2020; ISBN: 978-3-030-52158-5
Previous Article in Journal
Diverse Functions of Plant Zinc-Induced Facilitator-like Transporter for Their Emerging Roles in Crop Trait Enhancement
Previous Article in Special Issue
The Intragenesis and Synthetic Biology Approach towards Accelerating Genetic Gains on Strawberry: Development of New Tools to Improve Fruit Quality and Resistance to Pathogens
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Plants
Volume 11
Issue 1
10.3390/plants11010103
plants-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Profiling of N6-Methyladenosine (m6A) Modification Landscape in Response to Drought Stress in Apple (Malus prunifolia (Willd.) Borkh)

by
Xiushan Mao
1,
Nan Hou
2,3,
Zhenzhong Liu
2,3,* and
Jieqiang He
2,3,*
1
Shandong Transport Vocational College, 7369 Bohai Road, Weifang 261206, China
2
State Key Laboratory of Crop Stress Biology for Arid Areas, Yangling, Xianyang 712100, China
3
Shaanxi Key Laboratory of Apple, College of Horticulture, Northwest A&F University, Yangling, Xianyang 712100, China
*
Authors to whom correspondence should be addressed.
Plants 2022, 11(1), 103; https://0-doi-org.brum.beds.ac.uk/10.3390/plants11010103
Submission received: 18 November 2021 / Revised: 15 December 2021 / Accepted: 20 December 2021 / Published: 30 December 2021
(This article belongs to the Special Issue Research Progress on the Application of Biotechnology in Fruit Crop Breeding)
Browse Figures
Versions Notes

Abstract

:
Drought stress is a significant environmental factor limiting crop growth worldwide. Malus prunifolia is an important apple species endemic to China and is used for apple cultivars and rootstocks with great drought tolerance. N6-methyladenosine (m6A) is a common epigenetic modification on messenger RNAs (mRNAs) in eukaryotes which is critical for various biological processes. However, there are no reports on m6A methylation in apple response to drought stress. Here, we assessed the m6A landscape of M. prunifolia seedlings in response to drought and analyzed the association between m6A modification and transcript expression. In total, we found 19,783 and 19,609 significant m6A peaks in the control and drought treatment groups, respectively, and discovered a UGUAH (H: A/U/C) motif. In M. prunifolia, under both control and drought conditions, peaks were highly enriched in the 3′ untranslated region (UTR) and coding sequence (CDS). Among 4204 significant differential m6A peaks in drought-treated M. prunifolia compared to control-treated M. prunifolia, 4158 genes with m6A modification were identified. Interestingly, a large number of hypermethylated peaks (4069) were stimulated by drought treatment compared to hypomethylation. Among the hypermethylated peak-related genes, 972 and 1238 differentially expressed genes (DEGs) were up- and down-regulated in response to drought, respectively. Gene ontology (GO) analyses of differential m6A-modified genes revealed that GO slims related to RNA processing, epigenetic regulation, and stress tolerance were significantly enriched. The m6A modification landscape depicted in this study sheds light on the epigenetic regulation of M. prunifolia in response to drought stress and indicates new directions for the breeding of drought-tolerant apple trees.

1. Introduction

According to the central dogma, RNAs are the essential and fundamental components responsible for the transfer of genetic information from DNA to proteins. During this process, over 100 distinctly chemical modifications have been reported to modify various kinds of RNAs in all living species [1]. N6-methyladenosine (m6A) RNA methylation is a crucial internal modification and is found to occur in rRNA, mRNA, tRNA, miRNA, and long non-coding RNA [2,3,4,5]. Moreover, m6A modification accounts for 80% of all RNA methylation modifications [6]. In 1974, m6A was discovered as a dominant type of mRNA methylation in mammals for the first time [7]. Nowadays, m6A modifications have been widely reported in various species, such as viruses, plants, yeast, humans, and other mammals [3,5]. Among plants, m6A was identified in wheat (Triticum turgidum L.), oat (Avena sativa L.), and maize (Zeamays L.) about 40 years ago [8,9,10]. m6A is a dynamic and reversible modification process that requires three effectors: “writer”, “reader”, and “eraser” proteins. Writers carry out m6A modification, readers recognize the methylation, and erasers demethylate m6A modifications [11,12]. In the current study, writer protein complexes include METTL3, METTL14, WTAP, etc.; readers are mainly the YTH-domain-containing proteins; and erasers include FTO and ALKBH5 [13,14,15,16]. Additionally, studies have shown that m6A is involved in nuclear–cytoplasmic export, RNA stability, pre-mRNA splicing, primary microRNA processing, alternative polyadenylation site choice, and translation efficiency in mRNA metabolic processes [17,18,19,20,21,22].
Recently, with the development of m6A sequencing (m6A-seq) technology, an increasing number of comparative m6A methylome studies have been conducted to better understand its role in plant biological processes. In Arabidopsis thaliana, m6A regulates leaf morphology [23], trichome development [24], floral transition [25], and embryonic development [22]. In addition, m6A regulates microspore degeneration in rice [26] and is responsible for the fruit ripening of tomato [27] and strawberry [28]. Stress responses are also affected by m6A modification. In Arabidopsis, YTH-domain proteins evolutionarily conserved the C-terminal region 1 (ECT1) and ECT2 which interact with calcineurin B-like-interacting protein kinase1 (CIPK1) and mediate calcium signaling under various stresses [29]. Compared with wild-type Arabidopsis, the alkbh6 mutant plants exhibit low survival rates under abiotic stresses, including salt, drought, and heat stresses [30]. Increased levels of m6A methylation have been observed in rice plants’ response to viral infection [31]. In maize, m6A hypomethylation under drought stress has a favorable function in drought response [32]. The function of m6A modification in pak choi (Brassica rapa ssp. chinensis) under heat stress has also been investigated [33]. A recent study on apple found that the m6A reader YTH domain-containing RNA binding protein 2 (YTP2) regulates Mildew Locus O 19 (MdMLO19) mRNA stability and translation efficiency of antioxidant genes to confer powdery mildew resistance [34]. Although m6A has been reported in both biotic and abiotic stresses, its role in drought stress in non-model plants is currently unknown.
Extreme climate change causes frequent global droughts and high temperatures, which severely limit crop growth and yield [35]. Apple is one of the world’s popular fruits; however, its production and quality are frequently threatened by drought stress [36,37]. Apple propagation mainly relies on vegetative propagation via grafting and budding. M. prunifolia is a wild relative of apple with strong biotic and abiotic resistance to drought, cold, heat, and disease [38], which makes it one of the best rootstocks in northwest China, where apple cultivars are usually grafted to vigorous rootstocks, including M. prunifolia and Malus sieversii. In addition, M. prunifolia is commonly used as a parent in cross breeding studies for stress tolerance [39]. Despite the outstanding performance of M. prunifolia in improving apple drought tolerance, the molecular mechanism of M. prunifolia in response to drought is largely unclear.
In this study, we first performed a transcriptome-wide m6A modification profile of M. prunifolia seedlings and investigated changes in m6A modification after drought stress. We also performed an RNA-seq analysis and identified differentially expressed genes (DEGs) in M. prunifolia in response to drought. To investigate the potential relationship between m6A levels and gene expression levels in M. prunifolia in response to drought stress, we performed association analysis between differential m6A peaks and DEGs. Our data allowed the identification of some drought-responsive genes along with changes in gene expression by m6A modifications in M. prunifolia after drought stress, such as Heat shock protein 60 (HSP60), jasmonate-Zim-domain protein 3 (JAZ3), Scarecrow-Like 1 (SCL1), and ETHYLENE RESPONSE FACTOR1 (ERF1). Overall, our work illustrates the m6A modification landscape of M. prunifolia in response to drought stress and provides new insights into the molecular mechanisms operating in conditions of drought.

2. Materials and Methods

2.1. Plant Materials and Stress Treatment

One-year-old M. prunifolia seedlings were used as the plant materials. Seeds of M. prunifolia ‘Fupingqiuzi’ were collected from Fuping (Weinan, Shannxi, China) and stratified in wet sand at 4 °C for three months. Then, the seeds were sowed in a plant growth chamber with 8000 lux light intensity, 14 h light/10 h dark photoperiods at 25 °C. Three months later, the M. prunifolia seedlings were moved to a greenhouse at the Northwest Agriculture and Forest University, Yangling (34°20′ N, 108°24′ E), Shaanxi Province, China. The seedlings were transplanted into plastic pots (15 cm × 20 cm, ~1.3 L) filled with a mixture of garden soil and substrate (PINDSTRUP, Denmark) (1:1, v/v). The stress treatment was started a year later, when the seedlings were 1.5 m tall. Forty-two seedlings with uniform growth were chosen and every seventh seedling were used as a biological replicate. When the treatment began, all the seedlings were watered until saturated (control) and then water was withheld from half the plants until the relative soil water content reached approximately 40% (drought treatment). Mature leaves were collected from the middle of the trees for the following RNA extraction.

2.2. RNA-Seq Analysis

Leaves were collected from control and drought-treated M. prunifolia seedlings. Total RNA was extracted using the cetyltrimethylammonium bromide (CTAB) method according to a previously described procedure [40]. The RNA-seq library was constructed as previously reported by Xie et al. [41]. RNAs were subjected to sequencing on the Illumina HiSeq 4000 platform by Novogene (Beijing, China). Sequences were aligned to a recently released Malus × domestica genome sequence (GDDH13 version 1.1, https://iris.angers.inra.fr/gddh13/downloads/GDDH13_1-1_formatted.fasta.bz2, accessed on 18 November 2021.) [42] using HISAT2 v2.1.0. BAM conversion, sorting, and indexing were performed using SAMtools v1.9. Read counting within genes wase analyzed with HTSeq v0.12.4 using the gene annotation file (https://iris.angers.inra.fr/gddh13/downloads/gene_models_20170612.gff3.bz2, accessed on 18 November 2021.) [43]. Differences in gene expression were analyzed by DESeq2 v1.30.1 with a threshold of an adjusted p-value below 0.05 and |log2(fold change)| greater than 1 [44]. Length of genes were calculated by GenomicFeatures v1.42.3, and fragments per kilobase of transcript per million fragments mapped (FPKM) values were obtained by TBtools [45,46]. Heatmaps of gene expression levels were plotted using pheatmap v1.0.12 [47]. Gene Ontology (GO) enrichment analyses were performed using agriGO v2.0 and clusterProfiler v.3.18.1 [48,49].

2.3. m6A-Seq Analysis

mRNA m6A was sequenced by MeRIP-seq at Novogene (Beijing, China). Briefly, a total of 300 μg RNA was extracted from the leaves. The integrity and concentration of extracted RNAs were detected using an Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA, USA) and simpliNano spectrophotometer (GE Healthcare, Chicago, IL, USA), respectively. Fragmented mRNA (~100 nt) was incubated for 2 h at 4 °C with anti-m6A polyclonal antibody (Synaptic Systems, Göttingen, Germany) in the immunoprecipitation experiment. Then, immunoprecipitated mRNAs or Input was used for library construction with NEBNext ultra-RNA library prepare kit for Illumina (New England Biolabs, Ipswich, MA, USA). The library preparations were sequenced on an Illumina Novaseq platform with a paired-end read length of 150 bp according to the standard protocols. The sequencing was carried out with three independent biological replicates.
Raw reads from m6A-seq were trimmed to remove adaptor sequences and bases with a quality lower than 20 using Trimmomatic v0.39 and FastQC v0.11.9 [50,51]. The remaining reads were mapped onto the apple reference genome by HISAT2 v2.1.0 [52]. Post-processing was carried out by SAMtools v1.9 [53]. Peak calling was analyzed by exomePeak2 v1.2.0 with a p-value below 0.05 and a log2(fold change) greater than 1 [54]. The other running parameters of exomePeak2 were set as: fragment_length = 100, binding_length = 25, step_length = 25, peak_width = 50. The overlapping peaks of each biological replicate and a Venn diagram were generated by intervene v0.6.5 [55]. FindMotifsGenome.pl in HOMER v4.10.0 was employed to identify the m6A motifs [56]. Differentially methylated peaks were identified using exomePeak2 with a threshold of an adjusted p-value below 0.05 and |DiffModLog2FC| above 0.5. The CMRAnnotation tool in PEA v1.1 and bedtools were used to annotate the peaks’ different transcript distributions using the gene annotation file [57,58]. The visualization of m6A peaks was performed using the Integrative Genomics Viewer v2.10.2 [59].

3. Results

3.1. Transcriptome-Wide Mapping of m6A in Malus prunifolia Seedlings

To investigate whether m6A methylation participates in drought stress in apples, we constructed and sequenced a series of m6A-immunoprecipitation (IP) and matched input libraries to obtain the drought and control-treated M. prunifolia transcriptome-wide m6A maps. Each library was prepared with three biological replicates. Pearson correlation coefficient analysis among biological replicates showed reliable repeatability (Figure S1 in Supplementary Materials). As shown in Table S1, we generated a total of 23–28 million reads for each m6A-seq sample and 21–32 million reads for each input sample (Table S1). The proportion of clean and mapped reads in m6A-seq were around 52–77%. Transcriptome-wide m6A modification sites were identified using exomePeak2. After m6A peak calling analysis, we identified 19,783 and 19,609 common peaks in M. prunifolia under control and drought stress conditions, respectively (Table S2 and Table S3 in Supplementary Materials). The proportions of common peaks in M. prunifolia under control conditions were above 80% (84.23%, 85.64%, and 83.04%) (Figure 1a), while the proportions of common peaks in M. prunifolia under drought conditions were all around 73% (72.97%, 72.20%, and 73.81%) (Figure 1b). In order to estimate the accuracy of peaks, we randomly selected five m6A-containing genes from peaks of M. prunifolia under control and drought stress conditions by checking the read abundance in the Integrative Genomics Viewer (IGV) (Figure 1c). In M. prunifolia under control conditions, m6A modifications on MD00G1007600, MD00G1055400, and MD00G1055500 were modified in the 3′ untranslated region (UTR), the 5′ UTR, and the coding sequence (CDS), respectively. In M. prunifolia under drought conditions, m6A modifications on MD02G1049300 and MD00G1054500 were modified in the CDS and 3′ UTR (Figure 1d). These results indicated that the m6A-seq data are reliable.
To acquire a better understanding of the m6A distribution pattern in M. prunifolia under control conditions and drought stress, we evaluated the distribution of m6A in the whole transcriptome. The transcripts were divided into three non-overlapping regions: 5′ UTR, CDS, and 3′ UTR. As shown in Figure 1a,b, m6A modifications in M. prunifolia under control and drought conditions were mainly enriched in 3′ UTR, followed by CDS, with a small amount of enrichment in the 5′ UTR and intergenic regions (Figure 1a,b). The m6A distribution pattern in M. prunifolia under drought stress has a 3.17% greater proportion in the 3′ UTR than that under control conditions, as well as a 3.31% lower percentage in the CDS. Nevertheless, m6A maintained the same distribution trend in control- and drought-treated M. prunifolia, that is, m6A peaks were mainly distributed in the 3′ UTR, followed by the CDS.
Previous work has demonstrated that multiple different regions of a transcript may undergo m6A modification [28]; we therefore calculated the number of peaks in each transcript. As shown in Figure 1d, about 86% m6A-modified transcripts contained one m6A peak, about 11% contained two m6A peaks, and only a few contained more than three m6A peaks. The trends were almost identical in M. prunifolia under control and drought conditions, similar to those in strawberry and pak choi [28,33].
Furthermore, we used the HOMER software to investigate the m6A modification motif in M. prunifolia under control and drought conditions. Results showed that the UGUAH (H: A/U/C) sequence motif is the predominant and conserved sequence in M. prunifolia under control and drought conditions (Figure 1a,b). In the HOMER results, the motif containing UGUAH is ranked first in M. prunifolia under control and drought conditions with a p-value of 1e-245 and 1e-232, respectively (Figure 1a,b). The UGUAH motif was consistent with findings in Arabidopsis [24], tomato [28], and maize [32]. However, the conserved m6A modification motif RRACH (R: A/G; H: A/U/C) was also discovered in Arabidopsis [25,60,61], indicating that the pattern of m6A modification varies among species.
Since m6A methylation has been widely reported to be involved in regulating biological processes in plants [24,25,28,31,33,62], we performed gene ontology (GO) enrichment analyses of m6A-containing genes in control- and drought-treated M. prunifolia seedlings. As shown in Figure 2a,b, the m6A-containing genes in M. prunifolia under control and drought were significantly enriched in a number of pathways: (1) RNA processing: RNA processing, RNA splicing, ncRNA (metabolic) processing, and tRNA (metabolic) processing; (2) others: RNA 3′-end processing, mRNA transport, histone modification, regulation of gene expression, chromosome (chromatin) organization, and DNA methylation or demethylation; (3) development: flower development, fruit development, and post-embryonic development; (4) stress: response to abiotic stimulus, response to osmotic stress, response to heat, response to temperature stimulus, response to salt stress, response to stimulus, response to abiotic stimulus, response to metal ion, immune response, protein folding, and fatty acid metabolic process. Notably, the m6A-containing genes under drought conditions showed a higher proportion in response to abiotic stimulus (GO: 0009628) than those under control conditions in M. prunifolia (Figure 2c). These results indicated not only that m6A is widely involved in various biological processes in M. prunifolia but implied that m6A modification in M. prunifolia is responsive to abiotic stimulus after drought treatment.

3.2. Differential m6A Methylation between Control and Drought-Treated M. prunifolia Seedlings

To observe the changes of m6A methylation in M. prunifolia seedlings after drought treatment, we produced m6A distribution plots (Figure 3a) using GuitarPlot [63] and a histogram of UGUAH statistics in peaks (Figure 3b. The transcript features in m6A distribution plots included upstream 1 kb, 5′ UTR, CDS, 3′ UTR, and downstream 1 kb. As shown in Figure 3a, the densities of m6A peaks were mainly enriched in the 3′ UTR in M. prunifolia under control and drought treatment, though we noticed a slightly increased peak density in M. prunifolia under drought conditions (Figure 3a). After calculating the percentage of the subsequences of the UGUAH motif (UGUAU, UGUAA, and UGUAC) under control and drought conditions, we found that UGUAA and UGUAU were mainly enriched in M. prunifolia. Compared to control-treated M. prunifolia, the proportion of all three subsequences of UGUAH in drought-treated M. prunifolia showed a slight increase. Our results suggested that the modification patterns of m6A did not change in M. prunifolia after the drought treatment.
To gain better insight into the potential roles of m6A in regulating drought resistance in M. prunifolia, we next focused on the differential m6A peaks with thresholds of |log2 (fold change) | > 0.5 and adjusted the p-value < 0.05 by comparing the m6A methylome of M. prunifolia under drought and control conditions (Table S4 in Supplementary Materials). After drought treatment, 4069 peaks were up-regulated and 135 peaks were down-regulated, corresponding to 4026 and 135 transcripts, showing that more peaks were hypermethylated under drought stress in M. prunifolia (Figure 3c). We then assigned the differential m6A peaks to transcript features (Figure 3d). As shown in Figure 3d, the 4069 hypermethylated m6A peaks were highly enriched in the 3′ UTR (71.88%) and CDS (24.70%); similarly, the 135 hypomethylated peaks were also mainly distributed around the 3′ UTR (70.37%) and CDS (19.26%). These results are consistent with the increased m6A peak density in 3′ UTR in M. prunifolia under drought conditions (Figure 3a). To further explore the biofunctional aspects of these hypermethylated and hypomethylated genes, GO enrichment analysis was performed. For the hypermethylated genes, most genes were significantly enriched in stress and stimulus-related GO slims, such as response to abiotic stimulus, abscisic acid (ABA), heat, temperature stimulus, osmotic stress, and protein folding. Other GO slims that exhibited an association with epigenetic regulation included RNA processing and splicing, chromatin organization, and gene expression regulation. Additionally, hypomethylated genes were mainly enriched in the lipid metabolic process (GO:0006629), response to abiotic stimulus (GO:0009628), response to abscisic acid (GO:0009737), response to oxygen-containing compounds (GO:1901700), response to light stimulus (GO:0009416), and response to stimulus (GO:0050896). These data imply that the m6A levels of some genes in response to drought, including those responsive to ABA, may be influenced by the drought treatment.

3.3. Differential Gene Expression Analysis

M. prunifolia is known for its tolerance of harsh environments and is particularity adapted to drought stress [38]. To profile the gene expression changes regulated by drought stress in M. prunifolia, we performed differential gene expression analysis using RNA-seq data. The differentially expressed genes (DEGs) were identified with thresholds of |log2 (fold change)|> 1 and an adjusted p-value < 0.05 by comparing the reads of M. prunifolia under control and drought conditions using the DESeq2 package (Table S5 in Supplementary Materials) [44]. As shown in Figure 4a, the volcano plot showed that 6029 genes were up-regulated in M. prunifolia under drought stress compared with the control condition, while 8034 genes were down-regulated. The heatmap also displayed the same results using fragments per kilobase of exon model per million mapped reads (FPKM). GO enrichment analysis revealed that DEGs significantly concentrated in relation to three aspects of GO slims: (1) metabolic processes: the positive flavonoid metabolic process and fatty acid metabolic process; (2) hormones: the hormone-mediated signaling pathway, response to abscisic acid, and response to hormone; and (3) stress: (positive regulation of) response to stimulus, response to abiotic stimulus, response to osmotic stress, response to water (deprivation), response to oxidative stress, and immune response. These data showed that M. prunifolia underwent dramatic and significant changes in the expression levels of a large number of drought-related genes after drought stress, which may be related to the drought resistance of M. prunifolia.

3.4. Association Analysis of m6A Levels with Gene Expressions Involved in Apple Drought Tolerance

As a common regulatory mechanism, m6A modification regulates gene expression in a wide range of biological processes [64]. In order to estimate the relationship between the m6A modification and gene expression levels, we divided the genes into nine groups according to FPKM from low to high as well as into three categories based on transcript distribution, and calculated the fraction of m6A-containing genes in each group (Figure 5a). As shown in Figure 5a, the m6A peak fraction increased in the 5′ UTR, 3′ UTR, and CDS with increasing gene expression levels. The highest density was at the seventh group for 3′ UTR, the eight group for CDS, and the ninth group for 5′ UTR. Overall, the m6A peak fraction between control and drought treatment showed little variation in the CDS and 5′ UTR, but more m6A modifications were correlated with higher gene expression levels in 3′ UTR under drought conditions.
To investigate the potential relationship between m6A levels and gene expression levels in response to drought stress, we performed association analyses between differential m6A peaks and DEGs. As shown in Figure 5b, two volcano diagrams showed the overlapping of methylated genes and DEGs. In hypermethylated genes, 972 genes were up-regulated and 1238 genes were down-regulated in M. prunifolia under drought stress. In hypomethylated genes, 42 and 30 genes showed higher and lower expression levels in M. prunifolia under drought stress, respectively. We then extended this analysis to the entire transcriptome of all m6A-modified genes (Figure 5c,d). Genes bearing hypermethylation and hypomethylation exhibited no significant expression changes compared with the non-differential m6A-modified genes according to a Wilcoxon test (Figure 5d). Considering the transcript distribution characteristics in m6A methylation, we analyzed the gene expression changes in the whole transcriptome (Figure 5e). However, statistical analysis indicated that different transcript distributions did not significantly affect differential gene expression compared with non-differential m6A genes (Figure 5e). These data imply a complex relationship between m6A levels and expression levels.
To further understand the DEGs affected by m6A modifications, we performed GO enrichment analysis. The results revealed that these genes were significantly enriched in stress-related GO slims, such as lipid and pigment biosynthetic processes, the fatty acid metabolic process, protein folding, response to hormone, response to heat, response to osmotic stress, and response to temperature stimuli (Figure 5f). Chromatin organization, which plays a role in plant responses to drought, was also significantly enriched [65]. Based on this GO enrichment analysis, several DEGs affected by m6A modifications were exhibited in IGV [59]. Two up-regulated Heat shock protein 60 (HSP60) genes (MD05G1182500, MD10G1170700) were hypermethylated in M. prunifolia under drought stress (Figure 5h). HSPs are essential components of thermotolerance in plants and HSP60 is reported to be up-regulated in Arabidopsis after high temperature stress [65,66]. Two hypomethylated jasmonate-Zim-domain protein 3 (JAZ3) genes (MD14G1238100, MD16G1020800) were down-regulated in M. prunifolia under drought stress, consistent with results in poplar [67].
It has been reported that the establishment of stress acclimation and stress adaptation associates with changes in genome DNA methylation and may depend on small RNA pathways requiring Dicer-like 2 (DCL2) and DCL3 [68,69]. We also found down-regulated DCL3 with hypomethylation in M. prunifolia under drought stress. The Arabidopsis ETHYLENE RESPONSE FACTOR1 (ERF1)-overexpressing plants (35S:ERF1) are more tolerant to drought and salt stress compared with wild-type plants [70]. Under drought stress, ERF1 (MD13G1135000) in M. prunifolia showed hypermethylation and up-regulated gene expression. Additionally, as shown in Figure 5g, hypermethylated Scarecrow-Like 1 (SCL1) exhibited up-regulated gene expression in M. prunifolia after drought stress. Ct-SCL1 has been reported to play a key role in the survival of the cluster bean (Cyamopsis tetragonaloba L.) under drought stress by interacting with the SWITCH SUBUNIT 3B (SWI3B) protein through chromatin remodeling and stress-based epigenetic memory [71]. The above data suggested that m6A modifications are widely involved in the response to drought in M. prunifolia by affecting the expression of drought-related genes.

4. Discussion

The m6A modification pattern has been reported in many plant species. In different ripening stages of strawberry, m6A peaks are mainly modified in the 3′ UTR and stop codon, followed by the CDS, and there was no significant change in the ratio of the three regions [28]. In maize, m6A peaks frequently occur in the 3′ UTR (around 70%), followed by about 20% in the stop codon [32,72]. This trend is also observed in Arabidopsis [25,73]. During tomato fruit ripening, m6A peaks are mainly modified in the 3′ UTR and stop codon, and the percentage of 3′ UTR modifications increases with fruit ripening [27]. In pak choi (Brassica rapa ssp. chinensis), m6A peaks are mainly enriched in the 3′ UTR and CDS with approximate proportions which show no change after heat stress [33]. Therefore, m6A frequently occurs in th 3′ UTR or stop codon. In contrast, a recent report on apple shows m6A peaks mainly enriching the CDS (53.20%), followed by the 3′ UTR (36.80%) and 5′ UTR (10.00%) [34]. The transcript distribution annotations of peaks are based on the apple gene annotation file without region information for start and stop codons (GDDH13 version 1.1, https://iris.angers.inra.fr/gddh13/downloads/gene_models_20170612.gff3.bz2, accessed on 18 November 2021) [42]. Here, we extracted the region information for the 3′ UTR, CDS, and 5′ UTR from the gene annotation file and completed the transcript annotation using middle point of peaks (probably the strongest point of m6A modifications) [32,74]. Peaks that do not overlap at all with the three regions are defined as intergenic peaks. Our results showed that m6A peaks were mainly enriched in the 3′ UTR in M. prunifolia under control and drought conditions. The difference between the two results in apple may be caused by different annotation methods. As shown in Figure 5g, m6A modification of MD09G1253000 in the 3′ UTR might also occur in the stop codon region. Interestingly, we noticed that only about 60% of the transcripts had 3′ UTR and about 36% of the transcripts did not have either 5′ UTR and 3′ UTR information in the gene annotation file during the process of peak annotation. These data further indicate that incomplete information in the apple gene annotation file and different annotation methods may lead to different m6A distribution results in apple.
m6A is widely conserved among eukaryotes and tends to occur in the RRACH (R: A/G; H: A/U/C) consensus motif, which was identified in mammals in the 1970s [75]. In the early days of m6A research in Arabidopsis, RRACH was also identified [61]. However, a later study in Arabidopsis identified the conserved UGUAY (Y: C/U) motif [23]. In addition, in maize and tomato, the conserved m6A motif was identified as UGUAMM (M = A or C) and UGUAYY, respectively [27,32]. In our results, the UGUAH motif ranked first in the HOMER [56] results and was highly significant, similar to the “URUAY” motif reported by Guo et al. [34]. Moreover, our results also showed that UGUAU and UGUAA had the highest percentage (Figure 3b). The m6A conserved motif in M. prunifolia is relatively consistent with that in Arabidopsis, maize, and tomato. These data indicate not only that UGUA is the core sequence of m6A modification in plants but also show the complexity and bias of m6A modification among different species.
There have been several reports on m6A modifications and their effects on mRNA abundance in plants in response to biotic and abiotic stresses. In pak choi, more m6A modifications show hypermethylation after heat treatment, but m6A affected the same number of up- and down-regulated genes [33]. In apple, MhYTP2 overexpression enhanced apple powdery mildew (PM) resistance and triggered more hypermethylated peaks [34]. Similarly, our results showed that a large number of hypermethylated peaks appeared in M. prunifolia under drought stress. To further investigate the potential relationship between differential m6A peaks and DEGs, we considered the role of methylation types and transcript distribution in m6A modifications and performed an analysis mentioned only in a study of strawberry fruit ripening [28]. However, no significant or apparent correlation was found between the differential deposition of m6A in gene features and altered gene expression. We speculated that this may be related to the different organ or biological process in our study.
Abiotic stresses adversely affect plant growth and productivity. Numerous studies have been conducted to decipher the genetic and molecular mechanisms of plant drought stress tolerance [76,77]. Members of the heat shock protein (HSP) family remodel proteins and play various positive roles in plant responses to drought stress. HSP20 was up-regulated under high-temperature stress in pepper and grasses [78,79]. In tobacco, NtHSP70-1 was found to be an ABA-inducible gene, and the over-expressed NtHSP70-1 can confer drought stress tolerance [80]. In Arabidopsis thaliana, the overexpression of AtHsp90.2, AtHsp90.5, and AtHsp90.7 enhanced plant responses to drought and salt stresses, and cytosolic Hsp90 might be involved in plant stress responses in an ABA-dependent manner [81]. In our results for m6A-modified genes with expression changes in drought-treated M. prunifolia compared to control-treated M. prunifolia, we found a number of up-regulated HSP genes encoding HSP20, HSP70, HSP90.5, HSP88.1, HSP90.6, and HSP60. Further, it is well known that ABA mediates the drought stress response by regulating stomatal closure and stress-responsive gene expression [82,83]. SNF1-ralated protein kinases 2 (SnRK2s) are key regulators that manage the adaptive responses to osmotic stress, including drought stress [84]. SNRK2.6, a member of subclass III, plays essential roles in the positive regulation of ABA signaling and is strongly activated by ABA [85]. ERF1 integrates JA, ET, and ABA signaling through stress-specific gene regulation and plays a positive role in salt, drought and heat stresses [70]. In Arabidopsis, the overexpression of ERF1 enhanced drought and salt tolerance [70]. Similar to our results, hyper-methylated SNRK2.6 and ERF1 were up-regulated in M. prunifolia after drought treatment as compared to control-treated M. prunifolia. These m6A-modified genes with altered gene expression in HSP encoding genes and the ABA pathway of M. prunifolia after drought treatment establish a link between m6A modification and drought tolerance in apple.
In summary, we depicted a transcriptome-wide m6A profiling in M. prunifolia and investigated changes in m6A modification after drought stress. We found that drought stimulated hypermethylated m6A peaks in M. prunifolia. Several m6A-modified drought-responsive genes, including HSP60, JAZ3, SCL1, and ERF1, were presented. Our research provides new support for understanding the epigenetic regulatory mechanisms of apples in response to drought stress and the breeding of drought-tolerant apple trees.

Supplementary Materials

The following supporting information can be downloaded at: https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/plants11010103/s1, Figure S1: The Pearson correlation analysis of m6A-seq and RNA-seq data, Table S1: Summary of sequenced and mapped reads in m6A-seq and RNA-seq samples generated in this study, Table S2: m6A-modified genes in M. prunifolia seedlings under control conditions, Table S3: m6A-modified genes in M. prunifolia seedlings under drought conditions, Table S4: Differential m6A peaks in drought-compared to control-treated M. prunifolia, Table S5: Differentially expressed genes in drought-versus control-treated M. prunifolia.

Author Contributions

J.H. and Z.L. designed the study. N.H. prepared the plants. X.M. and J.H. analyzed the data. J.H. and X.M. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the “Central government-guided local science and technology development special project—Integrated research and demonstration of apple quality improvement and efficiency technology in Qingyang National Agricultural Science and Technology Park”.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The RNA-seq and m6A-seq data have been deposited with the NCBI with the dataset identifier PRJNA781274.

Acknowledgments

We thank Chuang Ma’s lab (College of Life Sciences, Northwest A&F University) for assistance with the m6A data analysis. We thank the High-Performance Computing (HPC) platform of Northwest A&F University (NWAFU) for providing computing resources. We also thank Novogene (https://www.novogene.com/) for assistance with the m6A-seq assay.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

  1. Boccaletto, P.; Machnicka, M.A.; Purta, E.; Piątkowski, P.; Bagiński, B.; Wirecki, T.K.; de Crécy-Lagard, V.; Ross, R.; Limbach, P.A.; Kotter, A.; et al. MODOMICS: A Database of RNA Modification Pathways. 2017 Update. Nucleic Acids Res. 2017, 46, D303–D307. [Google Scholar] [CrossRef] [PubMed]
  2. Wei, C.-M.; Gershowitz, A.; Moss, B. Methylated Nucleotides Block 5′ Terminus of HeLa Cell Messenger RNA. Cell 1975, 4, 379–386. [Google Scholar] [CrossRef]
  3. Wei, W.; Ji, X.; Guo, X.; Ji, S. Regulatory Role of N6-Methyladenosine (M6A) Methylation in RNA Processing and Human Diseases. J. Cell. Biochem. 2017, 118, 2534–2543. [Google Scholar] [CrossRef]
  4. Cantara, W.A.; Crain, P.F.; Rozenski, J.; McCloskey, J.A.; Harris, K.A.; Zhang, X.; Vendeix, F.A.P.; Fabris, D.; Agris, P.F. The RNA Modification Database, RNAMDB: 2011 Update. Nucleic Acids Res. 2011, 39, D195–D201. [Google Scholar] [CrossRef] [Green Version]
  5. Yue, H.; Nie, X.; Yan, Z.; Weining, S. N6-methyladenosine Regulatory Machinery in Plants: Composition, Function and Evolution. Plant Biotechnol. J. 2019, 17, 1194–1208. [Google Scholar] [CrossRef]
  6. Kierzek, E. The Thermodynamic Stability of RNA Duplexes and Hairpins Containing N6-Alkyladenosines and 2-Methylthio-N6-Alkyladenosines. Nucleic Acids Res. 2003, 31, 4472–4480. [Google Scholar] [CrossRef]
  7. Desrosiers, R.; Friderici, K.; Rottman, F. Identification of Methylated Nucleosides in Messenger RNA from Novikoff Hepatoma Cells. Proc. Natl. Acad. Sci. USA 1974, 71, 3971–3975. [Google Scholar] [CrossRef] [Green Version]
  8. Nichols, J.L.; Welder, L. Nucleotides Adjacent to N6-Methyladenosine in Maize Poly(A)-Containing RNA. Plant Sci. Lett. 1981, 21, 75–81. [Google Scholar] [CrossRef]
  9. Kennedy, T.D.; Lane, B.G. Wheat Embryo Ribonucleates. XIII. Methyl-Substituted Nucleoside Constituents and 5′-Terminal Dinucleotide Sequences in Bulk Poly(A)-Rich RNA from Imbibing Wheat Embryos. Can. J. Biochem. 1979, 57, 927–931. [Google Scholar] [CrossRef]
  10. Haugland, R.A.; Cline, M.G. Post-Transcriptional Modifications of Oat Coleoptile Ribonucleic Acids. 5′-Terminal Capping and Methylation of Internal Nucleosides in Poly(A)-Rich RNA. Eur. J. Biochem. 1980, 104, 271–277. [Google Scholar] [CrossRef]
  11. Saneyoshi, M.; Harada, F.; Nishimura, S. Isolation and Characterization of N6-Methyladenosine from Escherichia Coli Valine Transfer RNA. Biochim. Et Biophys. Acta (BBA) Nucleic Acids Protein Synth. 1969, 190, 264–273. [Google Scholar] [CrossRef]
  12. Meyer, K.D.; Jaffrey, S.R. Rethinking m6A Readers, Writers, and Erasers. Annu. Rev. Cell Dev. Biol. 2017, 33, 319–342. [Google Scholar] [CrossRef] [Green Version]
  13. Zheng, G.; Dahl, J.A.; Niu, Y.; Fedorcsak, P.; Huang, C.-M.; Li, C.J.; Vågbø, C.B.; Shi, Y.; Wang, W.-L.; Song, S.-H.; et al. ALKBH5 Is a Mammalian RNA Demethylase That Impacts RNA Metabolism and Mouse Fertility. Mol. Cell 2013, 49, 18–29. [Google Scholar] [CrossRef] [Green Version]
  14. Ping, X.-L.; Sun, B.-F.; Wang, L.; Xiao, W.; Yang, X.; Wang, W.-J.; Adhikari, S.; Shi, Y.; Lv, Y.; Chen, Y.-S.; et al. Mammalian WTAP Is a Regulatory Subunit of the RNA N6-Methyladenosine Methyltransferase. Cell Res. 2014, 24, 177–189. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Jia, G.; Fu, Y.; Zhao, X.; Dai, Q.; Zheng, G.; Yang, Y.; Yi, C.; Lindahl, T.; Pan, T.; Yang, Y.-G.; et al. N6-Methyladenosine in Nuclear RNA Is a Major Substrate of the Obesity-Associated FTO. Nat. Chem. Biol. 2011, 7, 885–887. [Google Scholar] [CrossRef]
  16. Patil, D.P.; Pickering, B.F.; Jaffrey, S.R. Reading M6A in the Transcriptome: M6A-Binding Proteins. Trends Cell Biol. 2018, 28, 113–127. [Google Scholar] [CrossRef]
  17. Meyer, K.D.; Patil, D.P.; Zhou, J.; Zinoviev, A.; Skabkin, M.A.; Elemento, O.; Pestova, T.V.; Qian, S.-B.; Jaffrey, S.R. 5′ UTR M6A Promotes Cap-Independent Translation. Cell 2015, 163, 999–1010. [Google Scholar] [CrossRef] [Green Version]
  18. Zhao, X.; Yang, Y.; Sun, B.-F.; Shi, Y.; Yang, X.; Xiao, W.; Hao, Y.-J.; Ping, X.-L.; Chen, Y.-S.; Wang, W.-J.; et al. FTO-Dependent Demethylation of N6-Methyladenosine Regulates MRNA Splicing and Is Required for Adipogenesis. Cell Res. 2014, 24, 1403–1419. [Google Scholar] [CrossRef]
  19. Wang, X.; Lu, Z.; Gomez, A.; Hon, G.C.; Yue, Y.; Han, D.; Fu, Y.; Parisien, M.; Dai, Q.; Jia, G.; et al. N6-Methyladenosine-Dependent Regulation of Messenger RNA Stability. Nature 2014, 505, 117–120. [Google Scholar] [CrossRef]
  20. Xiao, W.; Adhikari, S.; Dahal, U.; Chen, Y.-S.; Hao, Y.-J.; Sun, B.-F.; Sun, H.-Y.; Li, A.; Ping, X.-L.; Lai, W.-Y.; et al. Nuclear m6A Reader YTHDC1 Regulates MRNA Splicing. Mol. Cell 2016, 61, 507–519. [Google Scholar] [CrossRef] [Green Version]
  21. Kasowitz, S.D.; Ma, J.; Anderson, S.J.; Leu, N.A.; Xu, Y.; Gregory, B.D.; Schultz, R.M.; Wang, P.J. Nuclear M6A Reader YTHDC1 Regulates Alternative Polyadenylation and Splicing during Mouse Oocyte Development. PLoS Genet. 2018, 14, e1007412. [Google Scholar] [CrossRef]
  22. Fustin, J.-M.; Doi, M.; Yamaguchi, Y.; Hida, H.; Nishimura, S.; Yoshida, M.; Isagawa, T.; Morioka, M.S.; Kakeya, H.; Manabe, I.; et al. RNA-Methylation-Dependent RNA Processing Controls the Speed of the Circadian Clock. Cell 2013, 155, 793–806. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Arribas-Hernández, L.; Bressendorff, S.; Hansen, M.H.; Poulsen, C.; Erdmann, S.; Brodersen, P. An M6A-YTH Module Controls Developmental Timing and Morphogenesis in Arabidopsis. Plant Cell 2018, 30, 952–967. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Wei, L.-H.; Song, P.; Wang, Y.; Lu, Z.; Tang, Q.; Yu, Q.; Xiao, Y.; Zhang, X.; Duan, H.-C.; Jia, G. The m6A Reader ECT2 Controls Trichome Morphology by Affecting MRNA Stability in Arabidopsis. Plant Cell 2018, 30, 968–985. [Google Scholar] [CrossRef] [Green Version]
  25. Duan, H.-C.; Wei, L.-H.; Zhang, C.; Wang, Y.; Chen, L.; Lu, Z.; Chen, P.R.; He, C.; Jia, G. ALKBH10B Is an RNA N6-Methyladenosine Demethylase Affecting Arabidopsis Floral Transition. Plant Cell 2017, 29, 2995–3011. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Zhang, F.; Zhang, Y.-C.; Liao, J.-Y.; Yu, Y.; Zhou, Y.-F.; Feng, Y.-Z.; Yang, Y.-W.; Lei, M.-Q.; Bai, M.; Wu, H.; et al. The Subunit of RNA N6-Methyladenosine Methyltransferase OsFIP Regulates Early Degeneration of Microspores in Rice. PLoS Genet. 2019, 15, 1–19. [Google Scholar] [CrossRef] [Green Version]
  27. Zhou, L.; Tian, S.; Qin, G. RNA Methylomes Reveal the M6A-Mediated Regulation of DNA Demethylase Gene SlDML2 in Tomato Fruit Ripening. Genome Biol. 2019, 20, 156. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Zhou, L.; Tang, R.; Li, X.; Tian, S.; Li, B.; Qin, G. N6-Methyladenosine RNA Modification Regulates Strawberry Fruit Ripening in an ABA-Dependent Manner. Genome Biol. 2021, 22, 168. [Google Scholar] [CrossRef] [PubMed]
  29. Ok, S.H.; Jeong, H.J.; Bae, J.M.; Shin, J.-S.; Luan, S.; Kim, K.-N. Novel CIPK1-Associated Proteins in Arabidopsis Contain an Evolutionarily Conserved C-Terminal Region That Mediates Nuclear Localization. Plant Physiol. 2005, 139, 138–150. [Google Scholar] [CrossRef] [Green Version]
  30. Huong, T.T.; Ngoc, L.N.T.; Kang, H. Functional Characterization of a Putative RNA Demethylase ALKBH6 in Arabidopsis Growth and Abiotic Stress Responses. Int. J. Mol. Sci. 2020, 21, 6707. [Google Scholar] [CrossRef]
  31. Zhang, K.; Zhuang, X.; Dong, Z.; Xu, K.; Chen, X.; Liu, F.; He, Z. The Dynamics of N6-Methyladenine RNA Modification in Interactions between Rice and Plant Viruses. Genome Biol. 2021, 22, 189. [Google Scholar] [CrossRef]
  32. Miao, Z.; Zhang, T.; Qi, Y.; Song, J.; Han, Z.; Ma, C. Evolution of the RNA N6-Methyladenosine Methylome Mediated by Genomic Duplication. Plant Physiol. 2020, 182, 345–360. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Liu, G.; Wang, J.; Hou, X. Transcriptome-Wide N6-Methyladenosine (M6A) Methylome Profiling of Heat Stress in Pak-Choi (Brassica Rapa ssp. Chinensis). Plants 2020, 9, 1080. [Google Scholar] [CrossRef] [PubMed]
  34. Guo, T.; Liu, C.; Meng, F.; Hu, L.; Fu, X.; Yang, Z.; Wang, N.; Jiang, Q.; Ma, F. The M6A Reader MhYTP2 Regulates MdMLO19 MRNA Stability and Antioxidant Genes Translation Efficiency Conferring Powdery Mildew Resistance in Apple. Plant Biotechnol. J. 2021. [Google Scholar] [CrossRef]
  35. Nawaz, H.; Hussain, N.; Jamil, M.; Yasmeen, A.; Bukhari, S.A.H.; Aurangzaib, M.; Usman, M. Seed Biopriming Mitigates Terminal Drought Stress at Reproductive Stage of Maize by Enhancing Gas Exchange Attributes and Nutrient Uptake. Turk. J. Agric. For. 2020, 44, 250–261. [Google Scholar] [CrossRef]
  36. Li, X.; Chen, P.; Xie, Y.; Yan, Y.; Wang, L.; Dang, H.; Zhang, J.; Xu, L.; Ma, F.; Guan, Q. Apple SERRATE Negatively Mediates Drought Resistance by Regulating MdMYB88 and MdMYB124 and MicroRNA Biogenesis. Hortic. Res. 2020, 7, 98. [Google Scholar] [CrossRef]
  37. Radivojevic, D.; Milivojevic, J.; Pavlovic, M.; Stopar, M. Comparison of Metamitron Efficiency for Postbloom Thinning of Young ‘Gala’ and ‘Golden Delicious’ Apple Trees. Turk. J. Agric. For. 2020, 44, 83–94. [Google Scholar] [CrossRef]
  38. Tan, Y.; Li, M.; Yang, Y.; Sun, X.; Wang, N.; Liang, B.; Ma, F. Overexpression of MpCYS4, A Phytocystatin Gene from Malus Prunifolia (Willd.) Borkh., Enhances Stomatal Closure to Confer Drought Tolerance in Transgenic Arabidopsis and Apple. Front. Plant Sci. 2017, 8, 33. [Google Scholar] [CrossRef] [Green Version]
  39. Sokolov, V.V.; Savel’ev, N.I.; Goncharov, N.P.I.V. Michurin’S Work on Expansion of the Plant Horticulture Assortment and Improvement of Food Quality. Proc. Latv. Acad. Sci. Sect. B. Nat. Exact Appl. Sci. 2015, 69, 190–197. [Google Scholar] [CrossRef] [Green Version]
  40. Chen, P.; Yan, M.; Li, L.; He, J.; Zhou, S.; Li, Z.; Niu, C.; Bao, C.; Zhi, F.; Ma, F.; et al. The Apple DNA-Binding One Zinc-Finger Protein MdDof54 Promotes Drought Resistance. Hortic. Res. 2020, 7, 195. [Google Scholar] [CrossRef]
  41. Xie, Y.; Chen, P.; Yan, Y.; Bao, C.; Li, X.; Wang, L.; Shen, X.; Li, H.; Liu, X.; Niu, C.; et al. An Atypical R2R3 MYB Transcription Factor Increases Cold Hardiness by CBF-Dependent and CBF-Independent Pathways in Apple. New Phytol. 2018, 218, 201–218. [Google Scholar] [CrossRef] [PubMed]
  42. Daccord, N.; Celton, J.-M.; Linsmith, G.; Becker, C.; Choisne, N.; Schijlen, E.; van de Geest, H.; Bianco, L.; Micheletti, D.; Velasco, R.; et al. High-Quality de Novo Assembly of the Apple Genome and Methylome Dynamics of Early Fruit Development. Nat. Genet. 2017, 49, 1099–1106. [Google Scholar] [CrossRef]
  43. Anders, S.; Pyl, P.T.; Huber, W. HTSeq--a Python Framework to Work with High-Throughput Sequencing Data. Bioinformatics 2015, 31, 166–169. [Google Scholar] [CrossRef]
  44. Love, M.I.; Huber, W.; Anders, S. Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  45. Lawrence, M.; Huber, W.; Pagès, H.; Aboyoun, P.; Carlson, M.; Gentleman, R.; Morgan, M.T.; Carey, V.J. Software for Computing and Annotating Genomic Ranges. PLoS Comput. Biol. 2013, 9, e1003118. [Google Scholar] [CrossRef] [PubMed]
  46. Chen, C.; Chen, H.; Zhang, Y.; Thomas, H.R.; Frank, M.H.; He, Y.; Xia, R. TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data. Molecular Plant 2020, 13, 1194–1202. [Google Scholar] [CrossRef] [PubMed]
  47. Kolde, R. Pheatmap: Pretty Heatmaps. R Package Version 2012, 1, 726. [Google Scholar]
  48. Tian, T.; Liu, Y.; Yan, H.; You, Q.; Yi, X.; Du, Z.; Xu, W.; Su, Z. AgriGO v2.0: A GO Analysis Toolkit for the Agricultural Community, 2017 Update. Nucleic Acids Res. 2017, 45, W122–W129. [Google Scholar] [CrossRef]
  49. Wu, T.; Hu, E.; Xu, S.; Chen, M.; Guo, P.; Dai, Z.; Feng, T.; Zhou, L.; Tang, W.; Zhan, L.; et al. ClusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. Innovation 2021, 2, 100141. [Google Scholar] [CrossRef]
  50. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A Flexible Trimmer for Illumina Sequence Data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef] [Green Version]
  51. Fiancette, R.; Finlay, C.M.; Willis, C.; Bevington, S.L.; Soley, J.; Ng, S.T.H.; Baker, S.M.; Andrews, S.; Hepworth, M.R.; Withers, D.R. Reciprocal Transcription Factor Networks Govern Tissue-Resident ILC3 Subset Function and Identity. Nat. Immunol. 2021, 22, 1245–1255. [Google Scholar] [CrossRef]
  52. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Danecek, P.; Bonfield, J.K.; Liddle, J.; Marshall, J.; Ohan, V.; Pollard, M.O.; Whitwham, A.; Keane, T.; McCarthy, S.A.; Davies, R.M.; et al. Twelve Years of SAMtools and BCFtools. GigaScience 2021, 10, giab008. [Google Scholar] [CrossRef] [PubMed]
  54. Liu, L.; Song, B.; Ma, J.; Song, Y.; Zhang, S.-Y.; Tang, Y.; Wu, X.; Wei, Z.; Chen, K.; Su, J.; et al. Bioinformatics Approaches for Deciphering the Epitranscriptome: Recent Progress and Emerging Topics. Comput. Struct. Biotechnol. J. 2020, 18, 1587–1604. [Google Scholar] [CrossRef]
  55. Khan, A.; Mathelier, A. Intervene: A Tool for Intersection and Visualization of Multiple Gene or Genomic Region Sets. BMC Bioinform. 2017, 18, 287. [Google Scholar] [CrossRef] [Green Version]
  56. Heinz, S.; Benner, C.; Spann, N.; Bertolino, E.; Lin, Y.C.; Laslo, P.; Cheng, J.X.; Murre, C.; Singh, H.; Glass, C.K. Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities. Mol. Cell 2010, 38, 576–589. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Quinlan, A.R.; Hall, I.M. BEDTools: A Flexible Suite of Utilities for Comparing Genomic Features. Bioinformatics 2010, 26, 841–842. [Google Scholar] [CrossRef] [Green Version]
  58. Zhai, J.; Song, J.; Cheng, Q.; Tang, Y.; Ma, C. PEA: An Integrated R Toolkit for Plant Epitranscriptome Analysis. Bioinformatics 2018, 34, 3747–3749. [Google Scholar] [CrossRef] [Green Version]
  59. Robinson, J.T.; Thorvaldsdóttir, H.; Wenger, A.M.; Zehir, A.; Mesirov, J.P. Variant Review with the Integrative Genomics Viewer. Cancer Res. 2017, 77, e31–e34. [Google Scholar] [CrossRef] [Green Version]
  60. Luo, G.-Z.; MacQueen, A.; Zheng, G.; Duan, H.; Dore, L.C.; Lu, Z.; Liu, J.; Chen, K.; Jia, G.; Bergelson, J.; et al. Unique Features of the M6A Methylome in Arabidopsis Thaliana. Nat. Commun. 2014, 5, 5630. [Google Scholar] [CrossRef] [Green Version]
  61. Wang, Z.; Tang, K.; Zhang, D.; Wan, Y.; Wen, Y.; Lu, Q.; Wang, L. High-Throughput M6A-Seq Reveals RNA M6A Methylation Patterns in the Chloroplast and Mitochondria Transcriptomes of Arabidopsis Thaliana. PLoS ONE 2017, 12, e0185612. [Google Scholar] [CrossRef]
  62. Bhat, S.S.; Bielewicz, D.; Gulanicz, T.; Bodi, Z.; Yu, X.; Anderson, S.J.; Szewc, L.; Bajczyk, M.; Dolata, J.; Grzelak, N.; et al. MRNA Adenosine Methylase (MTA) Deposits m6A on Pri-MiRNAs to Modulate MiRNA Biogenesis in Arabidopsis Thaliana. Proc. Natl. Acad. Sci. USA 2020, 117, 21785–21795. [Google Scholar] [CrossRef] [PubMed]
  63. Cui, X.; Wei, Z.; Zhang, L.; Liu, H.; Sun, L.; Zhang, S.-W.; Huang, Y.; Meng, J. 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] [CrossRef] [Green Version]
  64. Zaccara, S.; Ries, R.J.; Jaffrey, S.R. Reading, Writing and Erasing MRNA Methylation. Nat. Rev. Mol. Cell. Biol. 2019, 20, 608–624. [Google Scholar] [CrossRef] [PubMed]
  65. Kim, J.-M.; Sasaki, T.; Ueda, M.; Sako, K.; Seki, M. Chromatin Changes in Response to Drought, Salinity, Heat, and Cold Stresses in Plants. Front. Plant Sci. 2015, 6, 114. [Google Scholar] [CrossRef] [Green Version]
  66. Grigorova, B.; Vaseva, I.; Demirevska, K.; Feller, U. Combined Drought and Heat Stress in Wheat: Changes in Some Heat Shock Proteins. Biologia Plant. 2011, 55, 105–111. [Google Scholar] [CrossRef]
  67. Jia, J.; Zhou, J.; Shi, W.; Cao, X.; Luo, J.; Polle, A.; Luo, Z.-B. Comparative Transcriptomic Analysis Reveals the Roles of Overlapping Heat-/Drought-Responsive Genes in Poplars Exposed to High Temperature and Drought. Sci. Rep. 2017, 7, 43215. [Google Scholar] [CrossRef] [Green Version]
  68. Boyko, A.; Blevins, T.; Yao, Y.; Golubov, A.; Bilichak, A.; Ilnytskyy, Y.; Hollander, J.; Meins, F.; Kovalchuk, I. Transgenerational Adaptation of Arabidopsis to Stress Requires DNA Methylation and the Function of Dicer-Like Proteins. PLoS ONE 2010, 5, e9514. [Google Scholar] [CrossRef]
  69. Yao, Y.; Bilichak, A.; Golubov, A.; Kovalchuk, I. Arabidopsis Thaliana SiRNA Biogenesis Mutants Have the Lower Frequency of Homologous Recombination. Plant Signal. Behav. 2016, 11, e1151599. [Google Scholar] [CrossRef] [Green Version]
  70. Cheng, M.-C.; Liao, P.-M.; Kuo, W.-W.; Lin, T.-P. The Arabidopsis ETHYLENE RESPONSE FACTOR1 Regulates Abiotic Stress-Responsive Gene Expression by Binding to Different Cis-Acting Elements in Response to Different Stress Signals. Plant Physiol. 2013, 162, 1566–1582. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  71. Manohara Reddy, B.; Kumari, V.; Anthony Johnson, A.; Jagadeesh Kumar, N.; Venkatesh, B.; Jayamma, N.; Sudhakar, C. Scarecrow like Protein 1,(Ct-SCL1) Involved in Drought Stress Tolerance by Interacting with SWI3B Component of Chromatin Modelling Complex in Cluster Bean, Cyamopsistetragonaloba (L.) Taub. Int. J. Res. Anal. Rev. 2018, 5. [Google Scholar]
  72. Luo, J.-H.; Wang, Y.; Wang, M.; Zhang, L.-Y.; Peng, H.-R.; Zhou, Y.-Y.; Jia, G.-F.; He, Y. Natural Variation in RNA m6A Methylation and Its Relationship with Translational Status. Plant Physiol. 2020, 182, 332–344. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  73. Shen, L.; Liang, Z.; Gu, X.; Chen, Y.; Teo, Z.W.N.; Hou, X.; Cai, W.M.; Dedon, P.C.; Liu, L.; Yu, H. N6-Methyladenosine RNA Modification Regulates Shoot Stem Cell Fate in Arabidopsis. Dev. Cell 2016, 38, 186–200. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  74. Miao, Z.; Zhang, T.; Xie, B.; Qi, Y.; Ma, C. Evolutionary Implications of the RNA N 6-Methyladenosine Methylome in Plants. Mol. Biol. Evol. 2021, msab299. [Google Scholar] [CrossRef]
  75. Wang, K.; Peng, J.; Yi, C. The m6A Consensus Motif Provides a Paradigm of Epitranscriptomic Studies. Biochemistry 2021, 60, 3410–3412. [Google Scholar] [CrossRef] [PubMed]
  76. Adams, H.D.; Guardiola-Claramonte, M.; Barron-Gafford, G.A.; Villegas, J.C.; Breshears, D.D.; Zou, C.B.; Troch, P.A.; Huxman, T.E. Temperature Sensitivity of Drought-Induced Tree Mortality Portends Increased Regional Die-off under Global-Change-Type Drought. Proc. Natl. Acad. Sci. USA 2009, 106, 7063–7066. [Google Scholar] [CrossRef] [Green Version]
  77. Šircelj, H.; Tausz, M.; Grill, D.; Batič, F. Biochemical Responses in Leaves of Two Apple Tree Cultivars Subjected to Progressing Drought. J. Plant Physiol. 2005, 162, 1308–1318. [Google Scholar] [CrossRef]
  78. Xu, Y.; Zhan, C.; Huang, B. Heat Shock Proteins in Association with Heat Tolerance in Grasses. Int. J. Proteom. 2011, 2011, 529648. [Google Scholar] [CrossRef] [Green Version]
  79. Huang, L.-J.; Cheng, G.-X.; Khan, A.; Wei, A.-M.; Yu, Q.-H.; Yang, S.-B.; Luo, D.-X.; Gong, Z.-H. CaHSP16.4, a Small Heat Shock Protein Gene in Pepper, Is Involved in Heat and Drought Tolerance. Protoplasma 2019, 256, 39–51. [Google Scholar] [CrossRef]
  80. Cho, E.K.; Hong, C.B. Over-Expression of Tobacco NtHSP70-1 Contributes to Drought-Stress Tolerance in Plants. Plant Cell Rep. 2006, 25, 349–358. [Google Scholar] [CrossRef]
  81. Song, H.; Zhao, R.; Fan, P.; Wang, X.; Chen, X.; Li, Y. Overexpression of AtHsp90.2, AtHsp90.5 and AtHsp90.7 in Arabidopsis Thaliana Enhances Plant Sensitivity to Salt and Drought Stresses. Planta 2009, 229, 955–964. [Google Scholar] [CrossRef] [PubMed]
  82. Li, C.; Tan, D.-X.; Liang, D.; Chang, C.; Jia, D.; Ma, F. Melatonin Mediates the Regulation of ABA Metabolism, Free-Radical Scavenging, and Stomatal Behaviour in Two Malus Species under Drought Stress. J. Exp. Bot. 2015, 66, 669–680. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  83. Takahashi, F.; Kuromori, T.; Urano, K.; Yamaguchi-Shinozaki, K.; Shinozaki, K. Drought Stress Responses and Resistance in Plants: From Cellular Responses to Long-Distance Intercellular Communication. Front. Plant Sci. 2020, 11, 556972. [Google Scholar] [CrossRef] [PubMed]
  84. Soma, F. Plant Raf-like Kinases Regulate the MRNA Population Upstream of ABA-Unresponsive SnRK2 Kinases under Drought Stress. Nat. Commun. 2020, 11, 1373. [Google Scholar] [CrossRef] [Green Version]
  85. Mizoguchi, M.; Umezawa, T.; Nakashima, K.; Kidokoro, S.; Takasaki, H.; Fujita, Y.; Yamaguchi-Shinozaki, K.; Shinozaki, K. Two Closely Related Subclass II SnRK2 Protein Kinases Cooperatively Regulate Drought-Inducible Gene Expression. Plant Cell Physiol. 2010, 51, 842–847. [Google Scholar] [CrossRef] [Green Version]
Plants 11 00103 g001 550
Figure 1. Transcriptome-wide m6A methylome in Malus prunifolia seedlings in response to drought. (a,b) Venn diagrams demonstrating the overlap of m6A peaks from three replicates (top), the enriched motif (bottom), and m6A peak distribution along transcripts (middle) under control (a) or drought conditions (b). The asterisks in the two motifs marked the positions that are modified. CDS, coding sequence; UTR, untranslated region; DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control conditions. (c) Integrative Genomics Viewer (IGV) tracks. (d) Percentage of m6A-containing transcripts with m6A peaks in DroMp and CKMp.
Figure 1. Transcriptome-wide m6A methylome in Malus prunifolia seedlings in response to drought. (a,b) Venn diagrams demonstrating the overlap of m6A peaks from three replicates (top), the enriched motif (bottom), and m6A peak distribution along transcripts (middle) under control (a) or drought conditions (b). The asterisks in the two motifs marked the positions that are modified. CDS, coding sequence; UTR, untranslated region; DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control conditions. (c) Integrative Genomics Viewer (IGV) tracks. (d) Percentage of m6A-containing transcripts with m6A peaks in DroMp and CKMp.
Plants 11 00103 g001
Plants 11 00103 g002 550
Figure 2. Gene ontology (GO) enrichment analysis of m6A-modified genes in Malus prunifolia seedlings in response to drought. (a,b) The “biological process” aspect of the GO enrichment analysis was divided into “RNA processing (epigenetic regulation)”, “others (epigenetic regulation)”, “development”, and “stress”. (c) The gene ratio of “response to abiotic stimulus (GO:0009628)” in Malus prunifolia seedlings in response to drought conditions. DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; FDR, false discovery rate.
Figure 2. Gene ontology (GO) enrichment analysis of m6A-modified genes in Malus prunifolia seedlings in response to drought. (a,b) The “biological process” aspect of the GO enrichment analysis was divided into “RNA processing (epigenetic regulation)”, “others (epigenetic regulation)”, “development”, and “stress”. (c) The gene ratio of “response to abiotic stimulus (GO:0009628)” in Malus prunifolia seedlings in response to drought conditions. DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; FDR, false discovery rate.
Plants 11 00103 g002
Plants 11 00103 g003 550
Figure 3. Differential m6A modifications in Malus prunifolia seedlings in response to drought. (a) Density of m6A peak distribution along transcripts of Malus prunifolia seedlings in response to drought conditions. The area circled in green shows the difference in density between the two groups. (b) Percentage of UGUAH (H: A/U/C). (c) Volcano plot showing hypermethylated and hypomethylated m6A peaks. (d) Pie charts showing the m6A peaks distribution within transcripts. (e) GO enrichment analysis of hypermethylated peak-related genes. (f) GO enrichment analysis of hypomethylated peak-related genes. DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; CDS, coding sequence; UTR, untranslated region; FDR, false discovery rate.
Figure 3. Differential m6A modifications in Malus prunifolia seedlings in response to drought. (a) Density of m6A peak distribution along transcripts of Malus prunifolia seedlings in response to drought conditions. The area circled in green shows the difference in density between the two groups. (b) Percentage of UGUAH (H: A/U/C). (c) Volcano plot showing hypermethylated and hypomethylated m6A peaks. (d) Pie charts showing the m6A peaks distribution within transcripts. (e) GO enrichment analysis of hypermethylated peak-related genes. (f) GO enrichment analysis of hypomethylated peak-related genes. DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; CDS, coding sequence; UTR, untranslated region; FDR, false discovery rate.
Plants 11 00103 g003
Plants 11 00103 g004 550
Figure 4. Differential gene expression in M. prunifolia in response to drought treatment. (a) Volcano plot showing up-regulated genes and down-regulated genes in M. prunifolia after drought treatment. (b) Heat map of differentially expressed genes (DEGs). (c) GO enrichment analysis of DEGs. Sig_Up, up-regulated genes; Sig_Down, down-regulated genes; DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; FDR, false discovery rate.
Figure 4. Differential gene expression in M. prunifolia in response to drought treatment. (a) Volcano plot showing up-regulated genes and down-regulated genes in M. prunifolia after drought treatment. (b) Heat map of differentially expressed genes (DEGs). (c) GO enrichment analysis of DEGs. Sig_Up, up-regulated genes; Sig_Down, down-regulated genes; DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; FDR, false discovery rate.
Plants 11 00103 g004
Plants 11 00103 g005 550
Figure 5. Correlation between m6A modification levels and mRNA abundance in M. prunifolia in response to drought treatment. (a) The ratio of m6A peaks in different transcript distributions to total transcripts in each subgroup was divided by the FPKM. (b) Volcano plots displaying the gene expression ratios of hypermethylated and hypomethylated transcripts. (c) Cumulative fraction of mRNA expression changes. (d) Box plot of gene expression ratios in hypermethylated, hypomethylated, and non-differential transcripts. Hypermethylated, all differentially expressed genes (DEGs) with hypermethylation; hypomethylated, DEGs with hypomethylation; non-differential, DEGs without m6A modification. (e) Box plot of gene expression ratios in CDS, 3′ UTR, 5′ UTR, and non-differential transcripts. None, DEGs without m6A modification. (f) GO enrichment analysis of the overlapping genes between DEGs and differentially modified m6A peak-related genes in M. prunifolia after drought treatment. (g,h) IGV tracks showing the m6A read distribution in drought-related genes from (b) and (c). DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; CDS, coding sequence; UTR, untranslated region.
Figure 5. Correlation between m6A modification levels and mRNA abundance in M. prunifolia in response to drought treatment. (a) The ratio of m6A peaks in different transcript distributions to total transcripts in each subgroup was divided by the FPKM. (b) Volcano plots displaying the gene expression ratios of hypermethylated and hypomethylated transcripts. (c) Cumulative fraction of mRNA expression changes. (d) Box plot of gene expression ratios in hypermethylated, hypomethylated, and non-differential transcripts. Hypermethylated, all differentially expressed genes (DEGs) with hypermethylation; hypomethylated, DEGs with hypomethylation; non-differential, DEGs without m6A modification. (e) Box plot of gene expression ratios in CDS, 3′ UTR, 5′ UTR, and non-differential transcripts. None, DEGs without m6A modification. (f) GO enrichment analysis of the overlapping genes between DEGs and differentially modified m6A peak-related genes in M. prunifolia after drought treatment. (g,h) IGV tracks showing the m6A read distribution in drought-related genes from (b) and (c). DroMp, M. prunifolia seedlings were treated with drought stress; CKMp, M. prunifolia seedlings were grown under control condition; CDS, coding sequence; UTR, untranslated region.
Plants 11 00103 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Mao, X.; Hou, N.; Liu, Z.; He, J. Profiling of N6-Methyladenosine (m6A) Modification Landscape in Response to Drought Stress in Apple (Malus prunifolia (Willd.) Borkh). Plants 2022, 11, 103. https://0-doi-org.brum.beds.ac.uk/10.3390/plants11010103

AMA Style

Mao X, Hou N, Liu Z, He J. Profiling of N6-Methyladenosine (m6A) Modification Landscape in Response to Drought Stress in Apple (Malus prunifolia (Willd.) Borkh). Plants. 2022; 11(1):103. https://0-doi-org.brum.beds.ac.uk/10.3390/plants11010103

Chicago/Turabian Style

Mao, Xiushan, Nan Hou, Zhenzhong Liu, and Jieqiang He. 2022. "Profiling of N6-Methyladenosine (m6A) Modification Landscape in Response to Drought Stress in Apple (Malus prunifolia (Willd.) Borkh)" Plants 11, no. 1: 103. https://0-doi-org.brum.beds.ac.uk/10.3390/plants11010103

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop