Comparative Analysis of MicroRNA Expression in Three Paulownia Species with Phytoplasma Infection

Paulownia witches’ broom (PaWB), caused by phytoplasma, is an important disease of Paulownia. To further identify the key miRNAs associated with the formation of PaWB symptoms, miRNA and degradome sequencing were performed to explore important miRNAs–target regulation in healthy and diseased Paulownia tomentosa, Paulownia fortunei, and P. tomentosa× P. fortunei seedlings, and the corresponding diseased seedlings treated with 75 mg L−1 dimethyl sulfate. A total of 212, 111, and 197 differentially expressed miRNAs (DEMs) were obtained in P. tomentosa, P. fortunei, and P. tomentosa × P. fortunei, respectively. Degradome sequencing detected 559, 251, and 568 target genes of the DEMs in P. tomentosa, P. fortunei, and P. tomentosa× P. fortunei, respectively. The expression patterns of selected miRNAs and the target genes were verified be qRT-PCR. Through analysis of the expression level of the DEMs in this study, combined with the results in our previous studies, as well as with those reported in other phytoplasma-infected plants, we concluded that miR156 is an important miRNA related to witches’ broom. According to the functions of the target genes of DEMs, we constructed a co-regulatory network of the DEMs-target genes interaction. These results will help to advance the understanding of the mechanism of PaWB.


Introduction
Paulownia are fast-growing trees native to China and have been cultivated for their environmental protection and commercial importance, but they are susceptible to Paulownia witches' broom (PaWB), which is one of the most devastating diseases caused by phytoplasma.The characteristics of infected Paulownia with witches' brooms include short internodes, yellowing or reddening of leaves, stunting and decline, phyllody and virescence of flowers, sterile flowers and necrosis, and altered volatile production [1,2].Phytoplasmas are obligate intracellular parasites with no cell walls that belong to the class Mollicutes.They have small genomes (about 530-1350 kb) with 23-29% GC content and can infect more than 1000 different plant species [3].Until now, no valid method has been developed to cultivate phytoplasma in vitro, which makes it difficult to investigate phytoplasma pathogenesis based on the pathogen.Therefore, researchers tend to focus on the molecular mechanisms in their host plants [4][5][6].Accumulating evidence has revealed the gene, protein, and metabolite changes in the host plants by differential methods [7][8][9], and miRNA sequencing is one of the methods used to research the plant-phytoplasma interaction [7].
MicroRNAs (miRNAs) are approximately 20-25 nt long and are important post-transcriptional regulators that modulate gene expression by targeting mRNAs for translational repression or cleavage [10,11].An increasing number of miRNAs have been identified and their sequences have been deposited in a publicly available database (http://www.mirbase.org/).The discovery of miRNAs as gene regulators has broadened the perspectives of the mechanisms involved in gene regulation, and revealed that miRNAs are associated with many biological processes such as plant development, protein degradation, cell proliferation and differentiation, and signal transduction, and play important roles in resistance to abiotic and biotic stresses [12,13].A number of studies of plant-phytoplasma interactions have shown that miRNAs are involved in gene regulation in many plants, including Mexican lime, Mulberry, Catharanthus roseus, Peanuts, and Ziziphus jujube [14][15][16].Some miRNA families are conserved among various plants; however, changes in their expression patterns in response to different phytoplasma infections were different in different plants.In Mexican lime infected with Candidatus Phytoplasma aurantifolia, a high expression of miR157 was shown to regulate changes in the leaf morphology of infected seedlings [15] In mulberry trees infected with Aster yellows phytoplasma, a high expression of miR393a was associated with changes in the morphology of infected mulberry seedlings, and upregulation of miR482a-5p was shown to be responsible for some symptoms in the infected plants [7,16].In Ziziphus jujube infected with Candidatus Phytoplasma ziziphi, overexpression of miR156b increased axillary branching, aberrantly expressed miR159 mediated MYB33, and the miR156-SPL-miR172 regulatory pathway was thought to lead to trolled leaf and shorter petioles and the development of green leaf-like structures instead of flowers, as well as the sterility of flowers [14].Recently, the differential expression of miR156 and miR172 was reported to lead to abnormal leaf shape and size, in addition to flower abortion in grapevine infected with Candidatus Phytoplasma asteris [17].
The functions of miRNAs in PaWB have also been investigated, and some miRNAs and their target genes involved in metabolism pathways have been reported [18][19][20][21][22][23], but most of the reported miRNAs were obtained from de novo transcriptome assemblies and were putative genes, because the lengths of unigene sequences obtained by transcriptome assembly are shorter than the reads obtained from the Paulownia reference genome sequence.This implies that the annotation information may be incomplete and a lot of time and effort may be required to validate their functions in the processes associated with PaWB occurrence.Further, individual miRNAs in the same family may be expressed differently in different paulownia tissues and have different functions in response to the same pathogen infection.Thus, phytoplasma-responsive miRNAs have not been fully explored and their roles in mediating gene expression regulation in response to phytoplasma are largely unknown.
In previous studies, we found that PaWB seedlings treated with 75 mg•L −1 dimethyl sulfate (DMS) recovered a healthy morphology and the phytoplasma 16SrRNA was not detected by Nested PCR [24][25][26].In this study, we explored key miRNA changes associated with the PaWB symptom forming process by small RNA sequencing, and used degradome sequencing to predict the target of miRNA.We compared gene expression changes in healthy and diseased P. tomentosa, P. fortunei, and P. tomentosa × P. fortunei, and in corresponding diseased seedlings treated with 75 mg•L −1 DMS at the genome level.The results will provide a novel platform to better understand paulowniaphytoplasma interactions.

Plant Materials
All plant materials used in this study were tissue cultured seedlings from the Laboratory of Forestry Biotechnology, Henan Agricultural University, Zhengzhou, Henan Province, China.The following nine samples were generated for use in this study: three healthy P. fortunei (PF), P. tomentosa (PT), and P. fortunei × P. tomentosa (PTF) seedling samples; three diseased PF, PT, and PTF Forests 2018, 9, 302 3 of 17 seedling samples named PFI, PTI, and PTFI, respectively; and three PFI, PTI, and PTFI seedling samples treated with 75 mg•L −1 DMS named PFI-75, PTI-75, and PTFI-75 respectively.The PF, PT, PTF, PFI, PTI, and PTFI seedlings were culture-grown in vitro on 1/2 Murashige-Skoog (MS) medium with 20 g•L −1 sucrose and 8 g•L −1 agar for 30 d. Subsequently, 1.5 cm terminal buds were collected from the PFI, PTI, and PTFI seedlings and soaked in 100-mL flasks with 0 mg•L −1 and 75 mg•L −1 DMS at 16 • C in a dark room, respectively.After soaking for 5 h, these seedlings were cultured on 1/2 MS medium with 20 g•L −1 sucrose and 8 g•L −1 agar in one 100-mL flask per sample.At the same time, the healthy PF, PT, and PTF seedlings were transferred to triangular flasks containing 1/2 MS culture medium without DMS as the control groups.For the control and treated seedlings, with three seedlings per flask, the total of each seedling was cultured in 20 100-mL flasks at 25 ± 2 • C under 130 µmol•m −2 •s −1 light intensity for 16 h•d −1 .After 30 day, the 1.5 cm buds were collected from the nine samples for further study.Each seedling was performed in triplicate.The excised samples were frozen immediately and stored at −80 • C.

Construction and Sequencing of Small RNA and Degradome Libraries
Total RNA was purified from the stem apexes of the nine samples (each sample was the mixture of three repeats), namely, healthy PF, PT, and PTF; diseased PFI, PTI, and PTFI; and PFI-75, PTI-75, and PTFI-75 seedlings using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocols.To construct the nine small RNA (sRNA) libraries, 18 to 30 nucleotides (nt) sRNAs were fractionated by polyacrylamide gel electrophoresis, and the fragments were then ligated with 5 and 3 RNA adapters by T4 RNA ligase.Reverse transcription and PCRs were performed to obtain sufficient single-stranded cDNA.The PCR procedure was conducted according to [23].An Agilent 21.00 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and ABI StepOnePlus Real-Time PCR System (Life Technologies, Grand Isle, NY, USA) were used to determine the quality and yield, respectively.Finally, the nine cDNA libraries were sequenced on an Illumina HiSeq™ 2000 platform.
To identify the potential target genes of the detected miRNAs, nine degradome libraries were also constructed for the nine samples.From the total RNAs, approximately 20 µg of total RNAs was used to prepare the degradome library.Briefly, poly (A) RNAs that possessed a 5 phosphate were extracted and ligated to an RNA adaptor containing a 3 MmeI recognition site using T4 RNA ligase.Then, the reverse transcription reaction and a short PCR were performed to obtain double-stranded DNA.Next, the PCR products were purified and digested with MmeI (NEB, Ipswich, MA, USA).A double-stranded DNA adapter was linked to the digested products and the ligated products were amplified by PCR and gel-purified to obtain the final degradome library.All the sRNA and degradome libraries were sequenced on an Illumina HiSeq™ 2000 system (Beijing Genomics Institute, Shenzhen, China).
After the conserved and novel miRNAs were identified in the nine libraries, miRNA expression was calculated as follows.The standard measure of reads per kilobase per million (RPKM) was used to normalize gene expression, where 0 was normalized as 0.01 to facilitate calculation.Normalized expression = (Actual sRNA sequencing read counts/Total clean read counts) × 1,000,000.Then, miRNA expression in different comparisons (PFI vs. PF, PTI vs. PT, PTFI vs. PTF, PFI-75 vs. PFI, PTI-75 vs. PTI, PTFI-75 vs. PTFI) was calculated, and miRNAs with |log 1.2 fold change| ≥ 0.263 and adjusted p-value < 0.05 were considered as significantly differentially expressed miRNAs (DEMs) [31,32].The fold-change was calculated according to the following equation: log1.2 fold changes = log 1.2 (normalized miRNA reads in treatment/normalized miRNA reads in control).
The p-value was obtained according to the calculations as follows: where N 1 is the total number of reads in the control library, N 2 is the total number of reads in the DMS treatment library, x is the number of reads for an miRNA in the control library, and y is the number of reads for an miRNA in the DMS treatment library.

Identification of miRNAs Related to PaWB
The abundance of miRNAs in the nine libraries was analyzed to obtain the DEMs related to PaWB.First, the DEMs in 1 PFI vs. PF, 2 PTI vs. PT, and 3 PTFI vs. PTF were compared.The same DEMs from these three comparisons were obtained 4 .Then, the DEMs in 5 PFI-75 vs. PFI, 6 PTI-75 vs. PTI, and 7 PTFI-75 vs. PTFI were compared.Similarly, the same DEMs in these three comparisons were generated 8 .Finally, 9 PaWB related to miRNAs were identified from the DEMs that intersected in 4 and 8 (Figure 1). the RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi)[27].Candidate miRNA sequences with a predicted secondary structure minimal fold energy index (MFEI) ≥0.9 were considered to be novel miRNAs [28][29][30].
After the conserved and novel miRNAs were identified in the nine libraries, miRNA expression was calculated as follows.The standard measure of reads per kilobase per million (RPKM) was used to normalize gene expression, where 0 was normalized as 0.01 to facilitate calculation.Normalized expression = (Actual sRNA sequencing read counts/Total clean read counts) × 1,000,000.Then, miRNA expression in different comparisons (PFI vs. PF, PTI vs. PT, PTFI vs. PTF, PFI-75 vs. PFI, PTI-75 vs. PTI, PTFI-75 vs. PTFI) was calculated, and miRNAs with |log 1.2 fold change| ≥ 0.263 and adjusted p-value < 0.05 were considered as significantly differentially expressed miRNAs (DEMs) [31,32].The fold-change was calculated according to the following equation: log1.2 fold changes = log 1.2 (normalized miRNA reads in treatment/normalized miRNA reads in control).
The p-value was obtained according to the calculations as follows: where N1 is the total number of reads in the control library, N2 is the total number of reads in the DMS treatment library, x is the number of reads for an miRNA in the control library, and y is the number of reads for an miRNA in the DMS treatment library.

Identification of miRNAs Related to PaWB
The abundance of miRNAs in the nine libraries was analyzed to obtain the DEMs related to PaWB.First, the DEMs in ① PFI vs. PF, ② PTI vs. PT, and ③ PTFI vs. PTF were compared.The same DEMs from these three comparisons were obtained ④.Then, the DEMs in ⑤ PFI-75 vs. PFI, ⑥ PTI-75 vs. PTI, and ⑦ PTFI-75 vs. PTFI were compared.Similarly, the same DEMs in these three comparisons were generated ⑧.Finally, ⑨ PaWB related to miRNAs were identified from the DEMs that intersected in ④ and ⑧ (Figure 1).

Analysis of Target Genes of the Differentially Expressed miRNA
To identify miRNA target genes, we constructed nine degradome libraries that corresponded to the nine sRNA libraries described above.The raw sequence reads were obtained using Illumina HiSeq™ 2000.After initial processing, the unique sequence signatures were mapped to the Paulownia reference genome (http://soap.genomics.org.cn/) using SOAP software (http://soap.genomics.org.cn/) to define the coverage rate.Perfectly matched reads were retained and extended to 31-nt long signatures by adding approximately 15 nt of the upstream sequence from the aligned Paulownia sequences.All the resulting reads were reverse complemented and aligned

Analysis of Target Genes of the Differentially Expressed miRNA
To identify miRNA target genes, we constructed nine degradome libraries that corresponded to the nine sRNA libraries described above.The raw sequence reads were obtained using Illumina HiSeq™ 2000.After initial processing, the unique sequence signatures were mapped to the Paulownia reference genome (http://soap.genomics.org.cn/) using SOAP software (http://soap.genomics.org.cn/) to define the coverage rate.Perfectly matched reads were retained and extended to 31-nt long signatures by adding approximately 15 nt of the upstream sequence from the aligned Paulownia sequences.
All the resulting reads were reverse complemented and aligned back to the corresponding identified miRNA sequence in the nine libraries.Alignments with scores <4 and no mismatches between the 10th and 11th nucleotides of the complementarity miRNA sequences were considered to be potential targets.T-plots were built based on the distribution of signatures along these transcripts.To better understand the functions of the predicted target genes, the best homologs were identified and used to assign Gene Ontology (GO) annotations (p-value < 0.05 was used as the threshold).

Verification of DEMs and Their Targets by qRT-PCR Analysis
To validate the results of the high-throughput sequencing, six DEM-target gene pairs were selected randomly for qRT-PCR validation.First-strand cDNA was synthesized for all the samples using an iScript cDNA synthesis kit (Bio-Rad, Hercules, CA, USA).The primers for the selected target genes were designed using the Beacon Designer software (version 7.7, Premier Biosoft International, Ltd., Palo Alto, CA, USA), and the specific stem-loop primers for the miRNAs were designed as reported previously [33].The PCR amplifications conditions were according to [23].The U6 snRNA and 18S rRNA served as the miRNA and target gene endogenous references, respectively.Relative changes in expression of the miRNAs and their target genes were calculated using the 2 −∆∆Ct method [34].Each gene was analyzed in three replicates.The primers sequences for the miRNAs and their target genes are listed in Table S1.back to the corresponding identified miRNA sequence in the nine libraries.Alignments with scores <4 and no mismatches between the 10th and 11th nucleotides of the complementarity miRNA sequences were considered to be potential targets.T-plots were built based on the distribution of signatures along these transcripts.To better understand the functions of the predicted target genes, the best homologs were identified and used to assign Gene Ontology (GO) annotations (p-value < 0.05 was used as the threshold).

Verification of DEMs and Their Targets by qRT-PCR Analysis
To validate the results of the high-throughput sequencing, six DEM-target gene pairs were selected randomly for qRT-PCR validation.First-strand cDNA was synthesized for all the samples using an iScript cDNA synthesis kit (Bio-Rad, Hercules, CA, USA).The primers for the selected target genes were designed using the Beacon Designer software (version 7.7, Premier Biosoft International, Ltd., Palo Alto, CA, USA), and the specific stem-loop primers for the miRNAs were designed as reported previously [33].The PCR amplifications conditions were according to [23].The U6 snRNA and 18S rRNA served as the miRNA and target gene endogenous references, respectively.Relative changes in expression of the miRNAs and their target genes were calculated using the 2 −ΔΔCt method [34].Each gene was analyzed in three replicates.The primers sequences for the miRNAs and their target genes are listed in Table S1.A total of 430 (PF), 473 (PT), and 477 (PTF) miRNAs were identified, and 267 (PF), 323 (PT), and 319 (PTF) of them were conserved miRNAs.These miRNAs are listed in (Table S3).Most of the conserved miRNAs in the three paulownia species were 21 nt long.The average length of the precursor sequences of the conserved miRNAs was approximately 130 nt and ranged from 51 to 255 A total of 430 (PF), 473 (PT), and 477 (PTF) miRNAs were identified, and 267 (PF), 323 (PT), and 319 (PTF) of them were conserved miRNAs.These miRNAs are listed in (Table S3).Most of the conserved miRNAs in the three paulownia species were 21 nt long.The average length of the precursor sequences of the conserved miRNAs was approximately 130 nt and ranged from 51 to 255 nt.According to the criterion used to predict novel miRNAs, a total of 163 (PF), 150 (PT), and 158 (PTF) sequences were predicted as potential novel miRNAs (Table S4).All the novel miRNAs had a MFEI > 0.85; the average MFEIs were 1.24 (PF), 1.22 (PT), and 1.20 (PTF).The novel miRNAs varied in length from 18 to 25 nt with a peak at 24 nt, which is similar to the length distribution of the conserved miRNAs.The lengths of the novel miRNA precursor sequences varied from 63 to nt, and the free energy varied from −15 kcal mol −1 to −140.6 kcal mol −1 .The stem-loop sequences of the novel miRNA precursors varied from 54 to 220 nt, with an average of 129 nt.The Pearson correlations between samples are shown in (Figure S1), which demonstrates the repeatability of the miRNA sequencing.

Differentially Expressed miRNAs in the Three Paulownia Species
To identify the miRNAs associated with PaWB, the expression levels of the known and novel miRNAs were calculated between the different libraries (PFI vs. PF, PTI vs. PT, PTFI vs. PTF, PFI-75 vs. PFI, PTI-75 vs. PTI, and PTFI-75 vs. PTFI).A total of 245, 260, and 224 significant DEMs were detected in the PF vs. PFI, PT vs. PTI, and PTF vs. PTFI comparisons, respectively, and a total of 197, 397, and 406 DEMs were detected in the PFI vs. PFI-75, PTI vs. PTI-75, and PTFI vs. PTFI-75 comparisons, respectively.Finally, we identified 111 (PF), 212 (PT), and 197 (PTF) DEMs (Table S5), among which 78 (PF), 166 (PT), and 163 (PTF) were conserved miRNAs and 33 (PF), 46 (PT), and 35 (PTF) were novel miRNAs.According to the comparison scheme of the miRNA associated with the PaWB, we obtained 26 common DEMs from the three paulownia species (Table 1).These results demonstrate that the number of DEMs varied in the different comparisons, and the miRNA families and individual miRNAs were not same in the three paulownia species, implying that different miRNA changes occurred in the three paulownia species during the phytoplasma infection process.S6), and 251 (PF), 559 (PT), and 568 (PTF) target genes were identified (Table S7).Finally, 196 genes were identified as targets of 76 common DEMs (35 common miRNA families) in the three paulownia species (Table S8).No target genes were detected for most of the DEMs.A possible reason for this could be that the target mRNAs were expressed at low levels, so the cut ends were not detected by degradome sequencing.Another possible reason is that some miRNAs may function by inhibiting translation rather than by mRNA digestion.

Functional Analysis of the Target Genes of the DEMs
In this study, we identified 76 common DEMs in the three Paulownia species by miRNA sequencing.Subsequently, a total of 196 target genes were predicted by degradome sequencing.To understand the functions of the target genes of the DEMs involved, we mapped the target genes to the GO database (Table S8).The results showed that the target genes were mainly involved in the nucleus (GO: 0005634), DNA binding (GO: 0003677), metal ion binding (GO: 0046872), regulation of transcription, DNA-dependent (GO: 0006355), and cell differentiation (GO: 0030154).Subsequently, according to the functional descriptions of the target genes in other transplant plants, we grouped the 196 target genes of the 76 DEMs into five major pathway categories: genes associated with photosynthesis; genes involved in plant hormones; genes involved in plant defense and related signal pathways; genes involved in the energy metabolism pathway; and genes involved in material metabolic pathways.

Genes Associated with Photosynthesis
In this section, there were three miRNAs (miR156c-5p, miR160a-5p, and miR156y) that medicated their target genes involved in photosynthesis.Two target genes of miR156c-5p encode phytochrome-interacting factor 3 (PIF3), magnesium chelatase subunit I (CHLI), and the target of miR156y encode monothiol glutaredoxin (Grxs); the target genes of miR160a-5p encode the RAV-like factor (RAV).In this study, miR156c-5p and miR156y were up-regulated in the phytoplasma-infected Paulownia seedlings and the two target genes were down-regulated.PIF3 is a key basic helix-loop-helix transcription factor, and Liu [35] reported that PIF3 associates with histone deacetylase HDA15 affected chlorophyll biosynthesis.CHLI is an important enzyme for chlorophyll a and chlorophyll b synthesis [36], and the down-regulation of CHLI might reduce the PaWB seedling photosynthesis.Vranish [37] confirmed that Grxs as the intermediate carrier plays a positive role in the Fe-S cluster biosynthetic pathway, while the down-regulated Grxs would hinder the Fe-S cluster biosynthetic pathway.MiR160a-5p was down-regulated in PaWB seedlings.The target of this miRNA encodes RAV that was up-regulated, and Zhao [38] showed overexpressing RAV reduced the chlorophyll content and photosynthetic rate in transgenic tobacco leaves.Therefore, this group of miRNAs and their target genes may associate with the reduction of photosynthesis.

Genes Associated with Plant Hormones
In this section, there were five DEMs (miR156c-5p, miR160a-5p, miR172a-3p, miR172d, and miR5368g-5p) that medicated their target genes regulation involved in plant hormones metabolism.The target genes of miR156c-5p encode PIF3; the targets of miR160a-5p encode RAV; the target gene of miR172a-3p encodes AP2-like factor (AP2) and the target gene of miR172d encode EREBP-like factor (EREBP); and the target gene of miR5368g-5p encodes zeaxanthin epoxidase (ZEP).In this study, miR156c-5p was up-regulated in the phytoplasma-infected Paulownia seedlings and the PIF3 target gene was down-regulated.The differently expressed PIFs are involved in brassinosteroids biosynthesis, and modulated plant growth and development [39].In PaWB seedlings, miR160a-5p, miR172a-3p, and miR172d were down-regulated and their targets RAV, AP2, and EREBP were up-regulated.AP2 and EREBP are members of the AP2/ERF superfamily, which plays vital roles in regulating various developmental and physiological processes in plants [40].In Arabidopsis, overabundance of the AP2/ERF transcription factor (WIND1) increased the B-type ARR-mediated cytokinin response, which also promoted cell dedifferentiation [41].RAV is also a member of the AP2/ERF transcription factor family that has been shown to be increased in CK transcript abundance.In transgenic tobacco, overexpressed GmRAV increased cytokinin signaling-related phenotypes including dwarfism and reduced apical dominance [42].miR5368g-5p was down-regulated in the PTI and PTFI, according to the negative regulation of miRNA and its target gene, the target gene ZEP was up-regulated, which functions in ABA biosynthesis [43].Therefore, this group of miRNAs and their target genes may determine some of the symptom of PaWB.

Genes Involved in Plant Defense and Related Signal Pathways
In this study, four main DEMs (miR172a-3p, miR172d, miR4414a-5p, and miR156y) medicated their target genes involved in plant defense and related signal pathways.The target gene of miR172 a-3p encode AP2 and the target gene of miR172d encode EREBP and glutathione peroxidase (GPX).In Arabidopsis thaliana, an AP2/ERF transcription factor (MACD1), was found to be involved in phytotoxin-triggered programmed cell death [44].Huang, et al. [45] reported that AP2/ERF transcription factors were reported to interact with the MAPK, and jasmonic acid, salicylic acid, ethene, and H 2 O 2 pathways to improve the ability of tomato plants to respond to yellow leaf curly virus.GPX is a redox sensor protein that maintains a steady-state of H 2 O 2 in plant cells and has a pivotal role in peroxide detoxification in eukaryotic organisms [46].For the DEM, like mi4414a-5p, its expression level is up-regulated in PTI and PTFI.One of the predicted targets that encodes the cyclic nucleotide gated channel (CNGC) was down-regulated, which can provide a pathway for Ca 2+ conductance across plasma membranes, and Ma and Berkowitz [47] reported that CNGC may participate in pathogen defense signal transduction cascades.One of the gene targets of miR156y encodes phospholipase D (PLD), which combined with PLD-derived phosphatidic acids, is an important signaling complex that plays vital roles in plant defense in plant-microbe interactions [48].Although the targets like CNGC and PLD were down-regulated in the PaWB seedlings, these signal cascade proteins may combine other related signal transduction pathways to regulate plant defense.

Genes Involved in Energy Metabolism Pathway
Energy metabolism is one of the unsolved mysteries in the plant-phytoplasma interaction.In this study, the target genes of miR477 b encode V-type H + -transporting ATPase subunit A(V-ATPase), and the target genes of miR5368f-3p encode F-type H + -transporting ATPase subunit gamma(F-ATPase) and the target genes of miR5368g-5p encode NADH dehydrogenase (ubiquinone) 1 alpha/beta subcomplex 1 (NADUFA/B1) involved in the energy metabolism.Fillingame [49] reported that both F-type and V-type ATP synthases can catalyze the synthesis of ATP in oxidative phosphorylation.NADH dehydrogenase (ubiquinone) 1 alpha/beta subcomplex 1 is an important enzyme for energy production in oxidative phosphorylation.In this study, miR5368g-5p and miR5368f-3p were down-regulated in the PTI and PTFI, according to the negative regulation of miRNA and its target gene, F-ATPase and NADUFA/B1were up-regulated, and the differential expressed miR477b medicated V-ATPase also changed the energy production in the oxidative phosphorylation pathway.It has been reported that phytoplasma lack the oxidative phosphorylation pathway [50].The energy source of the PaWB phytoplasma is still unknown.Our results assume that phytoplasma infection may affect energy metabolism in host paulownia in oxidative phosphorylation.However, how these miRNAs Forests 2018, 9, 302 9 of 17 mediate their target genes to regulate the energy metabolism in paulownia-phytoplasma interactions needs further elaboration.

Genes Involved in Material Metabolic Pathways
In this study, many target genes of DEMs took part in the material metabolic pathways.The target gene of miR5368h-5p encodes RNA methyltransferase (RNMT), and miR5368s-5p encodes RNA methyltransferase pyruvate kinase (PK), the two target genes of miR394a-5p encode L-arabinokinase (ARA) and glycogen synthase kinase 3 beta (GSK3B), the target genes of miR3630 encodes phosphoglycerate mutase (PGM).Besides, many target genes encoded transcription factors involved in transcription; for example, the target gene of miR477b encodes the transcription initiation factor TFIIF alpha subunit (TFIIF), the target gene of miR319, miR858a, and miR396a-3p encodes an myb proto-oncogene protein (MYB), and the target gene of miR167d-5p encodes the MADS-box transcription factor (MADS), which all play various roles in plant development.
In this study, according to the comparison scheme of the miRNAs related to PaWB in the material section, a total of 76 common DEMs related to PaWB were generated, and 196 target genes were predicted by degradome sequencing.Of which, most DEMs targeted more than one target gene, and the same target gene may be cleaved by different DEMs.Through an analysis of the functions of the target genes in other transplant plants, we found that each target gene was involved in several physiological processes.These results, together with the results of previous studies, indicate that an integrated co-regulatory network exists in the pathological process of phytoplasma infection [23].
Based on an integrated analysis of the functions of the target genes of the DEMs in this study and the previous ones [23,[51][52][53], we constructed a network that describes the regulation relationship of the DEM-mediated target genes in the manifestation of PaWB symptoms (Figure 3).

Expression Patterns of the DEMs and Their Target Genes by qRT-PCR
We selected six miRNAs from nine libraries and the corresponding six target genes for qRT-PCR assays.The results showed that DEMs and their target genes' expression (pau-miR156c-5p-PIP3, pau-miR319b-TCP2, pau-miR2118a-disease resistance protein RPP8, pau-miR156d-CHLI, pau-miR160c-5p-RAV) had a negative regulation pattern (Figures 4 and 5), and there was a strong correlation between the read abundances in the sequencing data and the expression levels obtained by qRT-PCR (Figure S2 and Figure 4).These results indicated that the miRNA expression profiles estimated from the Illumina sequencing data were quantitative and reliable.
expression levels obtained by qRT-PCR (Figures S2 and 4).These results indicated that the miRNA expression profiles estimated from the Illumina sequencing data were quantitative and reliable.In this picture, PIP3 was targeted by pau-miR156c-5p in Figure 4, transcription factor TCP2 was targeted by pau-miR319b-3p in Figure 4, disease resistance protein RPP8 was targeted bypau-miR2118a in Figure 4, CHLI was targeted by pau-miR156d in Figure 4, RAV was targeted by pau-miR160c-5p in Figure 4, and PE was targeted by pau-miR397a-5p in Figure 4. expression levels obtained by qRT-PCR (Figures S2 and 4).These results indicated that the miRNA expression profiles estimated from the Illumina sequencing data were quantitative and reliable.In this picture, PIP3 was targeted by pau-miR156c-5p in Figure 4, transcription factor TCP2 was targeted by pau-miR319b-3p in Figure 4, disease resistance protein RPP8 was targeted bypau-miR2118a in Figure 4, CHLI was targeted by pau-miR156d in Figure 4, RAV was targeted by pau-miR160c-5p in Figure 4, and PE was targeted by pau-miR397a-5p in Figure 4.In this picture, PIP3 was targeted by pau-miR156c-5p in Figure 4, transcription factor TCP2 was targeted by pau-miR319b-3p in Figure 4, disease resistance protein RPP8 was targeted bypau-miR2118a in Figure 4, CHLI was targeted by pau-miR156d in Figure 4, RAV was targeted by pau-miR160c-5p in Figure 4, and PE was targeted by pau-miR397a-5p in Figure 4.

MiR156 Is the Key Regulatory Factor Related to Witches' Broom
Many studies have shown that sRNAs play important roles in plant growth and development, defense against viruses and transposons, chromatin modifications, and responses to biotic and abiotic stresses.Although some progress has been made in understanding paulownia-phytoplasma interactions, various aspects still need to be verified [23,[51][52][53].Firstly, the accuracy of annotation information needs to be improved.Annotations based on transcriptome level data can easily generate putative or unclear annotation results.In the present study, all the annotation information for the three Paulownia species was based on the P. fortunei reference genome sequence, which can increase the accuracy of the annotations.Secondly, to verify the effectiveness of DEMs, we used the results of our previous researches [20,21,[24][25][26]54], which showed that the three PaWB seedlings could revert to healthy morphology after treatment with methyl methanesulfonate and DMS.The DEMs identified in previous researches were obtained based on morphological changes in seedlings treated with methyl methanesulfonate, whereas, in this study, miRNA changes were obtained in the three DMS-treated paulownia species.When the same miRNAs were detected in seedlings treated with both reagents, the miRNAs were identified as key miRNAs related to PaWB.Thence, we compared the common DEMs in the three paulownia species with the DEMs detected in our previous studies [23,[51][52][53][54], and found that miR156, miR398, miR408, and miR2118 were common DEMs in the healthy vs. phytoplasma-infected seedling comparisons.Thirdly, to confirm the key miRNAs related to witches' broom, we compared the four identified miRNAs with the miRNAs detected in phytoplasma-infected Mexican lime, Ziziphus jujuba, and mulberry [7,[14][15][16], and found that miR156 was the only common miRNA among the phytoplasma-responsive miRNAs in all these studies.We then analyzed the expression of miR156 in the three paulownia species and compared the expression level of miR156 of other phytoplasma-infected plants, and found their expression levels in most of the above plants were up-regulated in the phytoplasma-infected plants.Thus, we consider that miR156 may be an important miRNA that induces the symptom of witches' broom.

A miR156-miR160-miR5368 Cluster Regulates Leaf Yellowing in Phytoplasma-Infected Paulownia
To reveal the important factors for PaWB occurrence, in our laboratory, we have investigated the transcriptome, DNA methylation, metabolome, and proteome of the three Paulownia species [18,19,[55][56][57].Most of these researches have indicated that photosynthesis is arrested in PaWBinfected seedlings, which is similar to what was reported in other phytoplasma-infected plants [58,59].In our previous transcriptome research, we identified genes that were differentially expressed in photosystem I and photosystem II, and were associated with decreased photosynthesis [18,19,55].In the previous proteome study, we not only detected differentially abundant proteins associated with the above two photosynthesis units, but also found that the amounts of differentially abundant proteins changed in electron transfer and chloroplast development in response to phytoplasma infection [56,57].Although previous researches also investigated changes in miRNA expression in phytoplasma infected PT, PF, and PTF seedlings, the relationship between photosynthesis and changes of miRNAs has not been reported.In the present study, we identified several miRNAs and their targets that have not been previously described in paulownia, and that also took part in the photosynthetic pathway.
In this study, we detected three DEMs, miR156c-5p, miR160a-5p, and miR5368g-5p, that were related to chloroplast development in Paulownia leaves.MiR156c-5p was induced in the three healthy PT, PF, and PTF seedlings and the three recovered PTI-75, PFI-75, and PTFI-75 seedlings.This miRNA was predicted by degradome sequencing to target 94 genes, but only 17 of them were annotated at the genome level; among them, two encoded PIF3 and CHLI.PIF3 is a basic helix-loop-helix (bHLH) transcription factor that interacts with phyA and phyB and binds specifically to a cis-acting regulatory element (G-box) in the promoters of a variety of phytochrome-responsive genes, and evidence indicated that PIF3 was also a positive regulator of PHYB-mediated signal transduction in plants overexpressing PIF3 and in antisense PIF3 plants [60].It has also been reported that PIF3 mutants not only negatively induce chlorophylls and the assembly of photosynthetic complexes, but also act as negative regulators that affect carotenoid accumulation [61,62].In this study, PIF3 was down-regulated in the phytoplasma-infected PTI, PFI, and PTFI seedlings.Thus, the down-regulation of PIF3 might restrain the production of chlorophylls and carotenoids, which might act as a coordinately regulated response for the deetiolation of seedlings after phytoplasma infection.This result is in accordance with a previous report based on transcriptome data [19].Another predicted target of miR156c-5p encodes CHLI, which is an important enzyme for chlorophyll a and chlorophyll b synthesis [36].In transgenic tobacco, CHLI deficiency resulted in lower chlorophyll content and reduced antenna complex [63].In this study, CHLI was down-regulated in the three phytoplasma-infected Paulownia species, which indicated that chlorophyll biosynthesis may be reduced.The reduction of chlorophyll biosynthesis explains the yellowing of Paulownia leaves after phytoplasma infection.Another DEM, miR160a-5p, was down-regulated in PaWB seedlings.The target of this miRNA encodes RAV that was up-regulated, which contains AP2/ERF and B3 domains.The RAV encodes transcriptional regulators with a variety of functions related to developmental and physiological processes in plants [64,65].Transgenic tobacco overexpressing RAV showed reduced chlorophyll content and a reduced photosynthetic rate in its leaves [38].Another DEM, MiR5368g-5p, was down-regulated in PTI and PTFI seedlings.miR5368g-5p was predicted to target ZEP, which is the key enzyme responsible for the accumulation of antheraxanthin and violaxanthin [66].Yang, et al. [67] reported that a high expression of ZEP increased the stored carotenoid pigments in transgenic tomato plants.In this study, ZEP was up-regulated in the three PaWB seedlings, implying that carotenoids were abundantly produced in phytoplasma-infected Paulownia leaves.Hence, we concluded that a miR156c-5p-miR160a-5p-miR5368g-5p cluster mediated their target genes to increase carotenoids content and reduce chlorophyll content and antenna complexes, which can explain the yellowing of leaves in the phytoplasma-infected paulownia seedlings.

A miR160-miR172-miR397 Cluster Regulates Leaf Dwarf Morphology in Phytoplasma-Infected Paulownia
Changes in plant architecture, such as dwarfism, have been observed in many plants and are mostly caused by hormonal imbalances [58,68].MiRNAs play important roles in hormone medicated morphology development [69].In our previous proteome and miRNA sequencing studies of paulownia-phytoplasma interactions, we detected changes in gibberellins and brassinosteroids, and abscisic acid biosynthesis in response to phytoplasma infection [23,56].A previous transcriptome study identified two key cytokinin biosynthesis enzymes that were obviously induced [19].In this study, we identified four DEMs (miR160a-5p, miR172a-3p, miR172d, and miR397a-5p) and their target genes that regulated the dwarf morphology in phytoplasma-infected Paulownia leaves by regulating the cytokinin pathway.In PaWB seedlings, miR160a-5p was down-regulated and the target RAV was up-regulated.In transgenic tobacco, overexpression of a RAV-like orthologue increased cytokinin signaling, which reduced the apical dominance and produced seedlings with dwarfism morphology, whereas inhibition of the RAV-like orthologue produced the opposite phenotype in the seedlings [38,42].In this study, miR172a-3p and miR172d were down-regulated in the PaWB seedlings, the target gene of miR172a-3p encoded AP2, the target gene of miR172d encoded EREBP, the target genes were up-regulated, and AP2 and EREBP belonged to an AP2/ERF superfamily member.In Larix kaempferi, overexpression of the AP2/ERF transcription factor increased the number of branches and produced shorter and smaller rosette leaves [70].Another DEM, like miR397a-5p, was up-regulated in PTI and PTFI seedlings, but the target gene that encodes pectinesterase (EC.3.1.1.11)was also up-regulated in the three PaWB seedlings.Pectinesterase is the first enzyme that acts on pectin.Hasunuma, et al. [71] reported that the overexpression of pectinesterase in transgenic tobacco led to short internodes, small leaves, and dwarf symptoms.In this study, miR397a-5p-pectinesterase displayed imcomplete negative regulation, which demonstrated the complex regulation of miRNA-medicated target genes.

Conclusions
In this study, miRNA and degradome sequencing were performed to explore important miRNAs related to PaWB in nine miRNA and nine degradome libraries.The miRNA sequencing resulted in the identification of 430 miRNAs (267 known; 163 novel) (PF), 473 miRNAs (323 known; 150 novel) (PT), and 477 miRNAs (319 known; 158 novel) (PTF).Among them, 111 (PF), 212 (PT), and 197 (PTF) miRNAs were differentially expressed.The degradome sequencing detected 251 (PF), 559 (PT), and 568 (PTF) target genes for the DEMs.By comparing the DEMs detected in this study and our previous studies, as well as with the phytoplasma-responsive miRNAs reported previously, we identified miR156 as the key miRNA related to witches' broom.Based on the functions of the target genes of the DEMs, we constructed a network of the miRNA pathway involved in PaWB, and discovered a pau-miR156c-5p-miR160a-5p-miR5368g-5p cluster that regulated yellowing of leaves and a pau-miR160a-5p-miR172a-3p-miR172d-miR397a-5p cluster that regulated the dwarf morphology of phytoplasma-infected Paulownia.These results suggest that a co-regulatory network may be involved in the post-transcriptional regulation of genes associated with the pathological development of PaWB.In future work, we will focus on confirming the functions of these miRNAs by transgenic experiments, which will provide more information to understand the mechanism of PaWB occurrence.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/9/6/302/s1,Table S1: Primers of DEGs and target genes for qRT, Table S2: Read abundance of various classes of sRNAs sequences in nine libraries, Table S3: The conserved miRNA in the three Paulownia, Table S4: The new miRNA in the three Paulownia, Table S5: Differentially expressed miRNAs in the three paulownia species, Table S6: Data summary of degradome library, Table S7: Degradome sequencing of the three Paulownia species, Table S8

Figure 1 .
Figure 1.The comparative scheme of the miRNA related to PaWB.

Figure 1 .
Figure 1.The comparative scheme of the miRNA related to PaWB.

Figure 2 .
Figure 2. The length distribution of miRNAs generated in nine samples.

Figure 2 .
Figure 2. The length distribution of miRNAs generated in nine samples.

Figure 3 .
Figure 3.The potential regulatory network of the DEMs-target genes.

Figure 4 .
Figure 4. qPCR analysis of the expression level of the miRNAs.

Figure 5 .
Figure 5. qPCR analysis of the expression level of the target genes of the selected miRNAs.

Figure 4 .
Figure 4. qPCR analysis of the expression level of the miRNAs.

Figure 4 .
Figure 4. qPCR analysis of the expression level of the miRNAs.

Figure 5 .
Figure 5. qPCR analysis of the expression level of the target genes of the selected miRNAs.

Figure 5 .
Figure 5. qPCR analysis of the expression level of the target genes of the selected miRNAs.
: Differentially expressed miRNAs targets identified by degradome sequencing, Figure S1: Pearson correlation between samples, Figure S2: The results of miRNA sequencing of the selected miRNAs.Author Contributions: G.F. conceived, designed, and performed the experiments.Y.Z.analyzed the miRNA data.Z.C. and X.L. performed sample treatment and seedling inoculation.X.C. and X.Z.wrote the paper.All authors read the script, provided comments, and approved the final manuscript.Funding: This research was funded by Natural Science Foundation of Henan Province of China (162300410158).