Transcriptome-Wide N6-Methyladenosine (m6A) Methylome Profiling of Heat Stress in Pak-choi (Brassica rapa ssp. chinensis)

Background: In higher eukaryotes, N6-methyladenosine (m6A) is the most common internal form of messenger RNA modification. By mapping the m6A methyl genome in multiple species, the potential regulatory function of reversible m6A methylation on mRNA is revealed. Recent studies have shown that RNA m6A modification influences mRNA transcription, location, translation, stability, splicing, and nuclear export. However, there are not enough data on the m6A transcriptome-wide map and its potential biological role in the heat stress of Pak-choi (Brassica rapa ssp. chinensis). Methods: In this work, MeRIP-seq was used to obtain the first transcriptome-wide profiling of RNA m6A modification in Pak-choi. Meanwhile, the transcriptome data were obtained by analyzing the input samples’ sequencing data. Results: Our research indicated that with three replicates, there were 11,252 common m6A peaks and 9729 common m6A-containing genes identified in the normal (CK) and heat stress (T43) groups. It was found that m6A peaks were highly enriched in the 3′ untranslated region in both CK and T43 groups. About 80% of the genes have one m6A site. The consensus sequence of m6A peaks was also enriched, which showed as AAACCV (V: U/A/G). In addition, association analysis found that there is a certain relationship between the degree of m6A methylation and the transcription level, indicating that m6A plays a certain regulatory role in gene expression. Conclusion: This comprehensive map in the study may provide a solid basis for determining the potential function of RNA m6A modification in Pak-choi under normal (CK) and heat stress (T43) conditions.


Background
In higher eukaryotes, N 6 -methyladenosine (m 6 A) is the most common form of internal modification in long-noncoding RNAs and polyadenylated mRNAs, which was first detected in the 1970s [1]. Generally speaking, it is catalyzed by a multicomponent complex composed of two active methyltransferases (such as methyltransferase 3 and methyltransferase 14). To date, the most common of the more than 100 types of RNA modifications identified are m 6 A, m 5 C, and m 1 C RNA methylation [2,3]. Studies have shown that defects in m 6 A methylation or demethylation will lead to serious physiological consequences [4]. With continuous research in this area, it has become increasingly clear that m 6 A is essential for gene expression regulation and plays an important role in plant development.
Some previous studies have shown that RNA m 6 A plays an important role in regulating cellular metabolism and controlling the migration of embryonic stem cells in mammals [5]. In addition, RNA m 6 A can also regulate the totipotency of stem cells in mice [6] and the operation of stem cells in Arabidopsis [7]. In Arabidopsis thaliana, it was found that every 1000 nucleotides contained 0.5-0.7 sites, and each actively expressed transcript contained 0.7-1.0 sites [8,9]. In plants, m 6 A is a combination of m 6 A methyltransferase and consensus sequence RRm 6 ACH (R: G/A; H: U > A > C) [7]. In order to gain a deeper understanding of the potential biological function of RNA m 6 A modification, it is necessary to detect m 6 A modification sites within the transcriptome. With the continuous development of science and technology on m 6 A research in humans and mice, a new method has been developed to perform transcriptome-wide m 6 A localization analysis, that is, using next-generation sequencing technology (MeRIP-seq) for methylated RNA immunoprecipitation and analysis of transcriptome-wide m 6 A [10,11]. This new research method (MeRIP-seq) has also been used to obtain the first full heat-stress map in the m 6 A transcriptome range in Pak-choi. Using this method, more than two-thirds of the transcripts in Arabidopsis were detected with m 6 A modification [12]. Recent studies have shown that m 6 A is predominantly located near the stop codons and 3 UTR [12]. These findings also suggest that m 6 A modification is highly dynamic and plays a specific role in regulating plant development.
Pak-choi, which is widely grown in the world today, is one of the most important vegetables in China. However, the gradual increase in ambient temperatures has affected the normal growth and development of crops, leading to a reduction in crop yield and quality [13,14]. In recent years, there have been studies on the heat stress response, including research on heat signal transduction pathways, heat stress protein identification, and transcriptional regulatory factors [15,16]. Previous studies have found that the expression level of heat stress proteins (HSPs) is affected and regulated by m 6 A on RNA [17]. Further research has found that 5 UTR m 6 A at a single site mediates thermal stress-induced translation of HSP70 [10,18]. Currently, there is very little research on transcriptome-wide N 6 -methyladenosine methylome profiling of heat stress in Pak-choi and other plants. In this work, we acquire the first-ever m 6 A transcriptome-wide map of heat stress in Pak-choi. In order to further study the function of m 6 A and provide a basis for identification of m 6 A future research, we collected transcriptomes of normal (CK) and heat stress (T43) condition leaf tissues from Pak-choi. Here, we obtain the first-known m 6 A map of the transcriptome range in Pak-choi. We also compare and analyze the patterns of m 6 A distribution between CK and T43 conditions to obtain differentially methylated peaks and then analyze potential functions in gene expression regulation under high temperature stress.

Transcriptome-Wide Detection of m 6 A Modification in Pak-choi
Using Illumina Novaseq™ 6000, input and IP libraries of normal and heat stress conditions were sequenced. Statistical analysis and quality control were performed on the original data generated by RNA sequencing (Table 1). We acquired almost 80 million reads per library. After screening and quality control of the original data generated by sequencing, valid reads were mapped to the reference genome (Table 2). Among the valid reads, about 80% were uniquely mapped to the reference genome ( Table 2). We next analyzed the distribution of m 6 A in the whole transcriptome for CK and T43 groups. The reads were mapped and distributed along CDS, 5 UTR, and 3 UTR, as depicted in Figure 1. The read frequency in 3 UTR of m 6 A-IP samples was significantly higher than that of input samples. In Arabidopsis, studies have found that a dominant m 6 A peak near the top codon or 3 UTR is observed in most nuclear mRNAs [11,19]. Here, in Pak-choi, our research found that m 6 A might be predominantly located in the 3 downstream region.    The m 6 A peaks (actually identified as m 6 A modification sites) were identified based on a comparison of read distribution between the input and IP samples using the MeTDiff package. In the T43 and CK groups, 15,919 and 15,436 m 6 A peaks were identified, respectively ( Figure 2A). The average Plants 2020, 9, 1080 4 of 12 length of peaks was 401.04 and 378.02 bp, respectively. The minimum length of peaks was 135 bp. In addition, 12,392 and 12,363 genes contained m 6 A peaks in T43 and CK groups, respectively ( Figure 2B). The Veen diagram identifies 11,252 common m 6 A peaks and 9729 common m 6 A peak-containing genes ( Figure 2). The GO enrichment analysis showed that the biological processes of RNA processing and methylation, salt stress response, and protein folding were significantly enriched.

m 6 A Topological Patterns in Pak-choi
Over 90% of the methylated transcripts showed one or two m 6 A sites, whether in the CK or T43 group (Figure 3), and about 80% contained one m 6 A site, which was much higher than previously reported in Arabidopsis (only 27%). A previous study of m 6 A topology in Arabidopsis chloroplasts and mitochondria found that only one m 6 A site in these two organelles covered more than 27% of the methylated transcripts [20].
Plants 2020, 9, x FOR PEER REVIEW 5 of 14 mitochondria found that only one m 6 A site in these two organelles covered more than 27% of the methylated transcripts [20].
In order to further understand the position of m 6 A in the transcripts, the subgene profiles of the m 6 A peaks in the transcriptomes of CK and T43 groups were studied. The m 6 A peaks were assigned based on their location on different transcript segments: 5′ UTR, 3′ UTR, and exon. The results showed that the m 6 A peaks were markedly located in the 3′ untranslated region (3′ UTR; Figure 4).  In order to further understand the position of m 6 A in the transcripts, the subgene profiles of the m 6 A peaks in the transcriptomes of CK and T43 groups were studied. The m 6 A peaks were assigned To determine whether the identified m 6 A peaks share common sequence elements that are characteristic of m 6 A RNA modification, we performed an unbiased search for consensus motifs enriched in regions surrounding the m 6 A peak identified in the CK group ( Figure 4C). The conserved To determine whether the identified m 6 A peaks share common sequence elements that are characteristic of m 6 A RNA modification, we performed an unbiased search for consensus motifs enriched in regions surrounding the m 6 A peak identified in the CK group ( Figure 4C). The conserved motifs, like AAACCV (V: U/A/G; p = 8.4 × 10 −32 ), were also enriched ( Figure 4C). In order to estimate the relationships between the m 6 A peak and expression levels, the genes were equally divided into 10 groups according to expression level, from low to high, and the m 6 A density of each group was calculated separately. The m 6 A peak density increased with the increase of expression levels. We also Plants 2020, 9, 1080 6 of 12 used this information to estimate that the genes with the highest expression levels contain about 0.75 m 6 A peaks per gene in the CK and T43 groups, respectively ( Figure 5).
Plants 2020, 9, x FOR PEER REVIEW 7 of 14 motifs, like AAACCV (V: U/A/G; p = 8.4 × 10 −32 ), were also enriched ( Figure 4C). In order to estimate the relationships between the m 6 A peak and expression levels, the genes were equally divided into 10 groups according to expression level, from low to high, and the m 6 A density of each group was calculated separately. The m 6 A peak density increased with the increase of expression levels. We also used this information to estimate that the genes with the highest expression levels contain about 0.75 m 6 A peaks per gene in the CK and T43 groups, respectively ( Figure 5). In order to study the influence of heat stress on m 6 A modification, differential m 6 A peaks were identified between the CK and T43 groups using MeTDiff. A total of 2603 differential m 6 A peaks were identified by comparing the CK and T43 groups (Table S1). Analysis of GO enrichment of the differential m 6 A peak-containing genes showed that the biological processes of virus-induced gene silencing, freezing response, and post-translational protein modification were enriched.

Differentially Expressed Genes Analysis
To investigate the potential relationship of m 6 A and gene expression, we used FPKM (fragments per kilobase of exon model per million mapped reads) to quantify the gene expression in T43 and CK groups using the input sequencing data. As a result, 31,832 genes were detected in at least one sample (Table S2). By comparing the T43 group and the CK group, we found that 7593 genes were upregulated, while 9379 genes were downregulated, after high-temperature stress ( Figure 6A). These results were further demonstrated by using gene expression density maps ( Figure 6B). The differentially expressed genes (DEGs) were categorized according to GO annotations ( Figure 7A). The biological processes, such as response to heat, response to cold, and response to salt stress, were enriched. Their molecular functions were primarily related to protein-binding activity. The cellular components involving the chloroplast envelope and chloroplast stroma were enriched. Through KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis, we found that most differential genes were enriched in the biosynthesis of amino acids and the translation processing pathway ( Figure 7B). Research findings indicated that changes in m 6 A might help regulate many genes expressed under stress [17]. In fact, Zhou et al. [21] found that m 6 A changes due to heat shock In order to study the influence of heat stress on m 6 A modification, differential m 6 A peaks were identified between the CK and T43 groups using MeTDiff. A total of 2603 differential m 6 A peaks were identified by comparing the CK and T43 groups (Table S1). Analysis of GO enrichment of the differential m 6 A peak-containing genes showed that the biological processes of virus-induced gene silencing, freezing response, and post-translational protein modification were enriched.

Differentially Expressed Genes Analysis
To investigate the potential relationship of m 6 A and gene expression, we used FPKM (fragments per kilobase of exon model per million mapped reads) to quantify the gene expression in T43 and CK groups using the input sequencing data. As a result, 31,832 genes were detected in at least one sample (Table S2). By comparing the T43 group and the CK group, we found that 7593 genes were upregulated, while 9379 genes were downregulated, after high-temperature stress ( Figure 6A). These results were further demonstrated by using gene expression density maps ( Figure 6B). The differentially expressed genes (DEGs) were categorized according to GO annotations ( Figure 7A). The biological processes, such as response to heat, response to cold, and response to salt stress, were enriched. Their molecular functions were primarily related to protein-binding activity. The cellular components involving the chloroplast envelope and chloroplast stroma were enriched. Through KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis, we found that most differential genes were enriched in the biosynthesis of amino acids and the translation processing pathway ( Figure 7B). Research findings indicated that changes in m 6 A might help regulate many genes expressed under stress [17]. In fact, Zhou et al. [21] found that m 6 A changes due to heat shock stress and activates hsp70 mRNA translation. We hypothesized that some of these differential genes might be involved in the amino acid metabolism synthesis pathway in-vivo. These findings call for further analysis of the relationship between m 6 A modification and regulation of transcriptional expression. stress and activates hsp70 mRNA translation. We hypothesized that some of these differential genes might be involved in the amino acid metabolism synthesis pathway in-vivo. These findings call for further analysis of the relationship between m 6 A modification and regulation of transcriptional expression.

Association Analysis between Differentially Expressed Genes and Differential m 6 A Peaks
To further associate m 6 A modification and gene expression, codifferential gene analysis was conducted. We combined the differentially expressed genes and differential m 6 A peaks. To visually represent the link between gene expression and m 6 A methylation change, we constructed a fourimage map of differential expression genes and differential m 6 A methylation genes ( Figure 8). As a result, 1641 genes were identified as significant codifferential genes with the following screening criteria: fold change ≥1.5 in m 6 A peak and expression levels ( Figure 8, Table S3). We found that 516 genes satisfied the significant difference in the upregulation of m 6 A and significantly downregulated their expression, and 714 differential genes satisfied the significant difference in the upregulation of m 6 A and significantly upregulated their expression. Similarly, 320 genes satisfied the significant difference in the downregulation of m 6 A and significantly downregulated their expression, and 91 genes satisfied the significant difference in the downregulation of m 6 A and significantly upregulated their expression (Figure 8). Previous studies have emphasized that m 6 A modification has a certain effect on plant development and stress resistance [22]. Further research revealed that m 6 A can inhibit site-specific transcriptome cleavage in plants, and this mechanism is necessary to properly regulate the salt-stress response transcriptome [23]. These results provide a solid foundation for further study of m 6 A modification.

Association Analysis between Differentially Expressed Genes and Differential m 6 A Peaks
To further associate m 6 A modification and gene expression, codifferential gene analysis was conducted. We combined the differentially expressed genes and differential m 6 A peaks. To visually represent the link between gene expression and m 6 A methylation change, we constructed a four-image map of differential expression genes and differential m 6 A methylation genes (Figure 8). As a result, 1641 genes were identified as significant codifferential genes with the following screening criteria: fold change ≥1.5 in m 6 A peak and expression levels ( Figure 8, Table S3). We found that 516 genes satisfied the significant difference in the upregulation of m 6 A and significantly downregulated their expression, and 714 differential genes satisfied the significant difference in the upregulation of m 6 A and significantly upregulated their expression. Similarly, 320 genes satisfied the significant difference in the downregulation of m 6 A and significantly downregulated their expression, and 91 genes satisfied the significant difference in the downregulation of m 6 A and significantly upregulated their expression (Figure 8). Previous studies have emphasized that m 6 A modification has a certain effect on plant development and stress resistance [22]. Further research revealed that m 6 A can inhibit site-specific transcriptome cleavage in plants, and this mechanism is necessary to properly regulate the salt-stress response transcriptome [23]. These results provide a solid foundation for further study of m 6 A modification.
Plants 2020, 9, x FOR PEER REVIEW 10 of 14 Figure 8. A four-image map to analyze the relationship between differential genes and differential peaks.

Plant Material and Tissue Collection
An excellent variety "Suzhouqing" from Pak-choi seeds were grown in pots with a soil matrix (soil:vermiculite = 3:1) and then placed in an artificial climate room (day: 16 h at 25 °C, the light intensity was set to 150 μmol m −2 s −1 ; night: 8 h at 18 °C). When the seedlings reach the five-leaf stage, some seedlings were moved to the growth room without light for heat stress treatment (43 °C for 8 h); other seedlings were used as a control group (25 °C without light). Control plants and treated plants were then sampled (leaves). The fresh weight of each sample was kept at 200 mg; fresh leaves were collected, and large veins were avoided as much as possible. All leaf samples were immediately Figure 8. A four-image map to analyze the relationship between differential genes and differential peaks.

Plant Material and Tissue Collection
An excellent variety "Suzhouqing" from Pak-choi seeds were grown in pots with a soil matrix (soil:vermiculite = 3:1) and then placed in an artificial climate room (day: 16 h at 25 • C, the light intensity was set to 150 µmol m −2 s −1 ; night: 8 h at 18 • C). When the seedlings reach the five-leaf stage, Plants 2020, 9, 1080 9 of 12 some seedlings were moved to the growth room without light for heat stress treatment (43 • C for 8 h); other seedlings were used as a control group (25 • C without light). Control plants and treated plants were then sampled (leaves). The fresh weight of each sample was kept at 200 mg; fresh leaves were collected, and large veins were avoided as much as possible. All leaf samples were immediately frozen in liquid nitrogen, then stored at −80 • C in a refrigerator. Three samples of each treatment were collected from all the above plants to be used in three biological replicates. At the same time, three trials were repeated.

Library Construction and RNA Sequencing
We used Trizol reagent (Invitrogen, Carlsbad, CA, USA) to extract total RNA from the samples obtained. The quality and quantity of the total RNA obtained were analyzed by a Bioanalyzer 2100 and an RNA 6000 Nano LabChip kit (Agilent, Palo Alto, CA, USA) with RIN value >7.0. About 200 µg of total RNA was isolated from poly (A) mRNA using magnetic beads (Invitrogen) with poly-T oligonucleotides. Then, the cleaved RNA fragment and m6A specific antibody were supplemented with BSA (0.5 µg µL −1 ). Then, we incubated those mixtures with protein A beads by elution using elution buffer. The eluted RNA was finally subjected to precipitation analysis with 75% ethanol. According to the chain-specific library prepared by the dUTP method [24], the eluted fragment containing m 6 A and the untreated input control fragment were converted into a cDNA library. The results showed that the average insertion size of the paired library was about 100 ± 50 bp. Finally, we performed 2 × 150 bp paired-end sequencing analysis on the Illumina Novaseq™ 6000 platform according to the method recommended by the supplier.

Data Analysis
According to the method described by our predecessors, sequence data analysis was performed on preprocessing sequencing reads and read alignments [25]. First, the original data of IP RNA-seq and input RNA-seq were analyzed by Trim Galore (version 0.3.7), adapters and low-quality (Q < 25) data were deleted, and data with readings less than 50 bp were deleted [26]. We also deleted all readings that map to multiple genomic regions, Keeping only those reads that were uniquely mapped to the reference sequence for further analysis of the m 6 A-modified peak.

Biological Information Analysis
To get clean data, we first used cutadapt [27] and local perl scripts to remove low-quality and contaminated sequences, followed by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/) software to clean data and conduct quality inspection and control. Finally, the default parameters of bowtie [28] were used to compare reads to the reference genome (data unpublished) to read the statistics. By using MeTDiff [29] and ChIPseeker [30], peak calling, differential m 6 A peak, and peak annotation analysis were performed, and, at the same time, by using MEME [31] and HOMER [32], motif analysis was performed on the enriched area. The gene quantification software used was StringTie [33], and the normalization mode was set to FPKM. The R package edgeR [34] was used to identify and analyze gene differences.

Expression and Function Analysis of Multilayer Genes
Cufflinks software (version 2.2.1) was used to measure gene expression levels per million unique mapping reads of genes per kilobase exon model input sequence reads, followed by the Cuffdiff program (version 2.2.1) to identify the number of differentially expressed genes between different samples. These genes (p ≤ 0.05 and normalization factor ≥1.5) are seen as differentially expressed genes [35]. GO analysis of the hypergeometric distribution test was based on the GO consortium database and R, using DAVID bioinformatics [36,37].

Conclusions
For the first time, we provide an m 6 A map on the transcriptome after heat-stress treatment in Pak-choi. More importantly, our map reveals the characteristics of m 6 A distribution in the Pak-choi transcriptome and its possible functional meaning. We also found a certain relationship between the degree of m 6 A methylation and the level of transcription, further indicating that m 6 A might play an important role in the regulation of gene expression. The construction of a comprehensive map in this study may provide a solid basis for determining the functional role of heat stress on RNA m 6 A modification in Pak-choi.
Author Contributions: G.L. and J.W. completed most of the experiments and wrote the paper; they contributed equally to this work. The manuscript was revised and approved by X.H. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors are grateful to Ling Liu and Song Shi of Ouyi Biological Technology Co. Ltd. for the provided sequencing platform and biological information analysis.

Conflicts of Interest:
The authors declare no conflict of interest.

Availability of Data and Materials:
The RNA-seq data were deposited in the SRA database under NCBI BioProject ID PRJNA646007. All datasets supporting the manuscript conclusions are included in the article and its attached files.