Abstract Phytoplasmas induce diseases in more than 1000 plant species and cause substantial ecological damage and economic losses, but the specific pathogenesis of phytoplasma has not yet been clarified. N ^6‐methyladenosine (m^6A) is the most common internal modification of the eukaryotic Messenger RNA (mRNA). As one of the species susceptible to phytoplasma infection, the pathogenesis and mechanism of Paulownia has been extensively studied by scholars, but the m^6A transcriptome map of Paulownia fortunei ( P. fortunei ) has not been reported. Therefore, this study aimed to explore the effect of phytoplasma infection on m^6A modification of P. fortunei and obtained the whole transcriptome m^6A map in P. fortunei by m^6A‐seq. The m^6A‐seq results of Paulownia witches' broom (PaWB) disease and healthy samples indicate that PaWB infection increased the degree of m^6A modification of P. fortunei . The correlation analysis between the RNA‐seq and m^6A‐seq data detected that a total of 315 differentially methylated genes were predicted to be significantly differentially expressed at the transcriptome level. Moreover, the functions of PaWB‐related genes were predicted by functional enrichment analysis, and two genes related to maintenance of the basic mechanism of stem cells in shoot apical meristem were discovered. One of the genes encodes the receptor protein kinase CLV2 (Paulownia_LG2G000076), and the other gene encodes the homeobox transcription factor STM (Paulownia_LG15G000976). In addition, genes F‐box (Paulownia_LG17G000760) and MSH5 (Paulownia_LG8G001160) had exon skipping and mutually exclusive exon types of alternative splicing in PaWB‐infected seedling treated with methyl methanesulfonate, and m^6A modification was found in m^6A‐seq results. Moreover, Reverse Transcription–Polymerase Chain Reaction (RT‐PCR) verified that the alternative splicing of these two genes was associated with m^6A modification. This comprehensive map provides a solid foundation for revealing the potential function of the mRNA m^6A modification in the process of PaWB. In future studies, we plan to verify genes directly related to PaWB and methylation‐related enzymes in Paulownia to elucidate the pathogenic mechanism of PaWB caused by phytoplasma invasion. Keywords: m^6A‐seq, N ^6‐methyladenosine, Paulownia fortunei, Paulownia witches' broom, phytoplasma 1. INTRODUCTION N ^6‐methyladenosine is the most common posttranscriptional modification of mRNA. It is regulated by methyltransferases, demethylases, and m^6A‐binding proteins (Deng et al., [38]2015; Dominissini et al., [39]2012; Luo et al., [40]2014; Meyer et al., [41]2012; Schwartz et al., [42]2013). The regulation of m^6A modification involves a complex dynamic system including processes related to embryonic development, apoptosis, spermatogenesis, and circadian rhythm (Jia et al., [43]2011; Liu et al., [44]2014; Ping et al., [45]2014; Yang et al., [46]2018; Zhao et al., [47]2017; Zheng et al., [48]2013; Zhong et al., [49]2008). If components of the m^6A modification system are defective, diseases or disorders, including tumors, neurological diseases, and embryonic developmental delays, can occur (Li et al., [50]2017; Lin et al., [51]2016; Liu et al., [52]2013; Ping et al., [53]2014; Zhang et al., [54]2016; Zhao et al., [55]2017; Zheng et al., [56]2013; Zhong et al., [57]2008). In addition, m^6A is involved in multiple RNA metabolic processes, including mRNA expression (Dominissini et al., [58]2012), translation efficiency (Wang et al., [59]2015), alternative splicing (Liu et al., [60]2015), and degradation (Wang et al., [61]2013). High‐throughput sequencing analysis of m^6A targeted antibodies in Arabidopsis showed that more than 70% m^6A peak was detected in Arabidopsis Can‐0 and Hen‐16, indicating that m^6A is a highly conservative modification of mRNA in plants (Luo et al., [62]2014). Subsequently, differential m^6A methylation patterns among three organs of Arabidopsis were analyzed using transcriptome wide high‐throughput deep m^6A‐seq. The results showed that over 80% of m^6A modified transcripts were identified in leaves, flowers, and roots of Arabidopsis (Wan et al., [63]2015). The study of m^6A patterns in the chloroplast and mitochondrial transcriptome of Arabidopsis shows that the m^6A motif that undergoes RNA methylation is highly conserved between the nucleus and the organelle transcriptome (Wang et al., [64]2017). The m^6A map of sea buckthorn (Hippophae rhamnoides Linn.) transcriptome identified 13,287 different m^6A peaks between leaf under drought and in control treatment. In addition, it is reported that after the m^6A modification disorder, the leaves and roots of plants have serious dysplasia, and the time of reproductive development changes (Zhang et al., [65]2021). These findings indicate that m^6A has a regulatory role in plant gene expression. In recent years, with the deepening of research on m^6A writers, erasers, and readers in plants, it has been found that m^6A methyltransferase in Arabidopsis functions by forming complexes, and the accessory subunits FIP37 and VIR act as key subunits to maintain the function of m^6A methyltransferase complexes (Shen, [66]2023). At the same time, many m^6A methylation‐related enzymes that play a crucial role in plant growth and development have been identified. FIP37, a homolog of WTAP in Arabidopsis, is a core component of the methyltransferase complex that mediates the methylation of transcription factors STM and WUS and affects the development of aboveground organs (Shen et al., [67]2016). ALKBH10B is an RNA m^6A demethylase that regulates inflorescence transformation in Arabidopsis (Duan et al., [68]2017). The m^6A reader ECT2 recognizes the site of m^6A modification and regulates the morphogenesis of Arabidopsis trichomes (Wei et al., [69]2018). RNA‐binding protein FLK is a novel mRNA m^6A reader protein that directly binds to the m^6A site in the 3′‐untranslated region of FLC transcript, thus repressing FLC levels by reducing its stability and splicing, thereby regulated floral transition in Arabidopsis (Amara et al., [70]2023). M^6A‐related proteins were also identified in tobacco. In tobacco infected with Tobacco mosaic virus (TMV), the gene expression level of m^6A demethylase ALKBH5 increased, and the results showed that TMV reduced the m6A level (Li et al., [71]2018). Until now, m^6A modification has rarely been investigated in woody plants, so it is necessary to study m^6A in Paulownia. Paulownia is an important fast‐growing timber and ornamental tree species indigenous to China, but now planted worldwide. Paulownia witches' broom (PaWB), a plant disease caused by phytoplasma infection, seriously affects the growth and development of Paulownia and even leads death, resulting in serious economic loss and ecological damage. Phytoplasmas are parasitic prokaryotes without cell wall and are transmitted by insect vectors such as psyllids, planthoppers (Kosovac et al., [72]2018), and woodlice. Phytoplasma is the causative agent of more than 1000 plant diseases, including Jujube witches' broom, PaWB, and Aster yellow witches' broom, which are extremely harmful to plants (Geng et al., [73]2015). Advances in biotechnology have made it possible to study posttranscriptional mRNA modifications, especially m^6A modifications, which have become the focus of many studies. Given previous findings on m^6A modification, most of the neurological diseases, embryonic developmental delays, and tumorigenesis caused by m^6A modification disorders involve stem cell stability. Arbuscule symptoms in PaWB are also caused by stem cell disorder in the stem apical meristem. Therefore, it is necessary to investigate whether the occurrence of PaWB is affected by m^6A modification. Previous studies have shown that phytoplasma‐infected Paulownia seedlings restored to normal morphology after treatment with 60 mg·L^−1 methyl methanesulfonate (MMS), and no phytoplasma could be detected (Wang et al., [74]2018). To investigate the changes in m^6A modification during recovery of phytoplasma‐infected Paulownia seedlings after MMS reagent treatment and to screen PaWB‐related genes, the terminal buds of PaWB‐infected seedlings (PFI) and PaWB‐infected seedlings treated with 60 mg·L^−1 MMS (PFIM60) were collected, and the m^6A modification map in Paulownia was obtained by m^6A‐Seq, and the differentially expressed genes with m^6A modification were further screened, by combining with transcriptome analysis. In view of the function of m^6A in regulating alternative splicing, the genes modified by m^6A with alternative splicing were analyzed in order to understand the gene expression changes of Paulownia caused by phytoplasma infection. Understanding the changes in mRNA m^6A modification in Paulownia after phytoplasma infection will lay a solid scientific foundation for the discovery of PaWB‐related pathogenesis and epigenetic regulation mechanisms. 2. MATERIALS AND METHODS 2.1. Plant materials The experimental plant materials used in this study were obtained from the Forest Biotechnology Laboratory of Paulownia Research Institute, Henan Agricultural University, Zhengzhou, China, and were harvested with permission. The plant materials were grown on 1/2 MURASHIGE & SKOOG (MS) medium in 100‐mL Erlenmeyer flasks for 30 days. Then, 1.5‐cm terminal buds from the PFI plants were transferred into 1/2 MS medium containing 60 mg·L^−1 MMS (PFIM60 samples) (Figure [75]1b) or 0 mg·L^−1 MMS (PFI samples) (Figure [76]1a) (Wang et al., [77]2018). Two seedlings were cultured in each bottle, and 90 bottles were cultivated for each sample. After 30 days, the witches' broom symptoms of the PFIM60 samples returned to normal significantly compared with the PFI samples grown for the same 30 days. The obtained samples were frozen in liquid nitrogen and stored at −80°C until used. FIGURE 1. FIGURE 1 [78]Open in a new tab Morphological changes in phytoplasma‐infected Paulownia fortunei seedlings. (a) Phytoplasma‐infected P. fortunei seedling. The diseased seedlings have obvious symptoms of witches' broom. (b) Seedlings treated with 60 mg·L^−1 methyl methanesulfonate. 2.2. Antibody enrichment and sequencing library construction Total RNA was extracted using a RNAsimple Total RNA Kit (TIANGEN, Cat.#DP419). More than 200 μg total RNA was used to isolate poly(A) mRNA with poly(T) oligo‐attached magnetic beads (Invitrogen). After purification, the mRNA fractions were fragmented into approximately 100‐nt long oligonucleotides using divalent cations under elevated temperature. The cleaved RNA fragments were incubated for 2 h at 4°C with m^6A‐specific antibody (No. 202003, Synaptic Systems, Germany) in IP buffer (50‐mM Tris–HCl, 750‐mM NaCl, and .5% Igepal CA‐630) supplemented with bovine serum albumin (.5 μg·μL^−1). The mixture was incubated with protein‐A beads and eluted with elution buffer (1 × IP buffer and 6.7‐mM m^6A). The eluted RNA was precipitated with 75% ethanol. The eluted m^6A‐containing fragments (IP) and untreated input control fragments were converted to form the final cDNA library according to the strand‐specific library preparation dUTP method (Levin et al., [79]2010). We constructed four m^6A‐seq libraries (PFI_1_IP, PFI_2_IP, PFI‐M60_1_IP, and PFI‐M60_2_IP) and four RNA‐seq libraries (PFI_1_IP, PFI_2_IP, PFI‐M60_1_IP, and PFI‐M60_2_IP). The average insert size for the paired‐end libraries was approximately 100 ± 50 bp. Paired‐end 2 × 150 bp sequencing was performed on an Illumina Novaseq™ 6000 (Biotech ltd, Hangzhou, China) platform following the manufacturer's recommended protocol. 2.3. Alignment of reads to the reference genome and visualization of m^6A peaks Cutadapt (Martin, [80]2011) and perl in‐house scripts were used to remove reads that contained adaptor contamination, low quality bases, and/or undetermined bases. Sequence quality was verified using FastQC (Version 0.11.1) ([81]http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Bowtie (Langmead & Salzberg, [82]2012) (with default parameters) was used to map reads to the Paulownia fortunei reference genome. Mapped reads of the IP and input libraries were input into the exomePeak package in R (Meng et al., [83]2014) to identify m^6A peaks (in BED or BAM format) for visualization using the UCSC Genome Browser or Integrative Genomics Viewer software ([84]http://www.igv.org/). 2.4. Transcription level analysis StringTie (Pertea et al., [85]2015) software was used to determine the mRNA expression levels in the input libraries by calculating Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) (FPKM = [total exon fragments / mapped reads (millions) × exon length (kb)]). The differentially expressed mRNAs were selected using the edgeR package in R with a threshold of log2 (fold change) >1 or <−1 and p‐value < 0.05 (Robinson et al., [86]2010). 2.5. Discernment of m^6A topological patterns MEME (Bailey et al., [87]2009) and HOMER (Heinz et al., [88]2010) were used to find de novo and known motifs in the differential modified genes. Localization of the motifs with respect to peak summit was performed using in‐house perl scripts. Called peaks were annotated by intersection with gene architecture using ChIPseeker (Yu et al., [89]2015). 2.6. M^6A MeRIP‐reverse transcription PCR The m^6A methylated RNA immunoprecipitation (MeRIP) assay was adapted from a previously reported protocol (Dominissini et al., [90]2013). Briefly, approximately 300‐mg total RNA was extracted and isolated using an mRNA screening kit. Each biological replicate for m^6A‐seq started with 400 mg of total RNA, which yielded approximately 10 mg of double poly(A) selected mRNA. Then, the mRNA was treated with mRNA fragmentation reagent, and fragmentation of the 300‐mL poly(A) RNA solution was performed at 94°C for exactly 5 min using a thermocycler. The fragmentation reaction was stopped by adding 50 mL of stop buffer to a final volume of 350 mL and immediately put on ice. Approximately 100‐ng mRNA was taken out as the input sample, and 150 mL of pre‐equilibrated m^6A‐Dynabeads was added to the remaining fragmented RNA to a final volume of 500 mL. The fragmented RNA was allowed to bind to the m^6A‐Dynabeads at 4°C while rotating (tail‐over‐head) at seven rotations per minute for 2 h. Then, the tubes containing the samples were placed on a magnet, and the bead complexes were allowed to cluster until the solution became clear. The 500‐mL liquid phase or supernatant was discarded because it contained m^6A negative fragments that were not captured by the anti‐m^6A antibody. The m^6A‐Dynabeads–RNA complexes were resuspended in 500‐mL m^6A binding buffer and incubated for 3 min at room temperature, and the clear supernatant was removed after placing the beads on a magnet. This step was repeated with 500‐mL low salt buffer and then with 500‐mL high salt buffer. (The 3‐min incubation time should not be exceeded to prevent release of the RNA from the beads.) This step was then repeated twice with 500‐mL Tris–EDTA (TE) buffer. The obtained immunoprecipitated product was prepared by resuspending in 100 mL of proteinase K buffer, and the obtained input sample was thawed. Then, 107 mL of RIP wash buffer, 15 mL of 10% sodium dodecyl sulfate (SDS), and 18 mL of proteinase K were added to a total volume of 150 mL, and the mixture was incubated at 55°C for 30 min with shaking to digest the antibody from the magnetic beads. After incubation, the supernatant was removed, and 250 mL of RIP wash buffer was added to the supernatant. Next, 400 mL of phenol/chloroform was added to each tube, which was then vortexed for 15 s, and .1 times the volume was added. RNA precipitation was promoted by adding 3‐M NaAc and 2‐mL GlycoBlue at −20°C for 1 h. The tubes were centrifuged at 14,000 g for 20 min at 4°C, 75% ethanol precipitation was performed twice, and the pellet was blown dry and then dissolved in 10‐mL RNA‐free water and tested. Primer sequences are shown in Table [91]S1. 3. RESULTS 3.1. The m^6A modification landscape of P. fortunei PaWB‐infected P. fortunei (PFI) seedlings and PFI seedlings treated with 60 mg·L^−1 MMS (PFIM60) were used for m^6A‐seq, RNA‐seq (total fragmented RNA as the control for m^6A‐seq), and transcriptome analysis, with two replicates for each. Heatmap correlation coefficient analysis between the biological replicates confirmed their high repeatability (Figure [92]2). We obtained 44–56 and 73–83 million reads for the RNA‐seq and m^6A‐seq libraries, respectively (Table [93]S2), of which 23–28 and 47–50 million were uniquely aligned to the P. fortunei reference genome (Table [94]S3). Almost 94% of reads aligned to the reference genome were located in exons; the others were located in introns or intergenic regions (Figure [95]S1 and Table [96]S4). FIGURE 2. FIGURE 2 [97]Open in a new tab Heatmap correlation coefficient analysis between biological replicates. A m^6A modification map of P. fortunei was constructed, showing the distribution of m^6A on chromosomes, as well as the distribution of m^6A peak positions in PFI and PFIM60 on chromosomes (Figure [98]3). The number of m^6A‐modified genes on chromosomes ranged from 472 to 1014; chromosomes 9 (1014) and 2 (1002) have more modified genes than the other chromosomes. The highest number of m^6A peaks (1506) was on chromosome 2, and the lowest number (702) was on chromosome 14. The number of methylated genes and peaks on chromosomes was similar in PFI and PFIM60. FIGURE 3. FIGURE 3 [99]Open in a new tab M^6A modification map of Paulownia fortunei. The outermost circle is the distribution of chromosome, the second circle is the distribution of PFI m^6A peak positions on chromosomes, and the third circle is the distribution of PFIM60 m^6A peak positions on chromosomes. In PFIM60 samples, 20,201 peaks of 13,505 genes were identified, and in PFI samples, 20,568 peaks of 13,838 genes were identified. The degree of m^6A modification of PFI was slightly higher than that of PFIM60, which indicated that PaWB phytoplasma infection caused the change of m^6A modification of Paulownia (Dataset [100]S1). Notably, there was greater enrichment of peaks at the transcription start and stop sites of genes in both samples (Figure [101]4a), suggesting that m^6A modification may play a key role in regulating mRNA degradation and microRNA processing in Paulownia. The differences in m^6A peaks observed between PFI and PFIM60 were compared and analyzed. The results showed that there were 306 peaks of 319 genes with reduced m^6A modification in PFIM60 compared with their m^6A modification in PFI; 36.59% of them were located in 3′ Untranslated Regions (UTRs), 17.41% were in 5′ UTRs, 21.58% were in first exons, and 24.15% were in other exons. This distribution was consistent with m^6A enrichment of all the genes in the P. fortunei genome (Figure [102]4b). FIGURE 4. FIGURE 4 [103]Open in a new tab Peak enrichment area. (a) Peak enrichment of transcription start sites for four samples, and the peak distribution of the bondable region near the Transcription Start Sites (TSS) is shown in the form of a heatmap. (b) Peak enrichment of genetic elements. 3.2. Differential modified genes between PFI and PFIM60 To determine the sequence specificity of the hypermethylated m^6A peaks, we performed motif prediction based on differential peak analysis using HOMER and compared the predicted motif with microRNA sequences in miRBase (1e^−20 < p < 1e^−10). A total of 10 highly enriched motif sequences were identified. Among them, the motif RRACH (R = A/G, H = A/C/U) found in Paulownia has not been confirmed in mammals and yeast. The most significantly enriched sequence was UUGUUUUGUACU (Figure [104]5), which is similar to the UGUAYY (Y = C/U) motif that is specific to m^6A modifications in tomato, rice, corn, and other plants, indicating a high degree of confidence in the predictions. These motif sequences are recognized and bound by m^6A‐related enzymes, which in turn affect gene expression. FIGURE 5. FIGURE 5 [105]Open in a new tab Motifs with higher proportions: UUGUUUUGUACU. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed on the peak‐related genes with weakened m^6A modification detected in the PFIM60 samples. The results showed that auxin biosynthesis, hormone signal transduction, and Mitogen–Activated Protein Kinase (MAPK) signal pathways were some of the most abundant (Figure [106]6 and Dataset [107]S2). The auxin biosynthesis pathway involves three genes: Paulownia_LG2G000124, Paulownia_LG2G000125, and Paulownia_LG2G000128. Seven genes are involved in the plant hormone (auxin, abscisic acid, and ethylene) signal transduction pathway: Paulownia_LG11G000159, Paulownia_LG12G000675, Paulownia_LG17G000489, Paulownia_LG17G000559, Paulownia_LG2G001457, Paulownia_LG4G000193, and Paulownia_LG4G000193000451. The MAPK signaling pathway involves six genes: Paulownia_LG10G000799, Paulownia_LG11G000159, Paulownia_LG17G000559, Paulownia_LG2G001457, Paulownia_LG4G000451, and Paulownia_LG9G000737. FIGURE 6. FIGURE 6 [108]Open in a new tab Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of genes with reduced m^6A modification after 60 mg·L^−1 methyl methanesulfonate treatment. 3.3. Extent of m^6A modification versus gene transcript levels in PFI and PFIM60 First, the transcription levels of the mRNAs in the PFI and PFIM60 libraries were analyzed by calculating FPKM (FPKM = [total exon fragments / mapped reads (millions) × exon length (kb)]). Differentially expressed mRNAs (DEGs) were identified using the edgeR package in R (log2 [fold change] >1 or <−1 and p‐value < 0.05). A total of 5899 DEGs were detected between the PFI and PFIM60 libraries; 3130 were upregulated and 2769 were downregulated in PFIM60. Gene ontology enrichment analysis showed that DEGs were mainly involved in membrane, cell wall, DNA‐templated transcriptional regulation, secondary metabolite biosynthesis, DNA‐binding transcription factor activity, and oxidoreductase activity (Figure [109]7a). KEGG pathway enrichment analysis showed that DEGs were mainly involved in plant circadian rhythm, phenylpropanoid biosynthesis, flavonoid biosynthesis, and glycerolipid metabolism pathways (Figure [110]7b). FIGURE 7. FIGURE 7 [111]Open in a new tab Functional annotation and pathway enrichment of differentially expressed genes. (a) Gene ontology (GO) functional annotation of differentially expressed genes. (b) Statistics of Kyoto Encyclopedia of Genes and Genomes pathway enrichment. To examine the extent of m^6A modification versus gene transcript levels in the PFI and PFIM60 samples, we correlated m^6A‐seq data with transcriptome sequences. In the transcriptome, differences in transcript abundance were categorized as upregulation and downregulation, while in m^6A‐seq data, changes in peak abundance indicated m^6A modification levels of upregulated genes. A total of 315 differentially methylated genes were predicted to be significantly differentially expressed at the transcriptome level. Four types of relationships between m^6A modification levels and gene expression were identified: (1) both m^6A modification levels and transcription levels upregulated (133 genes); (2) m^6A modification levels upregulated and transcription levels downregulated (58 genes); (3) m^6A modification levels downregulated and transcription levels upregulated (44 genes); and (4) m^6A modification levels and transcription levels both downregulated (80 genes) (Figure [112]8 and Dataset [113]S3). Gene ontology enrichment analysis of these differentially methylated genes revealed that plasma membrane and chloroplast were enriched under the cellular component category, DNA‐based template for transcriptional regulation and response to light stimulus were enriched under the biological process category, and DNA‐binding transcription factor activity and sequence‐specific DNA‐binding functions were enriched under the molecular function category (Figure [114]9a). KEGG pathway analysis showed that these differentially methylated genes were mainly involved in signal pathways such as ABC transport, diterpenoid biosynthesis, plant hormone signal transduction, and phenylpropanoid biosynthesis (Figure [115]9b). FIGURE 8. FIGURE 8 [116]Open in a new tab Four relationships were detected between m^6A modification levels and gene expression. FIGURE 9. FIGURE 9 [117]Open in a new tab Functional annotations and pathways of differential peaks and differentially expressed genes. (a) Statistics of gene ontology (GO) enrichment. (b) Statistics of pathway enrichment. 3.4. M^6A modification affects alternative splicing Alternative splicing is involved in many physiological processes, including responses to abiotic and biotic stresses (Barbazuk et al., [118]2008). M^6A modification can affect alternative splicing by promoting retention of exons at the transcriptional level; therefore, we jointly analyzed m^6A modification and alternative splicing. In PFIM60, 282 genes had significant exon skipping‐type alternative splicing differences and 84 genes had significant mutually exclusive exon‐type alternative splicing differences compared with the corresponding genes in PFI. When m^6A modification was higher in PFIM60 than in PFI, alternative splicing was associated with m^6A modification in two genes, namely, F‐box (Paulownia_LG17G000760) and MSH5 (Paulownia_LG8G001160). To validate the alternative splicing of these two genes in PFI and PFIM60, RT‐PCR analysis was performed. The results indicated that for PFIM60, when the degree of m^6A modification increased, the bands of two genes were weakened compared with those for PFI, suggesting that the degree of m^6A modification affected alternative splicing events (Figures [119]10 and [120]S2). Primer sequences are shown in Table [121]S5. FIGURE 10. FIGURE 10 [122]Open in a new tab RT‐PCR validation of F‐box and MSH5 in PFI and PFIM60, respectively. 18s rRNA is the internal reference gene. 3.5. Identification of potential candidate genes associated with PaWB In the process of invading a host, intrinsic pathogen‐associated molecular patterns known as PAMPs or damage‐related molecular patterns in phytoplasmas are recognized by pattern recognition receptors in the host defense system, which induces PAMP‐triggered immunity and activates the downstream defense pathway. To counteract PAMP‐triggered immunity, phytoplasmas release effectors as signaling molecules to regulate signal transduction and disrupt host defense responses. This in turn triggers a series of physiological and biochemical reactions, which results in metabolic disruption in the host plant and development of disease symptoms. Some of the identified DEGs were related to plant–pathogen interactions, plant hormone signal transduction, and plant tillering. Among them, four significantly DEGs were involved in the plant hormone signal transduction pathways: Paulownia14027, encoding SHORT‐ROOT (LSH4); Paulownia26350, encoding LIGHT‐DEPENDENT SHORT HYPOCOTYLS 4 (LSH4); Paulownia_LG17G000277, encoding histidine‐containing phosphate transfer protein 1‐like; and Paulownia_LG2G000140, which encodes regulatory protein NPR5‐like. One DEG was involved in the plant–pathogen interaction pathway: Paulownia000609, which encodes disease resistance protein RPM1‐like. Two DEGs were involved in the carotenoid biosynthesis pathway: Paulownia_LG11G000307, which encodes phytoene synthase 2, and Paulownia_LG18G000559, which encodes carotenoid cleavage dioxygenase 4. The most prominent symptom of PaWB is axillary bud clustering, so we also focused on genes that affect branching and tillering in plants, namely, CLAVATA1 (CLV1), CLV2, CLV3 (Nikolaev et al., [123]2007), WUSCHEL (WUS) (Lenhard et al., [124]2002), Arabidopsis TOPLESS‐like transcriptional corepressor (ASP1), SHOOT MERISTEMLESS (STM), ARABIDOPSIS RESPONSE REGULATOR 5 (ARR5), ARR6, ARR7, and ARR15 (Leibfried et al., [125]2005). In Arabidopsis, these genes encode proteins involved in the WUSCHEL–CLAVATA negative feedback loop, which is the essential mechanism for maintaining stem cells in shoot apical meristem (Shen et al., [126]2016). Paulownia homologs of CLV2 (Paulownia_LG2G000076) and STM (Paulownia_LG15G000976) were discovered by mapping Arabidopsis genes to the Paulownia reference genome using the local BLAST program. These two genes showed significant differences in m^6A modification and expression levels after phytoplasma infection. We performed m^6A MeRIP quantitative reverse transcription PCR verification on these two candidate genes. The results were consistent with the sequencing results, which verified the accuracy of the sequencing results (Figure [127]11). FIGURE 11. FIGURE 11 [128]Open in a new tab Validation of m^6A methylated RNA immunoprecipitation quantitative reverse transcription PCR. In PFI samples infected with phytoplasma, the m^6A modification degree of STM and CLV2 was higher than that of healthy PFIM60 samples. 3.6. STM and CLV2 maintain the dynamic balance of shoot apical meristem regulated by m^6A modification Conserved domains were identified in the proteins encoded by the Paulownia homologs of STM and CLV2 (Table [129]1). Similar to STM, Paulownia_LG15G000976 encodes a homeobox transcription factor with four conserved domains: KNOX2, KNOX1, Homeobox_KN, and ELK. Similar to CLV2, Paulownia_LG2G000076 encodes a protein kinase that contains the domain [130]PLN00113 superfamily. Previous studies have shown that STM and WUS can activate the transcription and expression of CLV2 by binding with CLV2 promoter, while CLV3, CLV1, and CLV2 jointly maintain the number of undifferentiated stem cells in stem tip meristem and regulate the production of normal stem tip. Transcriptome data analysis showed that STM was more highly expressed in PFI than in PFIM60, while m^6A‐seq data analysis showed that CLV2 had higher methylation levels in PFI than in PFIM60. TABLE 1. Conserved domain analysis of two proteins related to Paulownia witches' broom. Gene Gene ID Significant m^6A Gene Incomplete Incomplete Incomplete Incomplete STM AT1G62360 — — — KNOX2 KNOX1 Homeobox_KN ELK Paulownia_LG15G000976 Yes Down Down KNOX2 KNOX1 Homeobox_KN ELK Paulownia_LG14G000617 Yes Down Down KNOX2 KNOX1 Homeobox_KN ELK Paulownia_LG7G001667 No Up Down KNOX2 KNOX1 Homeobox_KN ELK CLV2 AT1G65380 — — — [131]PLN00113 superfamily — — — Paulownia_LG2G000076 Yes Down Down [132]PLN00113 superfamily — — — [133]Open in a new tab 4. DISCUSSION 4.1. M^6A modification of P. fortunei before and after PaWB The m^6A modification, evident in various tissues, reportedly regulates the gene expression and executes corresponding biological functions by regulating the RNA metabolism, splicing, degradation, and translation (Fu et al., [134]2014). In recent years, RNA modification has become an important research field in the study of posttranscriptional gene expression regulation. Many studies related to m^6A have been conducted in the model plants such as Arabidopsis thaliana and tobacco, but there are less reports on the map of m^6A in woody plants. To our knowledge, this is the first comprehensive and high‐throughput study on RNA methylation in the Paulownia. Our data indicate that overall methylation levels were slightly higher in PFI samples compared with levels in PFIM60 samples. However, the differently methylated regions were similar in both samples, with nearly 40% of the modified regions located in the 3′ UTRs. Through further analysis, we found that m^6A modification may play an important role in the infection of PaWB phytoplasma. 4.2. M^6A modification regulates the expression of PaWB‐related genes Phytoplasma infection of P. fortunei releases effectors that disturb the plant's immune response. The effectors can interact directly with key plant immunity‐related genes or disrupt immune response by disrupting immune signal transduction. The association analysis of m^6A‐seq data and transcriptome sequences identified genes associated with plant–pathogen interaction, plant hormone signal transduction, and plant tillering. LSH4 (Paulownia26350) was predicted to be involved in plant–pathogen interactions and signal transduction pathways. Constitutive expression of LSH4 at the shoot apex inhibited leaf growth at the vegetative stage as well as the formation of extra shoots or shoot organs in flowers at the reproductive stage. Methylation modification and transcription levels of LSH4 were significantly upregulated in PFI compared with levels in PFIM60. We speculated that Paulownia phytoplasma infection increased the level of m^6A modification in the 3′ UTR region of LSH4, which enhanced the stability of LSH4 mRNA and increased its expression level, thereby promoting the occurrence of PaWB. The symptoms of axillary bud clustering caused by PaWB are closely related to the development of shoot apical meristem. In Arabidopsis, three master regulators, STM, WUS, and CLV3, are critical for stem cell specification in shoot apical meristem. STM interacts with WUS in vivo and recruits CLV3 to express shoot apical meristem (Brand et al., [135]2002). In Paulownia, WUS expression was relatively low and was similar in PFI and PFIM60, whereas STM expression in PFI was higher than in PFIM60. We speculated that the increased methylation level in the 3′ UTR region of STM in PFI maintained its stable expression. In Arabidopsis, the m^6A methylase FIP37 affected the expression of WUS and STM and regulated the development of shoot apical meristem (Shen et al., [136]2016). These results suggested that m^6A methylase may also be involved in regulating STM expression in Paulownia. However, this possibility needs to be further confirmed by molecular biology experiments. 4.3. The relationship between m^6A modification and alternative splicing of PaWB‐related genes Alternative splicing is a key regulatory mechanism that participates in many physiological processes, controls plant development, and increases the complexity of proteome and transcriptome. In addition, alternative splicing plays a crucial role in plant defense response and photosynthesis (Wang & Brendel, [137]2006). Compared with PFI, a total of 282 SE‐type alternatively spliced genes and 84 MXE‐type alternatively spliced genes were found in PFIM60. To determine the relationship between alternative splicing and PaWB, the genes associated with alternative splicing events in PaWB‐infected Paulownia were investigated. Combining m^6A modification and alternative splicing analysis, it was found that when m^6A modification increased, alternative splicing was associated with two m^6A modified genes, namely, F‐box (Paulownia_LG17G000760) and MSH5 (Paulownia_LG8G001160). The gene encoding F‐box is involved in alternative splicing events. The F‐box protein is a subunit of the Skp1–Cul1–F–box (SCF) complex and plays an important role in plant hormone signal transduction and regulation of plant development and growth (Gonzalez et al., [138]2017). Jasmonic Acid (JA) can regulate the growth of most plants, including their ability to adapt to adversity and to resist pathogen invasion. The F‐box protein COI1 is the receptor for JA, and the JA‐dependent response as the core is mediated by the Skp1/Cullin/F‐box E3 ubiquitin ligase complex containing COI1 (Gonzalez et al., [139]2017). The SCF complex can degrade JA, resulting in the release of MYC2 to express genes responsive to JA (Liu et al., [140]2017). In short, F‐box proteins may be involved in regulating Paulownia growth and improving its ability to adapt to stress. Due to the lack of reliable data, this speculation still needs further research to verify. MSH5 is a homolog of DNA mismatch recognition protein (MutS). DNA mismatch repair is a highly conserved biological pathway that can improve the fidelity of DNA replication and recombination. When the MutS protein recognizes small loops of mismatched and unpaired nucleotides, it initiates mismatch repair (Jiricny, [141]2013; Spampinato, [142]2017). Previous studies have shown that Arabidopsis and other plants encode MutS protein homologs (MSH) that are conserved in other eukaryotes. Plants lacking MSH7 show significantly faster and longer germination rates during the juvenile vegetative period. In primary roots, the number of stem leaves, axillary, and lateral inflorescences is higher than that of wild type (Chirinos‐Arias & Spampinato, [143]2020; Gonzalez & Spampinato, [144]2020). These findings suggest that the growth of mutant plants seems to be caused by damage to cell cycle checkpoints, which allow cells to divide without proper DNA repair, and its role in the process of phytoplasma infection of Paulownia remains to be further explored. 5. CONCLUSIONS In this research, m^6A‐seq was applied to determine methylation levels in phytoplasma‐infected Paulownia. Results indicated that m^6A modifications differed from those previously reported in model plants. We speculated that the Paulownia genome had a novel m^6A modification motif. DEGs associated with PaWB were detected by transcriptomic analysis. The effect of m^6A modification on alternative splicing was also analyzed. In future studies, we plan to verify genes directly related to PaWB and methylation‐related enzymes in Paulownia to elucidate the pathogenicity mechanism of PaWB caused by phytoplasma invasion. AUTHOR CONTRIBUTIONS GF conceived and designed these experiments. HY, YC, and YF provided suggestions on experimental design and analysis. PX and SH analyzed the data and wrote the paper. ZX, XL, HY, YC, and GF revised the manuscript. All authors read and approved the manuscript. CONFLICT OF INTEREST STATEMENT The authors declare that they have no competing interests. PEER REVIEW The peer review history for this article is available in the [145]Supporting Information for this article. Supporting information Figure S1. Locating of reads after aligned to the reference genome. Almost 94% of the reads that were aligned to the reference genome were locating in exons; the others were located in introns or intergenic regions. Figure S2. Gel imaging verified by RT‐PCR of F‐box and MSH5 in PFI and PFIM60, respectively. F‐box‐1, F‐box‐2, F‐box‐3 and MSH5–1, MSH5–2, MSH5–3 are the different degrees of exposure for RT‐PCR gel imaging of F‐box and MSH5, respectively. In each electrophoresis image, except for the Marker, the first three lanes are PFI samples, and the last three lanes are PFI‐M60 samples. Table S1. Primer sequence used for MeRIP‐qRT‐PCR in this study. Table S2. Sequencing sequence statistics and quality control. Table S3. Reference genome comparison to Reads statistics. Table S4. Regional distribution. Table S5. Primer sequence used for RT‐PCR in this study. [146]Click here for additional data file.^ (761.8KB, pdf) Dataset S1. M6A peaks identified in PFI and PFIM60 samples (PFI_peak). [147]Click here for additional data file.^ (2.6MB, xlsx) Dataset S2. Down or up peaks of m6A modification in PFI and PFIM60 (up). [148]Click here for additional data file.^ (121.5KB, xlsx) Dataset S3. Genes with up‐ or down‐regulated expression in PFI and PFIM60 (down). [149]Click here for additional data file.^ (2.8MB, xlsx) Dataset S3. Genes with up‐ or down‐regulated expression in PFI and PFIM60 (down). [150]Click here for additional data file.^ (64.1KB, docx) ACKNOWLEDGMENTS