Transcriptome and Small RNA Sequencing Analysis Revealed Roles of PaWB-Related miRNAs and Genes in Paulownia fortunei

Paulownia witches’ broom (PaWB) is an epidemic disease caused by phytoplasmas infection, which is responsible for large production and economic losses. The study of PaWB has made significant progress, but the specific molecular mechanisms associated with PaWB remain unclear. To clearly know the gene expression profiles of plantlets infected with phytoplasmas, in this study, we used high-throughput sequencing technology to generate an integrated analysis of the transcriptome and microRNAs (miRNAs) of Paulownia fortunei (seem.) Hemsl. plantlets, and to obtain a comprehensive resource for the relationship between vital miRNA-target gene pairs and PaWB. A total of 756 genes, and 45 conserved and 22 new miRNAs were identified associated with PaWB. In addition, 635 target genes were predicted for the 67 DERs (Differentially expressed miRNAs). An interaction network of these miRNAs and their target genes was constructed. Gene ontology (GO) and KEGG pathway analysis of these target genes indicated that genes encoding transcription factors (TFs), including auxin response factors (ARF), WRKY, NAC (NAM, ATAF1/2 and CUC2), and MYB (v-myb avian myeloblastosis viral oncogene homolog), and genes encoding superoxide dismutase (SOD), as well as alternative splicing were related directly or indirectly to PaWB. Our results shed light on the possible roles of genes and miRNAs in PaWB-infected plantlets, which will enhance the understanding of the PaWB mechanism in Paulownia plants.


Introduction
As an important member in Scrophulariaceae family, the genus Paulownia contains some fast growing tree species [1].Up to now, Paulownia has been introduced to many countries and areas, such as Japan, Southeast Asia, Australia, and Brazil [2].Taking advantage of its excellent characteristics such as high-quality, even texture, strong resistance to water stress and saline alkali, Paulownia was widely utilized in building construction, furniture, musical instruments, various craft products, and even for soil conservation, windbreak, and sand fixation [3,4].
Nevertheless, the production of Paulownia would suffer severe or even fatal loss when it is infected with phytoplasmas [5].Phytoplasmas are single-cell prokaryotic organisms without a cell wall, which can infect more than 1000 plant species such as jujube [6], Paulownia [7], mulberry [8], peanut [9] and sweet potato [10].As obligate parasites, they live and multiply not only in the intestinal tract, lymph, salivary glands and other tissues of insects, but also in the phloem cells of plants [11,12].The host plants of phytoplasmas exhibit mainly stem section shortening, leaf chlorosis, phyllody, and witches' broom [13,14].The symptoms in phytoplasmas-infected plants might be caused by nutrient consumption, which was directly related to the activity of effector proteins secreted by phytoplasmas [15].Since phytoplasmas have no cell wall, they can directly interact with their host plants by exposing their membrane proteins or effectors to plant cytoplasm [16].Previous research indicated that phytoplasmas can regulate the gene expression levels by secreting the effector proteins TENGU, SAP11, and SAP54 to induce typical symptoms in host plants [17,18].In order to defend pathogen invasion, the corresponding immune responses were subsequently launched in host plants.So far, two immune responses against pathogen have been reported in host plants, i.e., effector-triggered immunity (ETI) and pathogen-associated molecular pattern (PAMP)-triggered immunity (PTI) [19].During the past few decades, some proteins, metabolites, miRNAs and genes related to plant hormone signal transduction, secondary metabolism, and photosynthesis have been found to be involved in Paulownia witches' broom (PaWB) [20][21][22][23][24], while the interaction between phytoplasmas and infected plants is not yet clearly known.
The gene expression network is usually regulated by multiple factors at transcriptional levels.Among all these regulatory factors, microRNAs (miRNAs, a kind of small noncoding RNAs (21-24 (nucleotides) nt) that are derived from pre-miRNA (50-350 nt)) play vital roles in gene expression regulation.Generally, miRNAs inhibit the translation process by binding to the 3 untranslated regions of their target genes [25,26].In addition, miRNAs are also reported to influence many biological processes in biotic/abiotic stress responses [27,28].For instance, miRNA393 was found to help the host plant cells to recognize the characteristics of pathogens, and then triggered a series of defensive responses such as rapid gene expression changes, and phytohormone and metabolite induction [29,30].MiR156 and miR164 were induced in response to virus infection in Nicotiana tabacum Linn.[31] and Arabidopsis thaliana Heyn [32].At the same time, miR156, miR159 and miR172 were differentially expressed in phytoplasmas infection, which played critical roles in Jujube witches'-broom (JWB) disease [6].
Advances in high-throughput sequencing technologies have made the integrated analysis on small RNAs and transcriptome sequencing easier.Moreover, the availability of the Paulownia fortunei (P.fortunei) genome information has enhanced the accuracy of gene annotations in Paulownia.In this study, we combined transcriptome and small RNA (sRNA) sequencings to conduct a comprehensive analysis on the immunity mechanisms of PaWB in P. fortunei, and attempt to elucidate a miRNA-target gene interaction network related to PaWB.Overall, our results provide a platform to better understand the interaction of pathogen-host in Paulownia and other tree species.

Plant Materials
Healthy P. fortunei (PF) and PaWB-infected P. fortunei (PFI) plantlets were obtained from the Institute of Paulownia, Henan Agricultural University, China (34 • N, 113 • E).The tissue culture plants of PF and PFI were grown in 1/2 Murashige and Skoog (1/2 MS) medium in 350 mL plastic culture flasks for 30 days.The cultivation conditions in the phytotron were as follows: temperature, 25 ± 2 • C, light intensity, 130 µL•m −2 •s −1 , and photoperiod, 16 h light/8 h dark.The apical buds of PF and PFI (30 apical buds were mixed into one sample) were collected when they grew to 1.5 cm.The mixed samples were placed immediately in liquid nitrogen and stored at −80 • C for RNA extraction.Three biological replicates were used for transcriptome and sRNA sequencing of PF and PFI.Principal component analysis (PCA) [33] was utilized to perform biological replicates analysis.

Transcriptome Sequencing and Data Analysis
Total RNA of PF and PFI was isolated using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and the ribosomal RNAs were removed.RNA pollution and disintegration were evaluated by 1% agarose gel electrophoresis.The integrity and purity of the RNA samples were detected using a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) and RNA 6000 Nano LabChip Kit (Agilent).After purification, the poly (A)− or poly (A)+ RNAs were cleaved into fractions by fragmentation buffer.The cDNAs were synthesized according to the previous published report [34].Then, paired-end sequencing of PF and PFI was performed on an Illumina HiSeq 4000 platform (LC-BIO, Hangzhou, China) following the manufacturer's instructions.Low-quality reads that contained poly (N) and adapter sequences were removed from the raw reads using the Illumina Pipeline to obtain high-quality valid reads.The valid reads were mapped to the reference genome (http://paulownia.genomics.cn/page/species/index.jsp) by SOAP2 software (Version 2.21, http://soap.genomics.org.cn/).Then, the mapped reads were searched against NCBI's (National Center for Biotechnology Information) non-redundant protein sequence (NR) database (http://www.ncbi.nlm.nih.gov)[35], Swiss-Prot protein database (http://www.uniprot.org/),Clusters of Orthologous Groups (COG) database (http://www.ncbi.nlm.nih.gov/COG),gene ontology (GO) database (http://www.geneontology.org/), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (http://www.genome.jp/kegg)by BLASTX (e value < 10 -5 ) for functional annotation.
In order to identify the reads that were related to PaWB, the PF and PFI libraries were compared according to Audic et al., and fragments per kilo base of transcript per million mapped reads (FPKM) were used to measure the expression levels of the transcript [36].The critical p-values in multi-hypothesis testing and analysis were confirmed by the false discovery rate (FDR).Differentially expressed genes (DEGs) were identified by |log2 fold change| ≥ 1 and FDR ≤ 0.001 [37].The pathways were enriched by in-house software.

Construction, Sequencing, and Identification of sRNAs
Total RNA samples were used for the construction of PF and PFI sRNA libraries.The sRNA fractions were segregated by polyacrylamide gel electrophoresis (PAGE).The 18-30 nt fractions were used to produce the libraries, and T4 RNA ligase was used to connect 5 and 3 adapters.After a series of reverse transcription and PCR processes, single-end sequencing was performed on the PCR-amplified products using an Illumina HiSeq 2500 platform at the LC-BIO (Hangzhou, China).The adapter, pollution, low-quality reads, and the reads <18 nt were removed from the raw reads.The valid reads were mapped to the genome using the SOAP2 software [38].The mapped reads were further filtered using a proprietary pipeline script, ACGT101-miR v4.2 (LC Science, Houston, TX, USA), to remove rRNAs, tRNAs, and snRNAs.The remaining unique reads were matched to miRBase Release 21.0 (http://www.mirbase.org/)using BLAST to identify known miRNAs (aligned sequence should contain <2 mismatches).The remaining reads that cannot match to miRBase were examined using MIREAP (Version: 0.2, http://sourceforge.net/projects/mireap/)and RNAfold (Version: 3.8, http://rna.tbi.univie.ac.at/cgibin/RNAfold.cgi)[39] to predict the typical stem-loop hairpin structures to obtain candidate novel miRNAs.In addition, the criteria for novel miRNAs screening were formulated by Meyers et al. [40].

DERs Identification and Target Genes Prediction
DERs were predicted as Han et al. described [41].To better understand the molecular function of the miRNAs that related to PaWB, TargetFinder software (Copyright: 2007-2010, http://jcclab.science.oregonstate.edu/node/view/56334,computational target prediction algorithms) was used to identify the miRNA binding sites in the mRNAs to predict the target gene of the DERs.The BLASTX algorithm and GO analysis were used to search and annotate the target genes in the NCBI database (http:www.ncbi.nih.gov/gene/data/GO).

Construction of the miRNA-Target Gene Interaction Network
A miRNA-target gene interaction network was constructed based on the PaWB-related DERs and their target genes.The relationship between DERs and genes were identified using the STRING database v10.5 (https://string-db.org/).Finally, the interaction network of DERs and genes was constructed by Cytoscape software v2.8.3 (http://chianti.ucsd.edu/Cyto-2_8_3/)[42,43].

qRT (Quantificational Real-Time Polymerase)-PCR Validation
PCR was performed to verify the high-throughput sequencing and target gene predictions results.The RNA samples used for the qRT-PCR were the same as the samples that used for sequencing.The cDNAs of the miRNAs and mRNAs were synthesized with a SYBR ® PrimeScriptTM miRNA RT-PCR Kit (TaKaRa, Dalian, China) and a Superscript III First-Strand Synthesis system.The primers of miRNAs were designed basing on the corresponding mature miRNA sequences, and the primers of the target genes were designed using Primer 5 program (Applied Biosystems, Foster City, CA, USA).The qRT-PCR was performed with a SYBR ® Premix Ex TaqTM Kit (TaKaRa, Dalian, China) and a SYBR ® PrimeScriptTM miRNA RT-PCR Kit (TaKaRa, Dalian, China).The parameters in the PCR amplification program were set as: 95 • C for 1 min; 40 cycles of 95 • C for 10 s, and 55 • C for 15 s.18S rRNA and U6 were selected as the internal reference genes for the miRNAs and mRNAs, respectively.The primers used for the qRT-PCRs are listed in Table 1.
Table 1.Primers of quantitative PCR analysis.

Transcriptome Sequencing Analysis
To detect the genes that were related to PaWB, transcriptome sequencing (with three biological replicates) of PF and PFI libraries was performed.The PCA results of the three biological replicates showed the data have good reproducibility (Figure 1A).A total of 131,903,380 and 141,209,455 raw reads were obtained from the PF and PFI mRNA libraries, respectively.After filtering out low-quality reads, 118,789,530 and 134,197,596 valid reads were obtained (Table 2).The valid reads that can map to the P. fortunei reference genome were shown in Table 3.To investigate the expression levels of all the genes, a FPKM boxplot was drawn, and this also provided a preliminarily indication of the repeatability of the samples (Figure S1).In order to annotate the gene functions, NR, Swiss-Prot, KEGG, COG, and GO databases were searched (E < 10 −5 ).According to the DEG criteria, 756 DEGs that may be related to PaWB were detected (Table S1).To better understand the functions and classifications of the 756 DEGs, the GO and KEGG databases were used to perform gene annotation.The GO functional analysis indicated that most of the DEGs were mainly involved in these terms, such as the oxidation-reduction process, nucleus and protein binding.According to the KEGG pathway enrichment analysis, most of the DEGs were enriched in the glycolysis/gluconeogenesis, glycerolipid metabolism, and flavonoid biosynthesis pathways (Figures S2-S4).

miRNA Identification
The PF and PFI sRNA libraries were constructed to study the miRNA expression profiles in response to PaWB.The PCA of the three biological replicates was showed in Figure 1B.More than 10 million raw reads were obtained from each replicate.After filteration, 8,180,411 and 7,109,787 valid reads from PF and PFI were obtained, respectively (Table S2).The lengths of valid reads ranged from 18 nt to 24 nt (Figure 3), which is consistent with a previous report that suggested 24 nt and 21 nt sRNAs may play crucial roles in PaWB resistance.A total of 493 miRNAs in the PF and PFI libraries were identified (Table S3); 321 miRNAs belonging to 51 families were known miRNAs.Among them, miR156, miR159, and miR166 families occupied high percentages (Table S4).The remaining 172 miRNAs that cannot map to miRBase sequences were identified as new candidate miRNAs.Most of the new miRNAs were found in lower abundance than the conserved miRNAs, which is consistent with previous studies [46].

Functional Analysis of the PaWB-Related miRNAs and Their Target Genes
To identify the miRNAs which were involved in PaWB, 67 DERs (22 novel and 45 conserved) that satisfied two conditions, Fold Change (FC) ≥ 1.0 or ≤ −1.0 and p-values ≤ 0.05) were detected for analysis (Table S5).According to the heatmap, the vast majority of the DERs that belonged to the same family exhibited the same expression patterns.For instance, five, three, and four members of the miR319, miR156, and miR169 families were upregulated in the PFI library, and three, three, and two members of the miR395, miR6300 and miR171 families were downregulated.Among the new miRNAs, seven were upregulated and 15 were downregulated (Figure 4).In addition, 569 target genes of 39 conserved miRNAs were identified (Table S6).Furthermore, 66 target genes of seven new miRNAs were also predicted.Ten of the miRNAs targeted a single gene, while the remaining miRNAs targeted more than two genes.Moreover, Pf-miR8767 had 123 predicted target genes.
The 635 predicted target genes were annotated by GO functional classification to further determine whether they played vital roles in PaWB (Figure 5).The GO functional analysis indicated that transcription and transcription regulation processes were the most abundant terms under the biological process category.Under the cellular component category, the most abundant terms were nucleus and plasma membrane.Under the molecular function category, the most abundant term was protein binding.The functional analysis of these target genes indicated that TFs such as MYB3, MYB5, ARFs, NAC, and WRKY may be related to PaWB.In addition, GO and KEGG pathway enrichment were performed to further understand the biological functions of target genes.The results showed that some critical biological pathways may be involved in PaWB response (Figure 6A,B).GO terms related to disease resistance (GO: 0003700 transcription factor activity, sequence-specific DNA binding; GO: 0004364 glutathione transferase activity) and phytohormone (GO: 0009737 response to abscisic acid; GO: 0009741 response to brassinosteroid; GO: 0009734 auxin-activated signaling pathway) were enriched.The KEGG pathway analysis indicated that the 735 target genes were mainly involved in plant-pathogen interaction, amino sugar and nucleotide sugar metabolism, and endocytosis.For example, the target genes of pf-miR858, miR160, and miRNA164 families encode MYB-like, ARF, and NAC TFs (which were reported to be involved in plant defense, and growth and development processes).In this study, miR156, miR159, miR397, and miR403 families were significantly upregulated in the PFI library.Taken together, these results demonstrate that the DERs and their target genes may play vital roles in response to PaWB.

The miRNA-Target Gene Interaction Network That Related to PaWB
To further clarify the regulatory functions of the miRNAs and their target genes that involved in PaWB, a miRNA-target gene interaction network was constructed with Cytoscape.The network suggested that these DERs participated in multiple regulation process including morphogenesis, plant defense and plant hormones signal transduction.For instance, pf-miR858b, pf-miR2628 and pf-miR164b-5p mainly involved in regulation of plant defense.Among them, pf-miR858b played an essential role in the regulation layer of secondary cell wall biosynthesis.Pf-miR2628 might directly mediate the plant defense by regulating the expression of WRKY4 (or indirectly regulation of WRKY8).Pf-miR164b-5p might involve in plant defense process via activating the cascade response of protein kinases, which was essential for pathogen defense.Mildew resistance locus o 6 (MLO6) (target of Pf-miR164b-5p) mainly interacted with several proteins that participated in transmitting the signal to protein kinase, such as AT1G61360, AT1G74360, WAKL10, AT1G18390 and AT5G25930.Pf-miR5046 and pf-miR160h mainly participated in the signal transduction process of auxin (AUX) and abscisic acid (ABA), which were essential in plant development, and they seemed to have a more complicated mechanism to the hormone in the regulation layer.Furthermore, pf-miR2624 that targeted (SPL8) was mainly involved in the morphogenesis of leaf and floral organs.In summary, the network indicated that DERs may be involved in multiple developmental processes of plants in different ways to respond to PaWB-infection (Figure 7).

Verification of the Expression Patterns of miRNAs and Their Target Genes
To verify the reliability of the high-throughput sequencing results, 11 differentially expressed mRNAs, eight target genes and their corresponding DERs were randomly chosen for qRT-PCR verification.The results indicated that the expression patterns of the DERs and differentially expressed mRNAs were similar in high-throughput sequencing.Further analysis suggested that expression patterns of four miRNA-mRNA pairs (pf-miR160h-PAU000742.1,pf-miR858b-PAU027187.1, pf-miR169a-3p-PAU012221.1, and pf-miR172e-PAU025602.1) were negatively correlated, and the expression patterns of the other four miRNAs-mRNA pairs (pf-miR164b-5p-PAU003364.1, pf-miR167h-PAU012184.1, pf-miR156p-PAU012247.1, and pf-miR5048-PAU029039.1) were positively correlated (Figure 8).Similar regulation patterns of miRNA-mRNA pairs have been reported in cotton [47].These results indicate that the transcriptome data are sufficient to estimate miRNA-target gene interactions in phytoplasmas infected P. fortunei.

Discussion
The availability of reference genome sequences makes the annotation of important genes, proteins, lncRNAs (LncRNAs (long non-coding RNAs) are defined as having more than 200 nucleotides and little protein-coding potential, which are an important class of pervasive genes involved in a variety of biological functions [48,49].),metabolites, and miRNAs more precise compared with the transcriptome background, as well as providing meaningful insights into genome reorganization, gene evolution and comparative genomic analyses between multiple species [50].Thus, based on the reference genome, transcriptome and miRNAs were performed in phytoplasama infected plantlets, and obtained several key genes related to PaWB.Previous studies about genes related to PaWB have been found, but the mechanisms were still unclear.To get a deeper understanding of the interaction between Paulownia and phytoplasmas, we compared the data obtained in this study with others [51,52].The number and species of miRNAs or genes in the present study were different.The probable reasons for this difference are as follows: three biological replicates were performed in this study, making the data more reliable; the plant materials in previous reports were treated with reagent, which might lead to the changes in miRNAs and gene expression; and high-throughput sequencing was performed with the genome of P. fortunei as the reference.

PaWB Responsive DER-Target Gene Pairs Regulate Morphological Variations
Expression variations of miRNAs make it possible for responding to phytoplasmas invasion with reprogramming of gene expression.In this work, six DERs that may be related to plant defense or morphogenesis were selected to construct the miRNA-target gene interaction network based on gene function annotation.Previous studies demonstrated that downregulation of squamosa promoter binding protein-like 8 (SPL8) could increase branch development by facilitating the formation of axillary buds [53] and regulate the GA signaling pathway and biosynthesis in stem elongation [54].Gou et al. [53] indicated that GA signaling transduction was regulated by SPL8 from the upstream GA receptor (GID1) to the downstream responsive genes (GRASs).GID1 sensed and bound endogenous GA to induce the formation of GID1-GA-DELLA protein complex [53].There was increasing evidence that the complex caused rapid degradation of the DELLA protein [55] and DELLA as a nuclear transcription regulator could inhibit GA signaling and restrict plant growth [56].In addition, it was also reported that SPL regulate cell proliferation in leaf primordial [57].Leaf cell proliferation was positively controlled by functional transcriptional coactivators Ans and Grf-Interacting Factor 1 (GIF1) that affected the gene expression for leaf growth based on the interaction network [58,59].In this study, pf-miR2646 was upregulated in the phytoplasma-infected seedlings, which might decrease the expression level of its target gene (SPL8).Thereby, we hypothesize the downregulation of SPL8 inhibit the GA signal which may result in the stem section shortening in infected P. fortunei.Meanwhile, regulation of leaf cell proliferation by SPL may be involved in the occurrence of witches' broom.
Research has also shown that auxin responsive genes participated in biotic stresses to regulate the plant morphology.The involvement of ARFs in response to biotic stresses was also revealed in Arabidopsis thaliana [60], Gossypium [61] and Oryza sativa Linn [62].In addition, Fan et al. [63] suggested that the concentration of Indole-3-acetic acid (IAA) in PaWB-infected plants was significantly lower than that in healthy plants, while the concentration of cytokine (CK) in PaWB-infected plants was significantly higher than that in healthy plants.Moreover, the upregulation of CK/IAA could lead to the appearance of witches' broom symptoms.This is consistent with the results of that achieved in this study, i.e., ARFs (the target genes of pf-miR160h), are downregulated in the PFI library, implying that ARF TFs might play essential roles in response to PaWB-infection.

Regulation of Plant Defense Systems by DER-Target Genes Pairs
DERs modulate multiple biological processes not only in morphology, but also in the plant defense responsive to pathogen invasion.With the parallel evolution process of pathogen, the pattern of plant defense is constantly diversifying.Plants can establish a complete defense system to resist the invasion of pathogens by activating cascade response of protein kinases, strengthening the cell wall and building an salicylic acid (SA)-dependent immune response.As for strengthening the cell wall, the composition of the cell wall allows it to act as a major barrier to resist the invasion of pathogens [64].Zhong et al. [65] elucidated that over-expression of MYB can induce the gene expression changes of phenylpropanoid biosynthesis, which modulate the biosynthesis of lignin, one of the components of the cell wall [66,67].In Arabidopsis, the ectopic expression of MYB TFs AtPAP1 leads to increased gene expression of pheammonia-lyase, chalcone synthase and dihydroflavonol 4-reductase, which could increase the production of ligin and anthocyanin compounds [68].Evidence indicated that SECONDARY WALL-ASSOCIATED NAC DOMAIN PROTEIN1 (SND1) and its homologues such as NAC SECONDARY WALL THICKENING PROMOTING FACTOR1 (Nst1), Nst2, VASCULAR RELATED NAC-DOMAIN 7(VND7) and VND6 were the main switches activating a group of MYB TFs (including MYB20, MYB46, MYB85, MYB103), which in turn aroused the secondary wall synthesis process [65,69,70].In this study, pf-miR858a and pf-miR858b, whose target genes encoded MYB TFs (Arabidopsis MYB58 homolog), were downregulated, and the MYB genes displayed a higher expression level in PFI than that in PF plantlets.In consideration of the above-mentioned facts, we speculate that MYB and SND1 may work together on the synthesis of the cell wall to protect against pathogen invasion in P. fortunei.
As sessile organisms, plants have evolved unique cell surface receptors, which can sense chemical signals of intercellular communication to cope with various stresses.Receptor-Like Kinase (RLK) protein with a special structure (containing a signal sequence, an amino-terminal domain with a transmembrane region and a carboxyl-terminal kinase domain) is one of the dominant cell surface receptors that can receive signals of cell-to-cell communication.In Arabidopsis, leucine-rich repeatreceptor-like protein kinases (LRR)-RLKs containing an extracellular LRR and a Ser/Thr kinase domain were thought to be involved in plant defense processes via various signal cascade response [71,72].Recently, it has been reported that LRR domains of RLK could interact with multiple proteins resulting in signal transduction [73,74].MLO proteins are leucine-rich peptides and play key roles in the signal recognition process of plant defense [73,74].The binding of the endogenous plant elicitor peptides (PEPs) to their receptors (PEPRs) located on the plasma membrane resulted in an increase of cytosolic Ca 2+ , which was the early response of the signal cascade to activate the immune response [75,76].It also had been reported that the Ca 2+ increase played a critical role in MLO function [77].Peterhansel et al. [78] suggested that the development of MLO-mediated resistance required the involvement of two genes, ROR1 and ROR2, while MLO-mediated defense response might involve one or more small GTP-binding proteins of the ROP family [79].Therefore, Bhat's [80] research team believed that MLO, ROR2 and other proteins might develop a new pathogen-triggered microdomain at the site of biotic stress.Hence, among the interaction network when plant encounter stresses, the protein kinase family with leucine-rich repeat domain may interact with MLO to activate the signal cascade transduction to respond to phytoplasma invasion.In this study, the MLO6 (target of pf-miR164b-5p) is upregulated in PFI, indicating that the plant could trigger a series of defense responses when infected with phytoplasmas.
Salicylic acid (SA) plays an important role in protecting pathogen infection [81].As the key speed limiting enzyme in the biosynthesis of SA, Isochorismate synthase 1 (ICS1) was regulated by AtTCP8 [82].TCP8 showed positive regulation on the expression of ICS1 when interacting with the WRKY family [82].Moreover, in this interaction network, pf-miR2628 plays a key role in regulating the expression of WRKY4 (whose function was directly improved plant defense).In a word, during the process of PaWB-infection, pf-miR2628 and its targets may function as a key regulator to mediate the plant defense.
Taken together, our results revealed the significant roles of miRNA-target gene pairs in response to PaWB based on the interaction network, but further studies still need to be performed to explore their complex functions in Paulownia.

AS Events Related to PaWB Response in P. fortunei
AS is a vital regulatory mechanism at the post-transcriptional level in eukaryotes, which increases transcriptome complexity and protein diversity by splicing the same pre-mRNA [83].Research showed that AS can adjust the function of important plant stress-response components [84,85].In order to response to the invasion or attack from pathogens and pests, plant evolved defense mechanisms which involve multiple signal cascades and protein networks to provide specific and comprehensive defenses.
Serine/arginine-rich (SR) proteins belong to a conserved RNA-binding protein family, and play important roles in both constitutive splicing and AS.Recent studies have shown that alternative splicing of SR proteins play vital roles in stress response process [86][87][88].Xu et al. [86] have proved that MOS14 functioned as the nuclear import receptor of SR proteins, which play vital roles in RNA metabolism.MOS14 mutation changed the splicing patterns of resistance (R) gene SNC1 and RPS4, and then damaged the plant resistance, which was mediated by these two genes.As mentioned above, SR proteins can act as splicing factors for AS [89], and they also play a vital role in the plant immunity response.
In this study, the genes encoding SR proteins (PAU016740.2,PAU000297.1,and PAU019143.1)contained five AS categories (TSS, XSKIP, SKIP, AE, and XAE) in the PFI library.Taken together, these results suggest that AS events may influence the defense mechanisms of plants in response to pathogen invasion.We speculate that PaWB-infected plants could induce lots of SR proteins through AS events in order to respond pathogen invasion.In general, our results elucidate the crucial roles of AS events in answering to PaWB and provided new ideas to further investigate the transcriptome complexity

Figure 1 .
Figure 1.The PCA results of three biological replicates.(A) Represents three biological replicates of mRNA, while (B) represents miRNA.

Figure 3 .
Figure 3. Length distribution of small RNA reads.

Figure 5 .
Figure 5. Gene Ontology functional classification of identified target genes.

Figure 6 .
Figure 6.(A) Gene Ontology (GO) enrichment of the target genes; (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation for the target genes.

Figure 7 .
Figure 7.The interaction network of DERs and genes.Red diamonds represents DERs in PaWB-infected plants; the yellow circle represents the DEGs; the pink circle represents the gene that interacts with the target gene; the blue circle represents the downstream gene of the target gene; the green octagon represents the synergistic action with other genes; and the blue rectangle represents a specific function or tissue organ in plants.The thickness of gray lines between genes indicates the strength of the interaction; and as the black dotted and solid lines represent the indirect and direct effects, respectively.The arrows represent the positive regulation; and horizontal lines represent negative regulation.

Table 3 .
Mapped to the genome.

Table 4 .
Statistic of alternative splicing categories from mRNA data.