Long Non-Coding RNAs Responsive to Witches ’ Broom Disease in Paulownia tomentosa

Paulownia witches’ broom (PaWB) disease caused by phytoplasmas is a fatal disease that leads to considerable economic losses. Long non-coding RNAs (lncRNAs) have been demonstrated to play critical regulatory roles in posttranscriptional and transcriptional regulation. However, lncRNAs and their functional roles remain poorly characterized in Paulownia. To identify lncRNAs and investigate their roles in the response to PaWB phytoplasmas, RNA sequencing was performed for healthy Paulownia tomentosa, PaWB-infected P. tomentosa, and for healthy and PaWB-infected P. tomentosa treated with 100 mg L−1 rifampicin. A total of 28,614 unique mRNAs and 3693 potential lncRNAs were identified. Comparisons between lncRNAs and coding genes indicated that lncRNAs tended to have shorter transcripts and fewer exon numbers, and displayed significant expression specificity. Based on our comparison scheme, 1063 PaWB-related mRNAs and 110 PaWB-related lncRNAs were identified; among them, 12 PaWB-related candidate target genes that were regulated by nine PaWB-related lncRNAs were characterized. This study provides the first catalog of lncRNAs expressed in Paulownia and gives a revealing insight into the molecular mechanism responsible for PaWB.


Introduction
Phytoplasmas are minute bacteria without cell walls and are the causal agents of witches' broom disease in Paulownia.Witches' broom disease occurs in many countries and is the biggest threat to Paulownia production, causing serious economic losses [1,2].In Paulownia, witches' broom disease has been found in saplings and big trees, and leads to malformations of branch, leaf, stem, and flower, such as witches' broom, yellowing, phloem necrosis, and phyllody, respectively [3,4].
Paulownia trees are native to China and have been introduced into many countries [5].The superior traits of Paulownia, such as distinguished adaptive capacity to poor habitat, extremely fast growth, and high-quality timber, make it a popular tree species with important economic and ecological values [6].Accordingly, Paulownia witches' broom (PaWB) disease has been studied by biologists for many years.Vitamin C plays an important role in the resistance to PaWB [7], and plant hormones are known to be related to the morphogenesis of PaWB [8].The Sec protein translocation system [9], elongation factor EF-Tu [10], and two plasmids [11] of the PaWB phytoplasma have been analyzed.Transgenic technology has been used in the breeding of PaWB-resistant Paulownia [12].Rifampicin treatment can make PaWB-infected Paulownia recover to normal morphology [13].Furthermore, the rapid development of 'omics' has allowed for the investigation Forests 2017, 8, 348 2 of 15 of the transcriptomes [3,[14][15][16][17][18], microRNAs (miRNAs), degradomes [19][20][21][22], and proteomes [4,23] of Paulownia species, and the results have revealed changes after PaWB infection at transcriptional, post-transcriptional, and translational levels.
Long non-coding RNAs (lncRNAs) are generally defined as non-coding transcripts that are more than 200 nucleotides in length [24].Large numbers of lncRNAs have been identified in Arabidopsis thaliana [25][26][27], rice [28,29], wheat [30,31], maize [32], cotton [33,34], Brassica napus [35], and Populus [36,37], and have been found to play important roles in growth, development, differentiation, and stress responses.LncRNAs involved in the responses of plants to biotic stresses have also been identified.For example, in wheat, four lncRNAs involved in the response to stripe rust pathogen infection have been identified [31]; in Arabidopsis, lncRNAs responsive to Fusarium oxysporum infection were revealed and played an important role in antifungal immunity [38]; and in tomato, lncRNAs have been identified as endogenous target mimics for miRNAs in response to Tomato yellow leaf curl virus infection [39], and a lncRNA (lncRNA16397) that conveys resistance to Phytophthora infestans by co-expressing glutaredoxin has also been detected [40].However, in studies of PaWB, lncRNAs that play important roles in the Paulownia-phytoplasma interaction have not yet been identified.
In this study, we used next-generation deep-sequencing technologies to identify PaWB-related genes and lncRNAs among healthy Paulownia tomentosa (PT), PaWB-infected P. tomentosa (PTI), and PT and PTI treated with 100 mg L −1 rifampicin (PT-100 and PTI-100, respectively).Our results may help to clarify the molecular mechanisms of PaWB, and accelerate the process of preventing the disastrous effects of PaWB.

Plant Materials
All the biological materials used in this study were obtained from the Institute of Paulownia, Henan Agricultural University, China.The following four groups of P. tomentosa seedlings were set up: PT, PTI, PT-100, and PTI-100.The cultivation and rifampin treatment procedures were as described by Fan et al. [13].The terminal buds from three individual plants were combined to form one biological replicate, and at least three biological replicates were used for each treatment.

Library Construction and Sequencing
Total RNA was extracted from the terminal buds of P. tomentosa using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and treated with RNase-free DNase (Promega, Madison, WI, USA) according to the manufacturer's instructions.The quality and quantity of total RNA were measured using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).Total RNA was used to construct RNA sequencing (RNA-seq) libraries using a TruSeq Stranded Total RNA preparation kit with RiboZero Plant (Illumina Inc., San Diego, CA, USA) according to the manufacturer's instructions.The libraries were sequenced on an Illumina HiSeq 2000 platform by the Beijing Genomics Institute, Shenzhen, Guang dong, China.

Genome Mapping and Transcript Assembly
Reads with more than 10% N (unable to determine base information), with adapter sequence, or of low quality were removed from the raw reads to obtain clean reads.To remove ribosomal RNA reads, the clean reads were mapped to the SILVA rRNA database (http://www.arb-silva.de/)using the short read alignment software SOAP2 (BGI, Shenzhen, China).Finally, the clean reads were mapped to the Paulownia genome sequence (http://paulownia.genomics.cn)using TopHat2 [41].Cufflinks 2.0.0 was used to assemble the mapped reads [42].Transcripts that were no less than 200 nucleotides in length were selected for further analysis.

LncRNAs Identification
Based on the assembled transcripts annotated to the Paulownia genome sequence, transcripts encoding known proteins were identified and the unknown transcripts (<200 bp) were excluded.Non-coding transcripts were identified based on a coding potential score of less than 0, as calculated by the Coding Potential Calculator (CPC; http://cpc.cbi.pku.edu.cn).The candidate lncRNAs were classified into several categories, according to their genomic locations and previous descriptions of Sun et al. [43].

Genomic Characterization and Function Prediction of lncRNAs
We determined the distribution of the identified lncRNA on the Paulownia chromosomes according to the description of Li et al. [44].To explore lncRNA conservation, all the lncRNA sequences identified in this study were aligned against the CANTATAdb (http://yeti.amu.edu.pl/CANTATA) using BLASTN (E-value < 1×10 −5 ).To identify lncRNAs that may act as miRNA precursors, we aligned the lncRNA sequences to the miRNA precursor sequences in miRBase (Release 21, http://www.mirbase.org/).The secondary structures of the lncRNA and miRNA precursor sequences were predicted using RNAfold in the Vienna RNA package (http://rna.tbi.univie.ac.at/).To explore the homology of the lncRNAs based on their consensus RNA secondary structures, we classified the predicted lncRNAs into different non-coding RNA families using INFERNAL [45].
We used an algorithm to search for potential cis target genes within a 10-kb region upstream and downstream of the identified lncRNAs.Potential trans target genes were detected based on sequence complementarity and RNA duplex energy prediction to assess the impact of lncRNA binding on complete mRNAs.BLAST was used to determine sequence complementarity (identify ≥ 95%, E-value < 0.05), and RNAplex was used to calculated the energy of potential mRNA-lncRNA interactions.The RNAplex parameters were set as e < −30.The criteria used for the prediction of potential trans targets are described in previous studies [36,46,47].

Identification of PaWB-Related Genes and lncRNAs
All the sequencing data were uploaded to SRA database with accession number SRP117715.Expression levels of genes and lncRNAs were measured as fragments per kilobase of exon model per million mapped reads (FPKM).Differentially expressed genes (DEGs) and differentially expressed lncRNAs (DELs) were identified according to the methods described in previous studies [36,40].
To identify genes and lncRNAs related to PaWB, we made comparisons among the four libraries (Figure 1).The following comparison schemes were considered:

LncRNAs Identification
Based on the assembled transcripts annotated to the Paulownia genome sequence, transcripts encoding known proteins were identified and the unknown transcripts (<200 bp) were excluded.Non-coding transcripts were identified based on a coding potential score of less than 0, as calculated by the Coding Potential Calculator (CPC; http://cpc.cbi.pku.edu.cn).The candidate lncRNAs were classified into several categories, according to their genomic locations and previous descriptions of Sun et al. [43].

Genomic Characterization and Function Prediction of lncRNAs
We determined the distribution of the identified lncRNA on the Paulownia chromosomes according to the description of Li et al. [44].To explore lncRNA conservation, all the lncRNA sequences identified in this study were aligned against the CANTATAdb (http://yeti.amu.edu.pl/CANTATA) using BLASTN (E-value < 1×10 −5 ).To identify lncRNAs that may act as miRNA precursors, we aligned the lncRNA sequences to the miRNA precursor sequences in miRBase (Release 21, http://www.mirbase.org/).The secondary structures of the lncRNA and miRNA precursor sequences were predicted using RNAfold in the Vienna RNA package (http://rna.tbi.univie.ac.at/).To explore the homology of the lncRNAs based on their consensus RNA secondary structures, we classified the predicted lncRNAs into different non-coding RNA families using INFERNAL [45].
We used an algorithm to search for potential cis target genes within a 10-kb region upstream and downstream of the identified lncRNAs.Potential trans target genes were detected based on sequence complementarity and RNA duplex energy prediction to assess the impact of lncRNA binding on complete mRNAs.BLAST was used to determine sequence complementarity (identify ≥ 95%, E-value < 0.05), and RNAplex was used to calculated the energy of potential mRNA-lncRNA interactions.The RNAplex parameters were set as e < −30.The criteria used for the prediction of potential trans targets are described in previous studies [36,46,47].

Identification of PaWB-Related Genes and lncRNAs
All the sequencing data were uploaded to SRA database with accession number SRP117715.Expression levels of genes and lncRNAs were measured as fragments per kilobase of exon model per million mapped reads (FPKM).Differentially expressed genes (DEGs) and differentially expressed lncRNAs (DELs) were identified according to the methods described in previous studies [36,40].
To identify genes and lncRNAs related to PaWB, we made comparisons among the four libraries (Figure 1).The following comparison schemes were considered:    We applied nested-PCR to detect phytoplasma in seedlings according to the method described by Fan et al. [13].Total RNA was reverse transcribed into cDNA and used to measure the expression of lncRNAs by quantitative real-time PCR (qRT-PCR).The qRT-PCRs were performed on a DNA Engine Opticon machine (MJ Research, Waltham, MA, USA) using a LightCycler FastStar DNA master SYBR Green I kit (Roche Diagnostics, Mannheim, Germany).The primers were designed using Primer Express 3.0 (Applied Biosystems, Stockholm, Sweden) and the specificity of primer pairs was checked by sequencing the PCR products.All qRT-PCR amplifications were carried out in triplicate, with the standard reaction program.The specificity of the amplified fragments was checked using the generated melting curve.The generated real-time data were analyzed using the Opticon Monitor Analysis Software 3.1 (Bio-Rad, Hercules, CA, USA) tool and standardized to the levels of 18S RNA using the 2 −∆∆Ct method [48].The primers used for the qRT-PCRs are listed in Table S1.

Changes in Morphology of the Four Samples
Compared with PT (Figure 2a), PTI (Figure 2b) showed typical PaWB symptom, such as witches' broom, and yellowing and relatively small leaves.In the rifampicin treated-samples (PT-100 and PTI-100), PT-100 (Figure 2c) showed no significant differences compared with PT.PTI-100 (Figure 2d) showed asymptomatic morphology.We applied nest-PCR to detect phytoplasma DNA in the four samples, and found that phytoplasma DNA was detected only in PTI (Figure S1).represents the union of two groups, except their intersection.Same represents the intersection of two groups.

Nested-PCR and Quantitative RT-PCR
We applied nested-PCR to detect phytoplasma in seedlings according to the method described by Fan et al. [13].Total RNA was reverse transcribed into cDNA and used to measure the expression of lncRNAs by quantitative real-time PCR (qRT-PCR).The qRT-PCRs were performed on a DNA Engine Opticon machine (MJ Research, Waltham, MA, USA) using a LightCycler FastStar DNA master SYBR Green I kit (Roche Diagnostics, Mannheim, Germany).The primers were designed using Primer Express 3.0 (Applied Biosystems, Stockholm, Sweden) and the specificity of primer pairs was checked by sequencing the PCR products.All qRT-PCR amplifications were carried out in triplicate, with the standard reaction program.The specificity of the amplified fragments was checked using the generated melting curve.The generated real-time data were analyzed using the Opticon Monitor Analysis Software 3.1 (Bio-Rad, Hercules, CA, USA) tool and standardized to the levels of 18S RNA using the 2 −ΔΔCt method [48].The primers used for the qRT-PCRs are listed in Table S1.

Changes in Morphology of the Four Samples
Compared with PT (Figure 2a), PTI (Figure 2b) showed typical PaWB symptom, such as witches' broom, and yellowing and relatively small leaves.In the rifampicin treated-samples (PT-100 and PTI-100), PT-100 (Figure 2c) showed no significant differences compared with PT.PTI-100 (Figure 2d) showed asymptomatic morphology.We applied nest-PCR to detect phytoplasma DNA in the four samples, and found that phytoplasma DNA was detected only in PTI (Figure S1).

High-Throughput Sequencing and Differentially Expressed Genes Analysis
Four cDNA libraries were constructed using the samples, and sequenced on an Illumina HiSeq 2000 platform.Among the clean reads, we identified 26,155, 27,017, 26,334, and 26,826 unique mRNAs from the PT, PTI, PT-100, and PTI-100 libraries, respectively (Table 1 and Table S2).A comparative summary of the mRNAs in the four libraries is provided as a Venn diagram (Figure S2a); 24,282 unique mRNAs were common among the four libraries, whereas 367, 337, 229, and 320 unique mRNAs were only presented in PT, PTI, PT-100, and PTI-100 libraries, respectively.In the PT/PT-100, PTI/PTI-100, and PT/PTI comparisons, we detected 1491 DEGs (706 upregulated and 785 downregulated), 1802 DEGs (725 upregulated and 1077 downregulated), and 3744 DEGs (2208 upregulated and 1536 downregulated), respectively (Table S2).Based on our comparison scheme (Figure 1), 2481 specific DEGs from comparison 1 (PT/PT-100) and comparison 2 (PTI/PTI-100) were identified (comparison 3), and 1063 common DEGs from comparison 3 and comparison 4 (PT/PTI) were PaWB-related genes.(Figure S3a and Table S2).Based on our analysis, we identified 3434, 3504, 3520, and 3479 lncRNAs in the PT, PTI, PT-100, and PTI-100 libraries, respectively (Table S3).A comparative summary of the mRNAs in the four libraries is provided as a Venn diagram (Figure S2b); 3102 unique lncRNAs were common among the four libraries, whereas 10, 2, 1 and 3 unique lncRNAs were only presented in PT, PTI, PT-100, and PTI-100 libraries, respectively.A total of 3693 lncRNAs were revealed from the four libraries and classified into four categories: 41 into classcode "j", novel long noncoding isoforms with at least one splice junction shared with the reference genes; 73 into classcode "o", exonic overlap with a known transcript; 578 into classcode "x", exonic overlap with reference on the opposite strand; and 3001 into classcode "u", intergenic lncRNAs (Figure S4a).The INFERNAL analysis results showed that 493 of the unique lncRNAs belonged to 155 conserved lncRNAs families (Table S4).About half of the lncRNAs belonged to miRNA families, such as MIR159, MIR169, and MIR171, and five lncRNAs belonged to two small nucleolar RNA families (SNORD25 and SNORD27).Some non-plant-specific lncRNA families (X inactive specific transcript, small cajal body-specific RNA 8, fungal small nucleolar RNA, bacterial small transcript non-coding 150) were also found.Almost all the lncRNAs were distributed among the chromosomes and did not show significant chromosome location preferences (Figure 3a).The lengths of the lncRNAs were usually shorter than the lengths of the mRNAs, and more than 10% of all lncRNAs were ≤300, 300-400, or 400-500 nucleotides long, whereas about 10% of the mRNA were ≥3000 nucleotides in length (Figure 3b).Most lncRNAs have only exons, while most mRNAs have more than one exon (Figure 3c).All the unique lncRNAs were searched against the genomes of Arabidopsis thaliana, Oryza sativa, Glycine max, Selaginella, Chlamydomonas, Physcomitrella, Amborella, potato, Vitis vinifera, and Zea mays using BLAST.The results showed that only a small number of the P. tomentosa lncRNAs were conserved across the ten species; the highest number on lncRNAs matched V. vinifera lncRNAs, likely because Paulownia and V. vinifera are the most closely related in evolution (Table S5).

Prediction of lncRNA Candidate Target Genes
Because lncRNAs play important roles in regulating gene expression, identification and analysis of their candidate target genes could help understand their functions.A total of 6746 lncRNA-mRNA pairs were identified (4175 were cis-regulated and 2571 were trans-regulated, Table S6).In total, we detected 5333 candidate target genes for 3222 lncRNAs; among them, 3218 were cis target genes for 2290 lncRNAs, and 2401 were trans target genes for 2571 lncRNAs.We found that one lncRNA could have more than one target gene, and one target gene could be targeted by one or more lncRNAs.Among the 3222 lncRNAs, 1246 had a single candidate target gene, while nine candidate target genes were predicted for each of two lncRNAs (TCONS_00019479 and TCONS_00030009).Furthermore, among the 5333 candidate target genes, 4313 were targeted by one lncRNA, while the other two target genes (TCONS_00013387 and TCONS_00019818) were targeted by more than one lncRNA (each of the target genes were targeted by 11 lncRNAs) (Table S6).
We compared the 138 candidate targets of the 80 PaWB-related lncRNAs with the 1063 PaWB-related DEGs, and identified 12 PaWB-related candidate target genes among the DEGs that were regulated by nine PaWB-related lncRNAs (Table 2 and Figure 4).The KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis results indicated that the 12 candidate target genes were mapped to eight pathways, namely "Phenazine biosynthesis", "RNA transport", "Cysteine and methionine metabolism", "Amino sugar and nucleotide sugar metabolism", "Glycerolipid metabolism", "Monoterpenoid biosynthesis".

Quantitative RT-PCR Analysis of DEGs and DELs
To confirm the expression of the P. tomentosa lncRNAs and their response to PaWB, qRT-PCR analysis was performed to verify the results of the high-throughput RNA-seq.We randomly selected seven lncRNAs and six mRNAs for qRT-PCR validation.The results demonstrated that, except TCONS_00029942, the qRT-PCR results were consistent with those from the RNA-seq analysis (Figure 5).This correlation confirmed that the RNA-seq data were reliable.To confirm the expression of the P. tomentosa lncRNAs and their response to PaWB, qRT-PCR analysis was performed to verify the results of the high-throughput RNA-seq.We randomly selected seven lncRNAs and six mRNAs for qRT-PCR validation.The results demonstrated that, except TCONS_00029942, the qRT-PCR results were consistent with those from the RNA-seq analysis (Figure 5).This correlation confirmed that the RNA-seq data were reliable.

Discussion
LncRNAs play key roles in various biological pathways in animals and plants.There have been recent advances in RNA-seq which, when combined with genome-wide mapping, have resulted in the identification of lncRNAs in plants.The identification of plant lncRNAs has opened up the investigation of novel regulatory pathways in plants.For example, 7655 lncRNAs from control and gibberellin-treated Populus tomentosa have been identified, indicating that lncRNAs may participate in auxin signal transduction and synthesis of cellulose and pectin, and may influence growth and wood properties.In the present study, we investigated transcriptomic changes in PT, PTI, PT-100,

Discussion
LncRNAs play key roles in various biological pathways in animals and plants.There have been recent advances in RNA-seq which, when combined with genome-wide mapping, have resulted in the identification of lncRNAs in plants.The identification of plant lncRNAs has opened up the investigation of novel regulatory pathways in plants.For example, 7655 lncRNAs from control and gibberellin-treated Populus tomentosa have been identified, indicating that lncRNAs may participate in auxin signal transduction and synthesis of cellulose and pectin, and may influence growth and wood properties.In the present study, we investigated transcriptomic changes in PT, PTI, PT-100, and PTI-100 seedlings, and systematically identified genes and lncRNAs associated with PaWB.Our results represent the first comprehensive analysis of lncRNAs in Paulownia.The identified lncRNAs will be useful for other woody plants researchers and provide an important resource for future functional genomics and regulatory expression studies.
The Paulownia lncRNAs had fewer exons and were shorter than mRNAs in length, which is consistent with other studies [35,44].BLAST searches of the Paulownia lncRNAs against the genomes of species such as Arabidopsis thaliana, Oryza sativa, and Glycine max revealed that the majority of Paulownia lncRNAs were not conserved with lncRNAs in other species.The low conservation suggested that lncRNAs may undergo rapid evolution.It is not surprising that lncRNAs are not well conserved for various reasons.First, lncRNAs are not constrained by codon usage.Second, although lncRNAs may possess short conserved motifs, these short motifs are not easily identifiable by local alignment software such as BLAST.Third, some lncRNAs may interact directly with RNA-binding proteins through conserved secondary structures.
Xyloglucan (XG) is the main hemicellulose of the primary cell wall in dicotyledons.Xylosyl transferases (XT, EC:2.4.1.207)are a family of enzymes that catalyze xyloglucan endotransglucosylase and/or xyloglucan endohydrolase [49].XTs can cut xyloglucan chains and transfer the fragment with the new reducing end either to another XG or to water, and play an important role in cell wall remodeling, affecting plant growth and development [50].Transglycosylation can contribute to both cell wall assembly and cell wall loosening [51].In this study, XT was the predicted target gene of lncRNA (TCONS_00022242), and it was upregulated in PT/PTI.XT may be related to witches' broom and hyperplasia because of its roles in cell wall assembly and loosening.
sugar and nucleotide sugar metabolism (map00520), Neomycin, kanamycin and gentamicin biosynthesis (map00524), Glycerolipid metabolism (map00561), Zeatin biosynthesis (map00908), Biosynthesis of plant secondary metabolites (map01060), and Metabolic pathways (map01100).Therefore, changes in SQD1 expression may influence UDP-glucose content and the pathways in which UDP-glucose is involved.Plants constitutively synthesize a wide variety of secondary metabolites to aid fitness by preventing pathogen invasion [53].Among these metabolites, terpenoids are frequently described as natural products that are active against pathogens [54].Neomenthol dehydrogenase (MNR, PAU012048.1),which participates in the monoterpene synthesis of neomenthol and isomenthol, is a monoterpenoid dehydrogenase [55].In Capsicum annuum, the menthone reductase gene CaMNR1 regulates plant defenses against a broad spectrum of pathogens [56].In our study, MNR was the target of a PaWB-related lncRNA (TCONS_00026472), and was upregulated in PT/PTI.The change in MNR expression may be in response to PaWB-phytoplasmas infection.
DNA methylation plays an essential role in regulating plant development [57].DNA methyltransferase (Dnmt), an important enzyme for DNA methylation, is not only associated with DNA methylation, but is also linked to many important biological activities, including cell proliferation and senescence [58].In previous studies, we found that PaWB may be related to changes in DNA methylation [59,60].In this study, a PaWB-related gene, DNA (cytosine-5-)-methyltransferase (PAU029256.2,Dnmt1), was identified as the target gene of a PaWB-related lncRNA (TCONS_00032704), and was upregulated in PT/PTI.After infection, the change in Dnmt1 expression might influence DNA methylation and result in the phenotype of infected Paulownia.

Conclusions
Although the functions of most lncRNAs are still unknown, the evidence suggests that lncRNAs may play regulatory roles by interacting with RNA, DNA, and protein-coding genes.For example, it has been reported that lncRNAs can act as cisor trans-regulators of protein-coding gene expression in animals and plants.The lncRNAs detected in this study will lay the groundwork for future functional studies of lncRNAs in Paulownia.We revealed 26,155 and 27,017 genes, and 3434 and 3504 lncRNAs in healthy and PaWB-infected P. tomentosa, respectively.A total of 6746 lncRNA-mRNA pairs were identified (4175 were cis-regulated and 2571 were trans-regulated).In PT/PTI, 3744 DEGs (2208 upregulated and 1536 downregulated) and 750 (200 upregulated and 550 downregulated) were identified.According to our analysis, 1063 mRNAs and 110 lncRNAs were related to PaWB, in which 12 mRNAs were candidate target genes of nine lncRNAs.To understand the specific biological role of lncRNAs and their regulatory mechanisms in Paulownia, future research should include functional analyses of the lncRNA candidate target genes using overexpression, CRISPR-Cas9, or RNA interference gene silencing strategies.

( 1 )
DEGs and DELs in PT/PT-100 may be related to the influence of rifampicin; (2) DEGs and DELs in PTI/PTI-100 may be related to the influence of rifampicin and PaWB; (3) Differences between comparisons 1 and 2 may exclude DEGs and DELs related to the influence of rifampicin; (4) DEGs and DELs in PT/PTI may be related to PaWB; (5) The common DEGs and DELs between comparisons 3 and 4 may be related directly to PaWB.

( 1 )
DEGs and DELs in PT/PT-100 may be related to the influence of rifampicin; (2) DEGs and DELs in PTI/PTI-100 may be related to the influence of rifampicin and PaWB; (3) Differences between comparisons 1 and 2 may exclude DEGs and DELs related to the influence of rifampicin; (4) DEGs and DELs in PT/PTI may be related to PaWB; (5) The common DEGs and DELs between comparisons 3 and 4 may be related directly to PaWB.

Figure 1 .
Figure 1.Comparison schemes of the four samples.PT represents the healthy wild-type sample of P. tomentosa, PTI represents the sample of phytoplasma infected PT.PT-100 represents the sample of 100 mg L −1 rifampicin treated PT, PTI-100 represents the sample of 100 mg L −1 rifampicin treated PTI.Different represents differentially expressed genes/long non-coding RNAs (lncRNAs).Difference

Figure 1 .
Figure 1.Comparison schemes of the four samples.PT represents the healthy wild-type sample of P. tomentosa, PTI represents the sample of phytoplasma infected PT.PT-100 represents the sample of 100 mg L −1 rifampicin treated PT, PTI-100 represents the sample of 100 mg L −1 rifampicin treated PTI.Different represents differentially expressed genes/long non-coding RNAs (lncRNAs).Difference represents the union of two groups, except their intersection.Same represents the intersection of two groups.

Forests
. Nested-PCR and Quantitative RT-PCR

ForestsFigure 3 .
Figure 3. Characteristics of Paulownia lncRNAs.(a) Distribution of lncRNAs along each chromosome; (b) Transcript length distribution of lncRNAs and mRNAs; (c) Number of exons per transcripts for lncRNAs and mRNAs.

Figure 4 .
Figure 4. Characteristics of PaWB-related lncRNAs and mRNAs.(a) Expression profiling in PT, PTI PT-100, and PTI-100.The data were based on lg(FPKM); (b) LncRNA-mRNA networks.The blue polygons represent mRNAs, the yellow polygons represent lncRNAs.FPKM, fragments per kilobase of exon model per million mapped reads.

Figure 4 .
Figure 4. Characteristics of PaWB-related lncRNAs and mRNAs.(a) Expression profiling in PT, PTI PT-100, and PTI-100.The data were based on lg(FPKM); (b) LncRNA-mRNA networks.The blue polygons represent mRNAs, the yellow polygons represent lncRNAs.FPKM, fragments per kilobase of exon model per million mapped reads.

Figure 5 .
Figure 5. Expression pattern confirmation by qRT-PCR.(a) Changes in the relative expression levels as determined by qRT-PCR.Standard error of the mean for three technical replicates is represented by the error bars.Samples marked with various letters show a significant difference at p < 0.05; (b) Changes in the relative expression levels as determined by RNA-seq.The expression level of PT was set to 1.

Figure 5 .
Figure 5. Expression pattern confirmation by qRT-PCR.(a) Changes in the relative expression levels as determined by qRT-PCR.Standard error of the mean for three technical replicates is represented by the error bars.Samples marked with various letters show a significant difference at p < 0.05; (b) Changes in the relative expression levels as determined by RNA-seq.The expression level of PT was set to 1.

Table 1 .
Statistical data of the RNA-Seq reads.