All samples in this study were tissue culture seedlings, healthy P. tomentosa seedlings (HP) and P. tomentosa seedlings infected with PaWB phytoplasma (PIP), both of them were obtained from the Institute of Paulownia. These seedlings were cultured for 30 days on 1/2 MS medium before being clipped (Murashige and Skoog, 1962). After clipping, the terminal buds of 1.5-cm PIP were transferred into 100-mL triangular flasks containing 20 ml 1/2 MS culture medium supplemented with 0, 20, or 60 mg·L⁻¹ MMS (PIP, PIP-20, and PIP-60 respectively). The terminal buds of 1.5-cm HP were transferred into the same medium without MMS as the control. All the seedlings were cultured in a darkroom at 20°C for 5 days. After that, they were cultured at 25 ± 2°C under 130 μmol·m⁻² s⁻¹ light intensity with a 16/8 h light /dark photoperiod for 25 days. Then, the terminal buds of 1.5-cm seedlings in the four treatment groups were sheared, immediately frozen in liquid nitrogen, and stored at −80°C. For each treatment, 60 terminal buds were planted in 20 flasks. Each treatment was performed in triplicate.

The four samples (HP, PIP, PIP-20, and PIP-60) were ground to powder in liquid nitrogen, then these different powder were added the moderate lysis buffer with 1 mM PMSF and 2 mM EDTA, the details of the protein extraction method was described as Tang et al. (2016).

To label the peptides obtained from the HP, PIP, PIP-20, and PIP-60 samples with iTRAQ reagent, 100 μg total protein from each sample solution was digested with Trypsin Gold (Promega, Madison, WI, USA). The peptides in the four samples were labeled with isobaric tags, the details method of iTRAQ-labeled according to the method of Meng et al. (2014). The labeled peptide mixtures were pooled and dried by vacuum centrifugation. Strong cation exchange (SCX) chromatography was performed with a LC-20AB HPLC Pump system (Shimadzu, Kyoto, Japan). The labeled peptide mixture contained proteins extracted from the PIP, PIP-20, PIP-60, and HP samples, including two technical replicates of isotopic labeling. Each biological replicate consisted of 12 terminal buds of the same stage from four different samples.

The labeled peptide mixture was firstly resuspended and then the supernatant was loaded on a LC-20AD nano HPLC system (Shimadzu), then the data acquisition was performed with a TripleTOF 5600 System (AB sciex, concord, ON, Canada) fitted with a Nanospray III source (AB sciex) and a pulled quartz tip as the emitter (New Objectives, Woburn, MA). The mass spectrometer was operated with a RP of ≥30,000 full-width half-maximum for time-of-flight mass spectrometry scans. The method of the peptides were separated according to Qiao et al. (2012).

Raw data files were firstly converted into Mascot generic format files using Proteome Discoverer 1.2 software (Thermo Fisher Scientific, Bremen, Germany), then these data were matched against the transcriptome database containing 105,812 all-unigene sequences (Fan et al., 2015a,b). For protein identification, to reduce the probability of false peptide identification, only peptides at the 95% confidence interval, as defined by Mascot probability analysis, greater than “identity” were considered as identified, and each confident protein involved at least one unique peptide. For protein quantitation, it was required that a protein contained at least two unique spectra.

To further identify the functions of the identified proteins, the Blast2GO program was run against the NCBI non-redundant protein sequence database to assign them Gene Ontology (GO) terms. The Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Clusters of Orthologous Groups (COG) databases were used to classify and group the identified proteins. Quantitative protein ratios were weighted and normalized by the median ratio in Mascot. Proteins with a fold change ≥ 1.2 and p < 0.05 were determined as differentially abundant proteins.

In this section, the correlation analysis of the transcriptome and proteome data was carried out. Firstly, the proteins identified by iTRAQ were matched to 105,812 unigene sequences of the transcriptome. After that, the proteins which have homologs in the transcritome were screened out. Then we calculated the number of the correlated protein in four different comparisons, including the PIP vs. PIP-20 comparison, PIP-20 vs. PIP-60 comparison, HP vs. PIP-60 comparison and HP vs. PIP comparison, and the correlated proteins can be considered as the identified proteins which have been expressed at the transcription level. Subsequently, the quantitative analysis of the correlated proteins was performed, and the criterion of the quantitation protein is the unique peptide ≥2. According to the results of the quantitative analysis, we further picked out the differentially expressed proteins (fold change ≥ 1.2 and p < 0.05) in different comparisons. At last, the number of the genes that corresponding to the correlated proteins and the number of differential expressed proteins were figured out. The other important correlation analysis is to compare the abundance levels of the differentially expressed proteins and their corresponding genes from the transcriptome sequencing.

A comparison scheme of PaWB-related proteins was described by Liu et al. (2013). The proteins were searched against the NCBI non-redundant protein sequence database using the Blast2GO program to assign GO terms, and against the COG and KEGG databases to classify them to pathways and functional groups. To deeply analyze the function of PaWB-related proteins, protein BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to confirm the function of these proteins, then the functional classification of them was performed.

RNA samples from HP, PIP, PIP-20, and PIP-60 were extracted with Trizol (Sangon, Shanghai, China). First-strand cDNAs of the RNA samples were synthesized using an iScript cDNA synthesis kit (Bio-Rad, Hercules, CA, USA). From the transcriptome data (Fan et al., 2015a,b), we randomly chosen 12 unigene sequences that corresponded to the differentially abundant proteins for qRT-PCR, primers for each of the unigenes were designed using PrimerPremier 5.0 software (PREMIER Biosoft International, Palo Alto, Calif). The results were analyzed using the 2^−ΔΔCt method (Livak and Schmittgen, 2000). The expression level of each unigene was analyzed in three replicates. All the primers sequences used for the qRT-PCR are listed in Supplementary Table S1. Statistical analysis was performed using SPASS 19.0 (SPASS, Inc., Chicago, IL, USA).

We used an iTRAQ approach to analyze the total proteins in four P. tomentosa samples (HP, PIP, PIP-20, and PIP-60; Table 1). Protein homologs were identified against the transcriptome sequencing data from the same samples. We obtained 386,933 total spectra by combination analysis, and 21,948 of the spectra were matched using Mascot software. Among the 21,948 matched spectra, 17,811 were unique. The total number of detected proteins was 2,051 (Supplementary Table S2). The accuracy of the mass spectrometry was <2 ppm, and the database search peptides matching error was below 0.05 Da. We carried out two biological replicates to identify the total proteins and found that 1,597 of the proteins were identified in repeat 1, 1,657 were identified in repeat 2, and 1,203 were common to the two sets (Figure 1). The producibility of the proteomic analysis showed that the results for the duplicates were of high quality (Figure 2), which indicated that the proteome results were reliable.

The functions of the 2,051 identified proteins were classified using COG. The proteins were classified into 23 COG subgroups (Figure 3), among which “general function prediction only,” “post-translational modification, protein turnover, chaperones,” and “energy production and conversion” were significantly enriched and contained the largest numbers of proteins, while “RNA processing and modification,” “defense mechanisms,” and “cell motility” contained the smallest numbers of proteins (Supplementary Table S3). The 2,051 proteins were also classified into the three GO categories: biological process, molecular function, and cellular component (Figure 4). Under biological process, “metabolic process” contained the largest number of proteins; under cellular component, “cell” and “cell part” contained the largest numbers of proteins; and under molecular function, “catalytic activity” contained the largest number of proteins (Supplementary Table S4). In the KEGG pathways analysis (Supplementary Table S5), the 2,051 proteins were mapped to 117 biological pathways of which “metabolic pathways,” “biosynthesis of secondary metabolites,” and “ribosome” were the most enriched metabolism pathways, while “linoleic acid metabolism,” “monoterpenoid biosynthesis,” and “circadian rhythm-plant” were the least enriched pathways. Based on these results, we concluded that most of the annotated proteins may affect metabolic processes, post-translational modification, protein turnover, chaperones, and catalytic activity.

Differentially abundant proteins (DAPs) were defined as proteins having a 95% confidence level, p < 0.05, and fold changed >1.2. According to these criteria, the abundance of 1,172 of the 2,051 proteins did not differ significantly in the PIP vs. HP, PIP-60 vs. PIP-20, PIP-20 vs. PIP, and PIP-60 vs. HP comparisons, while the abundance of 879 of the proteins changed significantly (Table 2). Among the 879 DAPs, 223 were identified in PIP vs. HP, 204 were identified in PIP-60 vs. PIP-20, 278 were identified in PIP-20 vs. PIP, and 174 were identified in PIP-60 vs. HP. Further analysis indicated that most of the DAPs were correlated to several metabolism pathways, namely biosynthesis of secondary metabolites, photosynthesis, carbon fixation in photosynthetic organisms, carotenoid biosynthesis, nitrogen metabolism, porphyrin, and chlorophyll metabolism, glyoxylate and dicarboxylate metabolism, pyrimidine metabolism, porphyrin, and chlorophyll metabolism, and plant–pathogen interaction.

A correlation analysis of the transcriptome and proteome data was carried out, and the numerical relationships between the correlated proteins and the total proteins were calculated in the different comparisons. In the PIP vs. PIP-20 comparison, 2,050 proteins were correlated with the transcriptional data; among them, 1,046 proteins were quantitated and 278 of them were differentially abundant. In the PIP-20 vs. PIP-60 comparison, 2,048 proteins were correlated with the transcriptional data; among them 1,051 proteins were quantitated and 204 of them were differentially abundant. In the HP vs. PIP-60 comparison, 2,049 proteins were correlated with the transcriptional data; among them, 1,083 proteins were quantitated and 174 of them were differentially abundant. In the HP vs. PIP comparison, 1,073 proteins were quantitated and 223 of them were differentially abundant. All the numbers of these proteins corresponding genes were shown in Table 3.

Based on these results, we compared the abundance levels of the proteins and the expression levels of their corresponding genes. The following patterns were detected: (1) the same trend between the corresponding mRNA and protein levels; (2) the opposite trend between the corresponding mRNA and protein levels: (3) no change in protein levels, but changes in gene expression levels; (4) changes in protein levels, but no change in gene expression levels. These results show that the protein abundance level is not necessarily reflected in the expression level of the corresponding mRNA. Similar results have been reported in phytoplasma-infected Mexican lime (Monavarfeshani et al., 2013).

According to the comparison scheme of PaWB-related proteins, a total of 43 DAPs were obtained that changed their abundance levels in response to phytoplasma infection (Table 4), then the function of 43 DAPs were further classified in the COG, GO, and KEGG databases. The DAPs were classified into 13 COG categories (Supplementary Table S6), among which the three largest groups were “energy production and conversion,” “post-translational modification, protein turnover, chaperones,” and “general function prediction only.” The 43 DAPs were assigned 38 GO terms (Supplementary Table S7); 20 were under biological process, 11 were under cellular component, and five were under molecular function. In the KEGG pathway enrichment analysis, the 43 DAPs were mapped to 27 pathways (Supplementary Table S8), among which the two largest groups were “metabolic pathways” and “biosynthesis of secondary metabolites,” while the two smallest groups were “plant-pathogen interaction” and “phenylpropanoid biosynthesis,” indicating that changes in some metabolites may be important in the response of P. tomentosa seedlings to PaWB infection.

To confirm and dig out the function of the 43 DAPs, protein sequence BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to confirm the functional of the 43 DAPs in other plant. The result showed 40 DAPs functions have been confirmed in other plants (Table 4). Of these proteins, six proteins were classified in the section named as post-translational modification, protein turnover, chaperones, such as, 20S proteasome beta subunit, ubiquitin extension protein, subunit of chloroplasts chaperonin, peptidyl-Prolyl Isomerase, Co-chaperone GrpE family protein and fk506- and rapamycin-binding protein; six proteins participated in photosynthesis pathway, including chloroplast thylakoid lumen protein (TLP), photosystem I subunit E-2 (PSI-E2), plastid-specific ribosomal protein 4 (PSRP4), the component of the light harvesting complex (LHCI) and alpha carbonic anhydrase (CA), and the subunit E of photosystem I (PSI-E2); three proteins were associated with energy metabolism pathway, such as ATP synthase gamma subunit, NADP-malic enzyme (NADP-ME), and the subunit of glyceraldehyde-3-phosphate dehydrogenase (GAPC); two Ca²⁺-binding protein involved in cell signal transduction pathways; five proteins, included aldo/keto reductase, cinnamyl-alcohol dehydrogenase, fasciclin-like arabinogalactan protein, germin-like protein, and pectin methylesterase took part in plant defense. In addition these DAPs with functions, 3 of PaWB-related proteins were annotated as unknown function, their function still need to verify.

To confirm the results of the transcriptome sequencing data, 12 gene sequences corresponding to PaWB-related proteins were selected randomly for the qRT-PCR assays (Figure 5). The results showed that the relative expressions of six of the 12 DAP-related unigenes were significantly down-regulated in the process of the morphometric recovery, while five were significantly up-regulated. These 11 DAP-related expression patterns were consistent with the results of the iTRAQ LC-MS/MS analysis. For one of DAPs, no related mRNA expression was detected. This discrepancy may be attributed to post-transcriptional processing, post-translational processing and modification, or to different rates of degradation of mRNA and protein.

The onset of the yellowish leaves symptom in phytoplasma-infected plants is one of the most obvious symptoms of witches' broom. This symptom is closely related to changes in photosynthesis, and many studies have reported that photosynthesis is arrested in phytoplasma-infected plants, as evidenced by a decrease in pigment molecules such as chlorophyll b levels soon after infection (Scarpari et al., 2005). Accordingly, phytoplasma infection leads to a reduction in the expression of genes that encode photosystem components, further suppressing photosynthesis (Ji et al., 2009; Margaria et al., 2012; Nejat et al., 2015). Similar results have been reported in the interaction of Paulownia and phytoplasma, where several important protein components associated with photosynthesis were found to be significantly differentially abundant. These components include two proteins in photosystem I, six proteins in photosystem II, two proteins in light harvesting complex I, pyruvate kinase, and the pI6.8 24-kDa protein (Fan et al., 2003; Mou et al., 2013; Liu et al., 2013).

In this study, we identified six DAPs participated in photosynthesis, including TLP, two PSI-E2s, PSRP4, LHCI, and CA, which had reduced abundance levels in response to phytoplasma infection.

In the phytoplasma-infected seedlings, we detected one DAP, PSRP 4 was down regulated, which belongs to photosystem II. Tiller et al. (2012) reported that Arabidopsis RNAi PSRP 4 mutant appeared light-green phenotype, smaller mesophyll cells, and the chlorophyll content, chloroplast translation, plastocyanin content, and maximum quantum efficiency of photosystem II (F_V/F_M) significantly reduced, and the mutants also showed strongly reduced accumulation of three complexes (photosystem II, cytochrome b6f complex, photosystem I), which were responsible for photosynthetic electron transport.

There were three DAPs belong to Photosystem I (PSI), including TLP, LHCI, and PSI-E2. It is well known that thylakoid lumens in both integral and soluble membrane proteins, and TLP take part in light harvesting and electron transfer in the photosynthetic chain. In this study, the TLP was down-regulated in the diseased seedlings. Evidences showed that the reduction of TLP not only damaged the thylakoid electron flow, but also decreased chlorophyll a content and further slacked the photosynthesis (Wang et al., 2016; Zhu, 2016), this result is consistent with the previous RNA-seq results, which showed that a reduction in photosynthesis of phytoplasma-infected plants could be attributed to the loss of several thylakoid membrane proteins (Mou et al., 2013). At the same time, the expression of LHCI was reduced, which severely affected the light reactions in Photosystem I (PSI) (Bressan et al., 2016). Interestingly, we also detected a decreased abundance of PSI-E2. Ihnatowicz et al. (2007) demonstrated that the rate of thylakoid electron transfer were not affected in the absence of PSI-E2. Therefore, we speculate that the decreased abundance of the thylakoid lumenal protein in the PaWB-infected P. tomentosa plants might partly explain the yellowish leaves symptom.

Another DAP associated with photosynthesis was CA, which had reduced abundance level in the phytoplasma-infected seedlings. This protein plays a vital role in the early event of photosynthesis. Zhurikova et al. (2016) indicated that the reduction of the activity of CA affected the effective quantum yield of photosystem II, thus, its low abundance might weaken photosynthesis in the phytoplasma-infected plants.

The hormonal imbalance in plants is the most factors for the abnormal symptom. Hoshi et al. (2009) discovered that TENGU can interfere the expression of auxin-related gene, which was the main reason for the plant proliferation and dwarf symptom. In the interaction of paulownia and phytoplasma, Mou et al. (2013) indicated that the key enzymes of cytokinin biosynthesis including isopentenyl diphosphate isomerase and isopentenyltransferase were significantly induced in the phytoplasma-infected Paulownia; Liu et al. (2013) demonstrated that the auxin efflux carrier 5NG4 were down-regulated in the phytoplasma-infected seedlings, which resulted in auxin accumulation in the paulownia.

In our proteome analysis, the abundance of two proteins that related to gibberellin (GA) and brassinosteroid (BR) were changed after phytoplasma infection: for example, the abundance of the glycine rich protein GRP was increased in phytoplasma-infected seedlings. It has been shown that over-abundance of GRP reduced the activity of several intermediates of the GA biosynthetic pathway, and transgenic Arabidopsis (AtGRP) displayed reduced length of the vegetative stem (Löhr et al., 2014). The abundance of the membrane steroid binding protein (MSBP) also increased after phytoplasma infection. MSBP encodes a 220-amino acid protein that can bind to 24-epi-brassinolide and negatively regulates BR signaling, and overexpression of MSBP could inhibit cell elongation through downregulating cell elongation related genes, and result in the reduced cell elongation and shortened hypocotyl (Yang et al., 2005). Based on the results above, we concluded that the differentially abundant GRP and MSBP affected GA and BR biosynthesis, which were linked to the short internodes and dwarf morphology.

Phytoplasma genomes lack many metabolic pathways, making them unlikely to synthesize nucleotides, and they lack the key gene for ATP synthesis, they may need to import ATP from their host (Carle et al., 2011). Therefore, in phytoplasma-infected plants, the host must supply enough energy for intense phytoplasma growth. In phytoplasma-infected mulberry, the host's energy production was decreased (Ji et al., 2009). Oshima et al. (2001) demonstrated that, in onion, phytoplasma strongly depended on the host's glycolysis pathway to obtain more energy, while, in phytoplasma-infected Mexican lime, energy production by the TCA cycle was significantly enhanced (Monavarfeshani et al., 2013).

In this study, we observed that the abundance levels of some key proteins related to energy metabolism in the phytoplasma-infected P. tomentosa, such as ATP synthase gamma subunit, NADP-ME, and GAPC. In the phytoplasma-infected P. tomentosa plants, the ATP synthase γ subunit showed increased abundance. The ATP synthase γ subunit is closely related to the ATP production, and the up-regulation of ATP synthase γ subunit can enhance the photosynthetic rates in stressed plant (Budak et al., 2013). It is well known that the activity of chloroplast CF0-CF1-ATP synthase can be regulated by the light-dark. Kohzuma found that the modified γ subunit by mutating conserved D211V, E212L, and 226L acidic amino acids only altered the light induced regulation, but not metabolism induced regulation, and verified the function mutants of γ subunits affected photosynthetic ATP synthesis (Kohzuma et al., 2012, 2013). Therefore, in the photosynthesis pathway, the infected paulownia ATP might exhibit a tendency of increasing.

The other pathway of energy production is TCA cycle, an important pathway of energy production, in our study, the expression of DAP like NADP-ME showed decreased abundance, which plays a central role in the metabolite flux through the TCA cycle. In the Arabidopsis nadp-me knockout mutant, the contents of the TCA metabolites like 2-oxoglutarate and succinate were increased, while that of citrate and fumarate were decreased (Tronconi et al., 2008), implying that differential expression of NADP-ME had little impact on TCA cycle, and energy production by TCA cycle was not the necessary pathway for phytoplasma propagation in the host. This result was in disagreement with that in phytoplasma-infected Mexican lime (Monavarfeshani et al., 2013), the difference may be depend on the type of phytoplasmas and their hosts.

Another pathway of energy production is glycolytic pathway. GAPC, an important cytosolic enzyme which catalyzes a key reaction in glycolysis, was up-regulated in the phytoplasma infected seedlings. Evidence showed that GAPC play important roles in energy and carbohydrate metabolites, and Guo et al. (2014) reported that knockout or overexpression of GAPCs caused significantly changes in the content of intermediates of glycolytic pathway, the ratios of ATP/ADP, and NAD(P)H/NAD(P).

Compared with the DAPs functions of the three energy production pathways above, we supposed that glycolytic and photosynthesis may be the main energy provider for Paulownia response to the phytoplasma. However, details of the actual process need to be further researched.

The involvement of cell signal transduction pathways in plant—pathogen interactions directly or indirectly affected plant development (Ranjan et al., 2015). In Paulownia—phytoplasma interactions, many cell signal-related proteins have been identified from the transcriptome level, including plant hormones, calcium-dependent protein kinases, MAP kinases, receptor-like kinases, LRR receptor-like serine/threonine-protein kinases, L-ascorbate peroxidase, (S)-2-hydroxy-acid oxidase, and nitricoxide synthase (Mou et al., 2013; Liu et al., 2013; Cao et al., 2014a,b; Fan et al., 2015b). In our proteome analysis, we detected several DAPs mainly involved in Ca²⁺, ROS and plant defense signal pathway.

Calcium, an ubiquitous secondary messenger, plays an important role in all aspects of cell function, which has been regarded as versatile intracellular signal (Steinhorst and Kudla, 2014). Calcium-binding proteins is the component of the calcium-signaling pathway. In this study, two Ca²⁺-binding proteins were increased after phytoplasma infection. Zhou et al. (2013) reported that pathogen can manipulate the host Ca²⁺ signaling machinery to benefit their own life cycles. At the same time, evidence also has been documented that the calcium-binding protein could be as an effector taking part in plant defense (Ye et al., 2017), demonstrating that calcium-signaling pathway plays a central role in the interaction of paulownia and phytoplasma.

The abundance of the germin-like protein GLP, which belongs to a large ubiquitous family of plant glycoproteins, was also increased in the PaWB-infected seedlings. GLP plays a vital role in plant defense (Rietz et al., 2012). Further, it has been reported that high levels of GLP may initiate oxidative bursts in pathogen-infected plants and elevate the levels of endogenous reactive oxygen species (ROS) (Beracochea et al., 2015). Increasing evidence has shown that ROS (and redox signals) not only can induce cell damage, but also can act as reactive substrates to kill pathogens (Nejat et al., 2015). However, the ROS production was regulated by antioxidant protein disulfide isomerase (PDI), which was induced after phytoplasma infection. PDI, a major ER protein, usually acts as a molecular chaperone and component of signal-transduction pathways. It has been reported that PDI can limit potential cell damage by ROS generation after pathogen infection in plants (Stolf et al., 2011), which has been implicated in the complex interplay of defense-related signaling pathways.

Keeping the above views in mind, our results indicate that the response of P. tomentosa to PaWB infection involved several interconnected signaling pathways, including Ca²⁺ and ROS-mediated signaling, and plant defense signaling, which coordinate the plant's response to phytoplasma. Regulation of cell signal pathways need different proteins that are induced by interactions with susceptible as well as resistant hosts, which play positive or negative roles in the Paulownia response to phytoplasma, depending on the speed and intensity of the interaction responses.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.