The m6A 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 m6A-seq, RNA-seq (total fragmented RNA as the control for m6A-seq), and transcriptome analysis, with two replicates for each. Heatmap correlation coefficient analysis between the biological replicates confirmed their high repeatability (Figure 1). We obtained 44–56 and 73–83 million reads for the RNA-seq and m6A-seq libraries, respectively (Supplementary Materials: Tables S1), and 23–28 and 47–50 million of them uniquely aligned to the P. fortunei reference genome (Supplementary Materials: Tables S2). 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 (Supplementary Materials: Figure S1, Supplementary Materials: Tables S3).
We constructed a m6A methylation map of P. fortunei that shows the distribution of m6A on the chromosomes, and the distribution of the PFI and PFIM60 m6A peak positions on the chromosomes (Figure 2). The number of m6A-modified genes on the chromosomes ranged from 472 to 1,014; chromosome 9 (1,014) and chromosome 2 (1,002) had more modified genes than the other chromosomes. The highest number of m6A peaks (1,506) was on chromosome 2 and the lowest number (702) was on chromosome 14. The numbers of methylated genes and peaks on the chromosomes were similar in PFI and PFIM60.
We obtained 6,606 high-confidence m6A peaks located in 13,837 and 13,505 genes in PFI and PFIM60, respectively, using ChIPseeker; 12,554 of the genes were common to both samples (Supplementary Materials: Tables S4). Notably, there was greater enrichment of peaks at the transcription start and stop sites of genes in both samples (Figure 3a), suggesting m6A modification may play a key role in regulating mRNA degradation and microRNA processing in Paulownia. We compared and analyzed the differences in the observed m6A peaks between PFI and PFIM60. The results showed there were 306 genes with reduced m6A methylation in PFIM60 compared with their m6A methylation in PFI; 36.59% of them were located in 3′ UTRs, 17.41% were in 5′ UTRs, 21.58% were in first exons, and 24.15% were in other exons. This distribution is consistent with the m6A enrichment of all the genes in the P. fortunei genome (Figure 3b).
Differential modified genes between PFI and PFIM60
To determine the sequence specificity of the hypermethylated m6A peaks, we used HOMER for motif prediction based on the difference peak analysis, and compared the predicted motif to the 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), which was identified in Paulownia, has not yet been confirmed in mammals and yeast. The most significantly enriched sequence was UUGUUUUGUACU (Figure 4), which is similar to the UGUAYY (Y = C/U) motif that is specific to m6A modifications in tomato, rice, corn, and other plants, indicating that the prediction results have a high degree of confidence. These motif sequences are recognized and bound by m6A-related enzymes, which in turn affect gene expression.
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed on the peak-related genes with weakened m6A methylation modification detected in the PFIM60 samples. The results showed that the auxin biosynthesis, hormone signal transduction, and MAPK signal pathways were some of the most enriched pathways (Figure 5; Supplementary Materials: Tables S5). The auxin biosynthesis pathway involves three genes, Paulownia_LG2G000124, Paulownia_LG2G000125, and Paulownia_LG2G000128. The plant hormone (auxin, abscisic acid, and ethylene) signal transduction pathways involve eight genes, 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.
Extent of m6A methylation 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 the 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 5,899 DEGs were detected between the PFI and PFIM60 libraries; 3,130 were up-regulated and 2,769 were down-regulated in PFIM60. Gene ontology (GO) enrichment analysis showed that the DEGs were involved mainly in membrane, cell wall, transcriptional regulation using DNA as a template, secondary metabolite biosynthetic process, DNA binding transcription factor activity, and oxidoreductase activity (Figure 6a). The KEGG pathway enrichment analysis showed that the DEGs were involved mainly in plant circadian rhythm, phenylpropanoid biosynthesis, flavonoid biosynthesis, and glycerolipid metabolism pathways (Figure 6b).
To examine the extent of m6A methylation versus gene transcript levels in the PFI and PFIM60 samples, we correlated the m6A-seq data with the transcriptome sequences. In the transcriptome, differences in transcript abundance were categorized as up- and down-regulation, whereas in the m6A-seq data changes in peak abundance indicated the m6A methylation levels of up-regulated genes. A total of 315 differentially methylated genes were predicted to be significantly differentially expressed at the transcriptome level. Four types of relationships between the m6A methylation levels and gene expression were identified: 1) m6A modification levels and transcription levels both up-regulated (133 genes); 2) m6A modification levels up-regulated and transcription levels down-regulated (58 genes); 3) m6A modification levels down-regulated and transcription levels up-regulated (44 genes); and 4) m6A modification levels and transcription levels both down-regulated (80 genes) (Figure 7; Supplemental Data Table S6). The GO enrichment analysis of these differentially methylated genes showed 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 function were enriched under the molecular function category (Figure 8a). The KEGG pathway analysis showed that these differentially methylated genes were involved mainly in signal pathways such as ABC transport, diterpenoid biosynthesis, plant hormone signal transduction, and phenylpropanoid biosynthesis (Figure 8b).
m6A modification affects alternative splicing
Alternative splicing (AS) is involved in many physiological processes, including responses to abiotic and biotic stresses (Barbazuk et al. 2008). m6A modification can affect AS by promoting the retention of exons at the transcriptional level; therefore, we jointly analyzed m6A modification and AS. In PFIM60, 282 genes had significant exon skipping-type AS differences and 84 genes had significant mutually exclusive exon-type AS differences compared with the corresponding genes in PFI. When the degree of m6A methylation increased in PFIM60 compared with that in PFI, we found the AS of three genes was correlated with m6A modification, namely F-box (Paulownia_LG17G000760) and MSH5 (Paulownia_LG8G001160). To verify the AS of these two genes in PFI and PFIM60, we performed an RT-PCR analysis. We found that for PFIM60, when the degree of m6A methylation increased, the bands of two genes were weakened compared with those for PFI, indicating that the degree of m6A methylation affected AS events (Figure 9; Supplemental Data Figure S2).
Identification of potential candidate genes associated with PaWB
In the process of invading a host, inherent pathogen-associated molecular patterns called PAMPs or damage-related molecular patterns of the phytoplasma are recognized by pattern recognition receptors in the host defense system, which induces PAMP-triggered immunity and activates the downstream defense pathway. To resist PAMP-triggered immunity, phytoplasmas release effectors as signal molecules to regulate signal transduction and disrupt the host defense response. 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 pathway: Paulownia14027, which encoded SHORT-ROOT (LSH4); Paulownia26350, which encodes LIGHT-DEPENDENT SHORT HYPOCOTYLS 4 (LSH4); Paulownia_LG17G000277, which encodes histidine-containing phosphotransfer 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 (TPL)-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 that are involved in the WUSCHEL–CLAVATA negative feedback loop, which is the basic mechanism for maintaining stem cells in shoot apical meristem. We found Paulownia homologs of CLV2 (Paulownia_LG2G000076) and STM (Paulownia_LG15G000976) by mapping the Arabidopsis genes to the Paulownia reference genome using a local BLAST program. These two genes showed significant differences in m6A methylation and expression levels after phytoplasma infection.
We performed a m6A methylated RNA immunoprecipitation quantitative reverse transcription PCR (MeRIP-qRT-PCR) verification on these two candidate genes. The results were consistent with the sequencing results, which verified the accuracy of the sequencing results (Figure 10).
STM and CLV2 maintain the dynamic balance of shoot apical meristem regulated by m6A methylation 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 bind to the CLV3 promoter and activate its transcription, and CLV3, CLV1, and CLV2 together maintain the number of undifferentiated stem cells in the shoot apical meristem tissue and regulate the production of normal stem ends. The transcriptome data analysis showed that STM was more highly expressed in PFI than it was in PFIM60, and the m6A-seq data analysis showed that CLV1, CLV2, and CLV3 had higher methylation levels in PFI than in PFIM60.
Table 1
conserved domain analysis for two protein related to PaWB
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
|
-
|
-
|
-
|