N6-methyladenosine modification changes during the recovery processes for Paulownia witches' broom disease under the methyl methanesulfonate treatment

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⁶-methyladenosine (m⁶A) 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⁶A transcriptome map of Paulownia fortunei (P. fortunei) has not been reported. Therefore, this study aimed to explore the effect of phytoplasma infection on m⁶A modification of P. fortunei and obtained the whole transcriptome m⁶A map in P. fortunei by m⁶A-seq. The m⁶A-seq results of Paulownia witches' broom (PaWB) disease and healthy samples indicate that PaWB infection increased the degree of m⁶A modification of P. fortunei. The correlation analysis between the RNA-seq and m⁶A-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⁶A modification was found in m⁶A-seq results. Moreover, Reverse Transcription–Polymerase Chain Reaction (RT-PCR) verified that the alternative splicing of these two genes was associated with m⁶A modification. This comprehensive map provides a solid foundation for revealing the potential function of the mRNA m⁶A 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.

INTRODUCTION

N⁶-methyladenosine is the most common posttranscriptional modification of mRNA. It is regulated by methyltransferases, demethylases, and m⁶A-binding proteins (Deng et al., 2015; Dominissini et al., 2012; Luo et al., 2014; Meyer et al., 2012; Schwartz et al., 2013). The regulation of m⁶A modification involves a complex dynamic system including processes related to embryonic development, apoptosis, spermatogenesis, and circadian rhythm (Jia et al., 2011; Liu et al., 2014; Ping et al., 2014; Yang et al., 2018; Zhao et al., 2017; Zheng et al., 2013; Zhong et al., 2008). If components of the m⁶A modification system are defective, diseases or disorders, including tumors, neurological diseases, and embryonic developmental delays, can occur (Li et al., 2017; Lin et al., 2016; Liu et al., 2013; Ping et al., 2014; Zhang et al., 2016; Zhao et al., 2017; Zheng et al., 2013; Zhong et al., 2008). In addition, m⁶A is involved in multiple RNA metabolic processes, including mRNA expression (Dominissini et al., 2012), translation efficiency (Wang et al., 2015), alternative splicing (Liu et al., 2015), and degradation (Wang et al., 2013).

High-throughput sequencing analysis of m⁶A targeted antibodies in Arabidopsis showed that more than 70% m⁶A peak was detected in Arabidopsis Can-0 and Hen-16, indicating that m⁶A is a highly conservative modification of mRNA in plants (Luo et al., 2014). Subsequently, differential m⁶A methylation patterns among three organs of Arabidopsis were analyzed using transcriptome wide high-throughput deep m⁶A-seq. The results showed that over 80% of m⁶A modified transcripts were identified in leaves, flowers, and roots of Arabidopsis (Wan et al., 2015). The study of m⁶A patterns in the chloroplast and mitochondrial transcriptome of Arabidopsis shows that the m⁶A motif that undergoes RNA methylation is highly conserved between the nucleus and the organelle transcriptome (Wang et al., 2017). The m⁶A map of sea buckthorn (Hippophae rhamnoides Linn.) transcriptome identified 13,287 different m⁶A peaks between leaf under drought and in control treatment. In addition, it is reported that after the m⁶A modification disorder, the leaves and roots of plants have serious dysplasia, and the time of reproductive development changes (Zhang et al., 2021). These findings indicate that m⁶A has a regulatory role in plant gene expression. In recent years, with the deepening of research on m⁶A writers, erasers, and readers in plants, it has been found that m⁶A methyltransferase in Arabidopsis functions by forming complexes, and the accessory subunits FIP37 and VIR act as key subunits to maintain the function of m⁶A methyltransferase complexes (Shen, 2023). At the same time, many m⁶A 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., 2016). ALKBH10B is an RNA m⁶A demethylase that regulates inflorescence transformation in Arabidopsis (Duan et al., 2017). The m⁶A reader ECT2 recognizes the site of m⁶A modification and regulates the morphogenesis of Arabidopsis trichomes (Wei et al., 2018). RNA-binding protein FLK is a novel mRNA m⁶A reader protein that directly binds to the m⁶A 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., 2023). M⁶A-related proteins were also identified in tobacco. In tobacco infected with Tobacco mosaic virus (TMV), the gene expression level of m⁶A demethylase ALKBH5 increased, and the results showed that TMV reduced the m6A level (Li et al., 2018). Until now, m⁶A modification has rarely been investigated in woody plants, so it is necessary to study m⁶A 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., 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., 2015). Advances in biotechnology have made it possible to study posttranscriptional mRNA modifications, especially m⁶A modifications, which have become the focus of many studies. Given previous findings on m⁶A modification, most of the neurological diseases, embryonic developmental delays, and tumorigenesis caused by m⁶A 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⁶A modification.

Previous studies have shown that phytoplasma-infected Paulownia seedlings restored to normal morphology after treatment with 60 mg·L⁻¹ methyl methanesulfonate (MMS), and no phytoplasma could be detected (Wang et al., 2018). To investigate the changes in m⁶A 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⁻¹ MMS (PFIM60) were collected, and the m⁶A modification map in Paulownia was obtained by m⁶A-Seq, and the differentially expressed genes with m⁶A modification were further screened, by combining with transcriptome analysis. In view of the function of m⁶A in regulating alternative splicing, the genes modified by m⁶A 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⁶A modification in Paulownia after phytoplasma infection will lay a solid scientific foundation for the discovery of PaWB-related pathogenesis and epigenetic regulation mechanisms.

MATERIALS AND METHODS 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⁻¹ MMS (PFIM60 samples) (Figure 1b) or 0 mg·L⁻¹ MMS (PFI samples) (Figure 1a) (Wang et al., 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. 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.

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⁶A-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⁻¹). The mixture was incubated with protein-A beads and eluted with elution buffer (1 × IP buffer and 6.7-mM m⁶A). The eluted RNA was precipitated with 75% ethanol. The eluted m⁶A-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., 2010). We constructed four m⁶A-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.

Alignment of reads to the reference genome and visualization of m⁶A peaks

Cutadapt (Martin, 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) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Bowtie (Langmead & Salzberg, 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., 2014) to identify m⁶A peaks (in BED or BAM format) for visualization using the UCSC Genome Browser or Integrative Genomics Viewer software (http://www.igv.org/).

Transcription level analysis

StringTie (Pertea et al., 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., 2010).

Discernment of m⁶A topological patterns

MEME (Bailey et al., 2009) and HOMER (Heinz et al., 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., 2015).

M⁶A MeRIP-reverse transcription PCR

The m⁶A methylated RNA immunoprecipitation (MeRIP) assay was adapted from a previously reported protocol (Dominissini et al., 2013). Briefly, approximately 300-mg total RNA was extracted and isolated using an mRNA screening kit. Each biological replicate for m⁶A-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⁶A-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⁶A-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⁶A negative fragments that were not captured by the anti-m⁶A antibody. The m⁶A-Dynabeads–RNA complexes were resuspended in 500-mL m⁶A 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 S1.

RESULTS The m⁶A modification landscape of P. fortunei

PaWB-infected P. fortunei (PFI) seedlings and PFI seedlings treated with 60 mg·L⁻¹ MMS (PFIM60) were used for m⁶A-seq, RNA-seq (total fragmented RNA as the control for m⁶A-seq), and transcriptome analysis, with two replicates for each. Heatmap correlation coefficient analysis between the biological replicates confirmed their high repeatability (Figure 2). We obtained 44–56 and 73–83 million reads for the RNA-seq and m⁶A-seq libraries, respectively (Table S2), of which 23–28 and 47–50 million were uniquely aligned to the P. fortunei reference genome (Table S3). Almost 94% of reads aligned to the reference genome were located in exons; the others were located in introns or intergenic regions (Figure S1 and Table S4).

FIGURE 2. Heatmap correlation coefficient analysis between biological replicates.

A m⁶A modification map of P. fortunei was constructed, showing the distribution of m⁶A on chromosomes, as well as the distribution of m⁶A peak positions in PFI and PFIM60 on chromosomes (Figure 3). The number of m⁶A-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⁶A 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. M6A modification map of Paulownia fortunei. The outermost circle is the distribution of chromosome, the second circle is the distribution of PFI m6A peak positions on chromosomes, and the third circle is the distribution of PFIM60 m6A 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⁶A modification of PFI was slightly higher than that of PFIM60, which indicated that PaWB phytoplasma infection caused the change of m⁶A modification of Paulownia (Dataset S1). Notably, there was greater enrichment of peaks at the transcription start and stop sites of genes in both samples (Figure 4a), suggesting that m⁶A modification may play a key role in regulating mRNA degradation and microRNA processing in Paulownia. The differences in m⁶A peaks observed between PFI and PFIM60 were compared and analyzed. The results showed that there were 306 peaks of 319 genes with reduced m⁶A modification in PFIM60 compared with their m⁶A 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⁶A enrichment of all the genes in the P. fortunei genome (Figure 4b).

FIGURE 4. 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.

Differential modified genes between PFI and PFIM60

To determine the sequence specificity of the hypermethylated m⁶A peaks, we performed motif prediction based on differential peak analysis using HOMER and compared the predicted motif with microRNA sequences in miRBase (1e⁻²⁰ < p < 1e⁻¹⁰). 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 5), which is similar to the UGUAYY (Y = C/U) motif that is specific to m⁶A 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⁶A-related enzymes, which in turn affect gene expression.

FIGURE 5. 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⁶A 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 6 and Dataset 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. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of genes with reduced m6A modification after 60 mg·L−1 methyl methanesulfonate treatment.

Extent of m⁶A 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 7a). KEGG pathway enrichment analysis showed that DEGs were mainly involved in plant circadian rhythm, phenylpropanoid biosynthesis, flavonoid biosynthesis, and glycerolipid metabolism pathways (Figure 7b).

FIGURE 7. 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⁶A modification versus gene transcript levels in the PFI and PFIM60 samples, we correlated m⁶A-seq data with transcriptome sequences. In the transcriptome, differences in transcript abundance were categorized as upregulation and downregulation, while in m⁶A-seq data, changes in peak abundance indicated m⁶A 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⁶A modification levels and gene expression were identified: (1) both m⁶A modification levels and transcription levels upregulated (133 genes); (2) m⁶A modification levels upregulated and transcription levels downregulated (58 genes); (3) m⁶A modification levels downregulated and transcription levels upregulated (44 genes); and (4) m⁶A modification levels and transcription levels both downregulated (80 genes) (Figure 8 and Dataset 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 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 9b).

FIGURE 8. Four relationships were detected between m6A modification levels and gene expression.

FIGURE 9. Functional annotations and pathways of differential peaks and differentially expressed genes. (a) Statistics of gene ontology (GO) enrichment. (b) Statistics of pathway enrichment.

M⁶A modification affects alternative splicing

Alternative splicing is involved in many physiological processes, including responses to abiotic and biotic stresses (Barbazuk et al., 2008). M⁶A modification can affect alternative splicing by promoting retention of exons at the transcriptional level; therefore, we jointly analyzed m⁶A 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⁶A modification was higher in PFIM60 than in PFI, alternative splicing was associated with m⁶A 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⁶A modification increased, the bands of two genes were weakened compared with those for PFI, suggesting that the degree of m⁶A modification affected alternative splicing events (Figures 10 and S2). Primer sequences are shown in Table S5.

FIGURE 10. RT-PCR validation of F-box and MSH5 in PFI and PFIM60, respectively. 18s rRNA is the internal reference gene.

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., 2007), WUSCHEL (WUS) (Lenhard et al., 2002), Arabidopsis TOPLESS-like transcriptional corepressor (ASP1), SHOOT MERISTEMLESS (STM), ARABIDOPSIS RESPONSE REGULATOR 5 (ARR5), ARR6, ARR7, and ARR15 (Leibfried et al., 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., 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⁶A modification and expression levels after phytoplasma infection.

We performed m⁶A 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 11).

FIGURE 11. Validation of m6A methylated RNA immunoprecipitation quantitative reverse transcription PCR. In PFI samples infected with phytoplasma, the m6A modification degree of STM and CLV2 was higher than that of healthy PFIM60 samples.

STM and CLV2 maintain the dynamic balance of shoot apical meristem regulated by m⁶A modification

Conserved domains were identified in the proteins encoded by the Paulownia homologs of STM and CLV2 (Table 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 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⁶A-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	m6A	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	—	—	—	PLN00113 superfamily	—	—	—
	Paulownia_LG2G000076	Yes	Down	Down	PLN00113 superfamily	—	—	—

DISCUSSION M⁶A modification of P. fortunei before and after PaWB

The m⁶A 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., 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⁶A have been conducted in the model plants such as Arabidopsis thaliana and tobacco, but there are less reports on the map of m⁶A 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⁶A modification may play an important role in the infection of PaWB phytoplasma.

M⁶A 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⁶A-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⁶A 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., 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⁶A methylase FIP37 affected the expression of WUS and STM and regulated the development of shoot apical meristem (Shen et al., 2016). These results suggested that m⁶A methylase may also be involved in regulating STM expression in Paulownia. However, this possibility needs to be further confirmed by molecular biology experiments.

The relationship between m⁶A 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, 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⁶A modification and alternative splicing analysis, it was found that when m⁶A modification increased, alternative splicing was associated with two m⁶A 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., 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., 2017). The SCF complex can degrade JA, resulting in the release of MYC2 to express genes responsive to JA (Liu et al., 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, 2013; Spampinato, 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, 2020; Gonzalez & Spampinato, 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.

CONCLUSIONS

In this research, m⁶A-seq was applied to determine methylation levels in phytoplasma-infected Paulownia. Results indicated that m⁶A modifications differed from those previously reported in model plants. We speculated that the Paulownia genome had a novel m⁶A modification motif. DEGs associated with PaWB were detected by transcriptomic analysis. The effect of m⁶A 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.

ACKNOWLEDGMENTS

This work was supported by the Academic Scientist fund for Zhongyuan scholars of Henan Province (Grant No. 2018[99]) and Leading Talents of Zhongyuan Science and Technology Innovation (224200510010). We also thank Margaret Biswas, PhD, from Liwen Bianji (Edanz) (www.liwenbianji.cn/) for editing the English text of a draft of this 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 Supporting Information for this article.

DATA AVAILABILITY STATEMENT

The datasets generated and/or analyzed during the current study are available in the [NCBI] repository (SRA accession: PRJNA783736 [m⁶A-Seq data] and PRJNA624264 [RNA-Seq data]). The P. fortunei used in this study came from the Forest Biotechnology Laboratory of the Institute of Paulownia, Henan Agricultural University, Zhengzhou, China, and permission has been granted to collect P. fortunei.