Lasting consequences of psyllid (Bactericera cockerelli L.) infestation on tomato defense, gene expression, and growth

The tomato psyllid, Bactericera cockerelli Šulc (Hemiptera: Triozidae), is a pest of solanaceous crops such as tomato (Solanum lycopersicum L.) in the U.S. and vectors the disease-causing pathogen ‘Candidatus Liberibacter solanacearum’. Currently, the only effective strategies for controlling the diseases associated with this pathogen involve regular pesticide applications to manage psyllid population density. However, such practices are unsustainable and will eventually lead to widespread pesticide resistance in psyllids. Therefore, new control strategies must be developed to increase host-plant resistance to insect vectors. For example, expression of constitutive and inducible plant defenses can be improved through selection. Currently, it is still unknown whether psyllid infestation has any lasting consequences on tomato plant defense or tomato plant gene expression in general. In order to characterize the genes putatively involved in tomato defense against psyllid infestation, RNA was extracted from psyllid-infested and uninfested tomato leaves (Moneymaker) 3 weeks post-infestation. Transcriptome analysis identified 362 differentially expressed genes. These differentially expressed genes were primarily associated with defense responses to abiotic/biotic stress, transcription/translation, cellular signaling/transport, and photosynthesis. These gene expression changes suggested that tomato plants underwent a reduction in plant growth/health in exchange for improved defense against stress that was observable 3 weeks after psyllid infestation. Consistent with these observations, tomato plant growth experiments determined that the plants were shorter 3 weeks after psyllid infestation. Furthermore, psyllid nymphs had lower survival rates on tomato plants that had been previously psyllid infested. These results suggested that psyllid infestation has lasting consequences for tomato gene expression, defense, and growth.


Background
The tomato psyllid (or potato psyllid), Bactericera cockerelli Šulc (Hemiptera: Triozidae), is a major pest of solanaceous crops such as tomato (Solanum lycopersicum L.) and potato (S. tuberosum) in the U.S. [8]. The psyllid is native to the Southwestern U.S. and Northern Mexico [12,49,55,64] but has only recently become an important agricultural pest when it was discovered that B. cockerelli vectors the disease-causing pathogen 'Candidatus Liberibacter solanacearum' (Lso) [43]. Lso is a fastidious bacterial pathogen associated with zebra chip disease in potato as well as other diseases in solanaceous crops [37,41]. Today, Lso is considered a major pathogen of crops worldwide [20,63]. Currently, the only effective strategies for controlling the diseases associated with Lso involve calendar application of insecticide [8,42]. However, these strategies are unsustainable. Multiple reports indicate neonicotinoid resistance is increasing in certain B. cockerelli populations [5,45,50]. Since vector-borne disease systems are faced with the rapid evolution of pesticide resistance, major efforts have been made to develop novel solutions based on selectively breeding plants for improved host-plant resistance or genetically manipulating plants and insects for the purpose of disrupting disease transmission [3,4,34,35,57,66]. For example, disease transmission can be disrupted by manipulating the host or vector's genes associated with key molecular pathways that facilitate the movement of pathogens from host to vector and vice versa [1,32]. Such genetic manipulations can be accomplished through direct transformations or artificial selection, but these toolkits require certain a priori genomic information. Therefore, in order to pursue psyllid control strategies that manipulate the host plant's molecular pathways, the current study identifies the genes involved in the transcriptomic response of tomato plants to psyllid infestation.
The current study focuses on an insect-plant relationship, however the experiments described are informed by Lso disease development. Specifically, diseases caused by Lso are characterized by long latent periods. Indeed, symptoms in tomato and potato typically start developing 3 weeks after infection [33,40,43,59]. Logically, studies of Lso infection are conducted a few weeks or even months after plants are infested with psyllids and subsequently infected with Lso. To avoid the confounding effects of psyllid herbivory, some studies entirely divorce the effect of vector infestation by transmitting the pathogen from one host-plant to another via grafting [13,59]. Furthermore, the rate of Lso infection and disease development are independent of psyllid density [52]. Thus, the long-term effects of psyllid infestation on tomato plant biology and gene expression are divorced from Lso research and are still unknown. This is important knowledge gap considering psyllids are known to cause phenotypic changes in solanaceous crops under heavy infestation (≥100 insects per plant), a condition called 'psyllid yellows' [7,60]. Typically, studies of Lso infection have involved a single control group of plants that have not been exposed to either the psyllid vector or the Lso pathogen. Then, controls will be compared to plants exposed to both the psyllid and Lso. This practice has been acceptable because psyllid-responsive expression changes in plants are expected to be relatively unimportant compared to Lso challenge. Although this experimental design has been invaluable for characterizing Lso disease severity and psyllid transmission efficacy, an unintended consequence is the knowledge gap regarding the lasting consequences of psyllid infestation on tomato plant health. The molecular interaction between host plant and insect vector is especially important because plants have several long-term responses to insect damage that can impact their lifetime health, reproduction, and defense.
Plants undergo physiological, transcriptomic, or epigenetic changes which allows them to mount a stronger and faster responses to secondary challenges by previously perceived threats. This is called defense 'priming' [10,21,30,39]. Priming is a common phenomenon that has been studied in several plant species in response to bacteria, fungi, and chewing insects [11,24,61,68]. Furthermore, plants can remain immunologically primed for the rest of their lives or even across generations [47,53,62]. Therefore, it is reasonable to hypothesize that tomato plants deploy similar long-term defenses against psyllids postinfestation and that these changes have lasting consequences for tomato survival, growth, and development. In fact, the lasting the consequences of uninfected psyllid infestation were previously observed (but not quantified) in a study by Mendoza Herrera et al. [40].
The current study evaluated the persistent transcriptomic and physical responses of tomato plants to psyllid infestation. This was accomplished by comparing the transcriptomes of uninfested plants to plants that had been infested 3 weeks prior. Second, tomato plant growth was tracked across time to test the relationship between plant growth/development and immune response to psyllid infestation. This experimental design allows for the identification of genes involved in the tomato plant's response to psyllid infestation and whether these genes were associated with improved defense against psyllids. Third, psyllid populations were monitored for the number of eggs laid and nymphal survival when reared on previously uninfested tomato plants (controls) compared to psyllids reared on previously infested plants.

1-Transcriptomic analysis
Illumina sequencing of tomato cDNA libraries produced 95.2 million reads that met FastQC quality control criteria (i.e., Phred quality scores > 35). The average number of reads obtained from uninfested plants (17.4 ± 0.6 million) did not significantly differ from psyllid-infested ones (18.0 ± 0.4 million) (t-value = − 0.68; P = 0.25). HISAT2 alignment analysis showed that 96.3 ± 0.1% of all reads from uninfested plants and 96.2 ± 0.3% of all reads from psyllid-infested plants mapped to vSL3.0 of the S. lycopersicum genome (Supplementary Table 2); these alignment rates did not significantly differ (tvalue = 0.14; P = 0.45). The Ballgown analysis identified 362 differentially expressed genes (DEGs) between control and psyllid-infested plants (q-value < 0.01). These DEGs represented the pattern of systemic tomato plant gene expression following psyllid infestation. Gene expression patterns were visualized with a heatmap comparing the fold change (Z-Score) for each gene between samples (Fig. 1); Z-scores based on deviations from the average fpkm (fragments per kilobase per million read) value for a given gene. Additionally, a dendrogram (Fig. 1) and a principal component analysis (PCA, Fig. 2) comparing fpkm values across genes and samples were used to visualize relative similarities in gene expression across samples. Both the dendrogram and the PCA geometries suggested that the overall pattern of gene expression was consistent within each treatment, where per-gene fpkm values were most similar within treatment and most different between treatments. Furthermore, the PCA showed that the first principal component strongly separated the fpkm values of psyllid-infested plants from uninfested plants and accounted for 84.1% of the total variance in fpkm values, meaning the greatest differences in gene expression between samples were the differences between infested and uninfested plants. for details). Tomato plant DEGs were assigned to one or more of the following broader categories: Defense response to biotic or abiotic stress (55 DEGs), transcription/translation (50 DEGs), photosynthesis (35 DEGs), molecular signaling (33 DEGs), molecular transport (31 DEGs), reproduction (27 DEGs), protein phosphorylation/ubiquitination (26 DEGs), cellular turnover (23 DEGs), sugar metabolism (20 DEGs), ion transport/ homeostasis (16 DEGs), auxin signaling (9 DEGs), and cell wall biosynthesis/metabolism (6 DEGs) (Tables 1, 2, 3 and 4). RT-qPCR corroborated the relative expression levels in tested genes: Results showed that the unchanged PIP2-4 (Solyc06g011350.2) was expressed at similar levels in both uninfested (1.13 ± 0.01) and psyllid-infested plants (1.12 ± 0.01; t-value = 0.69, P = 0.26). The upregulated DRIP2 (Solyc06g084040.2) was expressed at significantly lower level in control (1.15 ± 0.02) compared to psyllid infested (1.36 ± 0.03; t-value = − 6.54, P < 0.01). The downregulated LON2 (Solyc04g080860.1) was expressed at significantly higher levels in control (1.45 ± 0.11) compared to psyllid infested (1.01 ± 0.06; t-value = 4.04, P < 0.01). Lastly, the downregulated D27 (Solyc08g008630.2) was expressed at significantly higher levels in control (1.26 ± 0.08) compared to psyllid infested (0.83 ± 0.08; t-value = 4.10, P < 0.01).   (CC). The left axis represents the -log 10 (p adj ) likelihood that a given MF, BP, or CC would be associated with a random selection of Arabidopsis genes. Circle sizes represent the relative number of times a given MF, BP, or CC appears among analyzed genes. In general, expression changes occurred throughout the cell and were most likely to be involved with cellular processes, metabolism, photosynthesis, response to stimulus, and biological regulation. Labels above, connected to arrows, or adjacent to circles describe specific the MF, BP, or CC associated with each circle; some labels have been omitted due to redundancy Table 1 The 55 tomato plant DEGs associated with defense response to abiotic and biotic stress. DEGs were sorted by log2-fold change (log2FC). These DEGs were identified in the transcriptome analysis comparing psyllid-infested and uninfested tomato plants 3 weeks after infestation (P < 0.01). NCBI Blast searches were used to identify Gene IDs and protein products in tomatoes as well as their homologs in other species. Specifically, the expression changes in 44 genes (80%, in bold) would have resulted in improvements to plant defense pathways. These DEGs were related to defense against insect damage, microbial infection, programmed cell death, salt stress, and drought. Simultaneously, 11 DEGs, especially those related to the hypersensitive response, underwent expression changes that would have resulted in impairments to plant defense pathways          Table 2 The 50 tomato plant DEGs associated with transcription and translation. DEGs were sorted by log2-fold change (log2FC). These DEGs were identified in the transcriptome analysis comparing psyllid-infested and uninfested tomato plants 3 weeks after infestation (P < 0.01). NCBI Blast searches were used to identify Gene IDs and protein products in tomatoes as well as their homologs in other species. Specifically, the expression changes in 44 genes (88%, in bold) would have resulted in improvements to transcription/translation pathways. These DEGs were related to post-translational modifications, miRNA processing, and gene silencing (Continued)  Table 2 The 50 tomato plant DEGs associated with transcription and translation. DEGs were sorted by log2-fold change (log2FC). These DEGs were identified in the transcriptome analysis comparing psyllid-infested and uninfested tomato plants 3 weeks after infestation (P < 0.01). NCBI Blast searches were used to identify Gene IDs and protein products in tomatoes as well as their homologs in other species. Specifically, the expression changes in 44 genes (88%, in bold) would have resulted in improvements to transcription/translation pathways. These DEGs were related to post-translational modifications, miRNA processing, and gene silencing (Continued)  Table 3 The 35 tomato plant DEGs associated with molecular signaling. DEGs were sorted by log2-fold change (log2FC). These DEGs were identified in the transcriptome analysis comparing psyllid-infested and uninfested tomato plants 3 weeks after infestation (P < 0.01). NCBI Blast searches were used to identify Gene IDs and protein products in tomatoes as well as their homologs in other species. Specifically, the expression changes in 28 genes (85%, in bold) would have resulted in improvements to molecular signaling pathways. These DEGs were related to protein phosphorylation and mobilization to the vacuole

Psyllid development experiments
The psyllid development experiments showed that psyllids laid a statistically similar number of eggs on plants that had been previously infested (36.6 ± 13.4, n = 28) and uninfested plants (48.8 ± 12.1, n = 27) (t-score = − 0.71, P = 0.24). Also, the rate of egg hatching was similar between psyllids raised on previously infested plants (88.3 ± 6.7%) compared to psyllids raised on uninfested plants (89.1 ± 2.8%) (n = 55; t-score = 0.04, P = 0.48). In contrast, the same experiments showed that nymphs had a significantly lower survival rate when reared on previously psyllid-infested plants (71.9 ± 6.0%) compared to nymphs reared on uninfested plants (85.4 ± 3.7%) (tscore = − 1.89, P = 0.03). These differences, though, were Table 3 The 35 tomato plant DEGs associated with molecular signaling. DEGs were sorted by log2-fold change (log2FC). These DEGs were identified in the transcriptome analysis comparing psyllid-infested and uninfested tomato plants 3 weeks after infestation (P < 0.01). NCBI Blast searches were used to identify Gene IDs and protein products in tomatoes as well as their homologs in other species. Specifically, the expression changes in 28 genes (85%, in bold) would have resulted in improvements to molecular signaling pathways. These DEGs were related to protein phosphorylation and mobilization to the vacuole (Continued)  Table 4 The 33 tomato plant DEGs associated with photosynthesis. DEGs were sorted by log2-fold change (log2FC). These DEGs were identified in the transcriptome analysis comparing psyllid-infested and uninfested tomato plants 3 weeks after infestation (P < 0.01). NCBI Blast searches were used to identify Gene IDs and protein products in tomatoes as well as their homologs in other species. Specifically, the expression changes in only 7 genes (20%, in bold) would have resulted in improvements to photosynthesis. Simultaneously, 28 DEGs, especially those related to response to light stimulus and photorespiration, underwent expression changes that would have resulted in impairments to photosynthesis  Table 4 The 33 tomato plant DEGs associated with photosynthesis. DEGs were sorted by log2-fold change (log2FC). These DEGs were identified in the transcriptome analysis comparing psyllid-infested and uninfested tomato plants 3 weeks after infestation (P < 0.01). NCBI Blast searches were used to identify Gene IDs and protein products in tomatoes as well as their homologs in other species. Specifically, the expression changes in only 7 genes (20%, in bold) would have resulted in improvements to photosynthesis. Simultaneously, 28 DEGs, especially those related to response to light stimulus and photorespiration, underwent expression changes that would have resulted in impairments to photosynthesis (Continued)      only apparent after nymphs had spent 3-5 days on previously-psyllid infested plants. These results suggest that tomato plants responded to psyllid infestation by mounting an immune response that made them less suitable hosts for psyllid nymphs 3 weeks after the first infestation (Fig. 5).

Discussion
Transcriptomic analysis of S. lycopersicum leaves showed that 362 genes were differentially expressed in tomato plants 3 weeks after psyllid infestation, suggesting that a week-long infestation by a small number of B. cockerelli had lasting consequences for gene expression in tomato plants ( Figs. 1 and 2). Homologs of the DEGs were associated with 1) defense against abiotic and biotic stress, 2) transcription/translation, 3) molecular signaling, and 4) photosynthesis (Tables 1, 2  were consistent with the results of the transcriptome analysis by demonstrating that psyllid infestation had lasting consequences for tomato plant growth (Fig. 4) and defense (Fig. 5). Specifically, the growth experiments demonstrated that tomato growth was stunted by psyllid infestation while the psyllid development experiments demonstrated that tomato plants that had been previously challenged by psyllids were less suitable hosts for nymphs. Among the DEGs identified in the transcriptome analysis, 55 were homologs of genes associated with defense against biotic and abiotic stress (Table 1). For example, regulatory protein NPR3 (NPR3; Solyc02g069310.2) is a substrate-specific adapter of an E3 ubiquitin-protein ligase complex which mediates the ubiquitination and subsequent proteasomal degradation of target proteins, and consequently regulates the basal defense response to pathogens [69]. Since expression of NPR3 was significantly up-regulated (P = 0.001) in tomato plants 3 weeks after psyllid infestation, its associated defensive pathway was likely increased. Furthermore, NPR3 is involved in defense against insects, therefore its up-regulation may have been a consequence of plant defensive priming and/or the crosstalk between the jasmonic acid and salicylic acid pathways [16,46]. Recently, a study performed in citrus plants showed that exposure to Asian citrus psyllids for 14 and 150 days resulted in induction of NPR1 and a delay in plant growth compared to the unfed plants. This effect was not detected after 7 days. The authors concluded that the prolonged exposure (~150 days) of citrus to Asian citrus psyllid feeding suppressed plant immunity and inhibited growth, probably through the salicylic acid signaling pathway [28]. Based on the functional characterization of Arabidopsis homologs, the expression changes observed in 80% stress-related DEGs would have likely coincided with increased responsiveness to abiotic and biotic stressors (see Table 1 for citations).
A subset of 50 DEGs were homologs of genes involved in transcription and/or translation (Table 2). For example, RNA-binding KH domain-containing protein RCF3 (RCF3; Solyc03g034200.2) is a negative regulator of osmotic stressinduced gene expression [67]. Since the expression of RCF3 was down regulated in tomato plants 3 weeks after psyllid infestation (P = 0.001), stress responsive gene expression would have increased. This interpretation is supported by the up regulation of genes such as homeobox-leucine zipper protein ATHB-12 (ATHB-12; Solyc01g096320.2), phospholipase D alpha 4 (PLDALPHA4; Solyc03g121470.2), and inactive poly [ADP-ribose] polymerase RCD1 (RCD1; Solyc08g005270.2). Furthermore, the expression profile changes observed in 88% of DEGs related to transcription/translation likely coincided with increased transcription/translation (see Table 2 for citations). Similarly, a subset of 35 genes were homologs of genes that function in molecular signaling (Table 3). In fact, the most common functional categories associated with DEGs were cellular processing and intracellular signaling ( Fig. 3;  Supplementary Figure 3). Together, these results suggest tomato plants were still active in responding to the psyllid threat 3 weeks after psyllids were last sensed by the plant.
A set of 33 DEGs were homologs of genes involved in photosynthesis (Table 4). For example, RNA polymerase sigma factor sigA (SIGA; Solyc03g097320.2) controls the transcription of the psaA gene and modulates photosystem stoichiometry, meaning its down regulation in tomato plants would have likely led to impaired photosynthesis after psyllid infestation [14,19]. Furthermore, the expression changes in 26 (80%) DEGs related to photosynthesis would have likely also coincided with impaired photosynthesis. In support of this observation, the long-term, deleterious effects of psyllid infestation on tomato plant growth were evidenced by the experiments that tracked tomato plant stem length after psyllid infestation. These experiments showed the growth rate in tomato plant stems slowed after psyllid infestation (Fig. 4). These results were consistent with our previous study that observed stunted growth in tomato plants after psyllid infestation [40]. In addition to stunted stem growth, other developmental processes were likely impacted by psyllid infestation. For example, 6 DEGs were homologs of genes involved in auxin signaling. Since auxin-related signaling has several effects on plant growth and orientation, expression changes in these genes may be related to the stunting observed in tomato plants after psyllid infestation. Changes to plant growth, development, and photosynthesis post-herbivory may be related to the molecular crosstalk that takes place between plant defensive pathways and plant growth/development pathways [23,26,54].
Although 251 DEGs were homologs of genes for which published characterizations were available, 111 DEGs (30.7%) lacked any supporting information. This means nearly a third of the lasting consequences of psyllid infestation on tomato gene expression remain unknown. Of these DEGs, 78 (70.3%) were up-regulated in psyllidinfested plants relative to controls, consistent with the general pattern observed across DEGs. Therefore, it is reasonable to hypothesize that many of these expression changes would also be related to stress response, translation/transcription, molecular signaling, and/or photosynthesis.
In conclusion, the results of this manuscript are the first to report the long-lasting effects of psyllid herbivory on plant gene expression and health. The transcriptomic and growth experiments demonstrated that tomato plants underwent expression changes that likely repressed growth and developmental pathways in favor of promoting the expression of a select number of genes which are likely involved in defense against psyllid challenge. The DEGs that improved defense may constitute the genes directly involved in the tomato's long-term response to psyllid challenge. This hypothesis is supported by the psyllid development experiments which showed psyllid nymphs had lower survival rates on psyllidinfested plants relative to uninfested plants (Fig. 5). The results presented in the current research showed that short exposures to small numbers of phloem feeding insects can have significant and lasting consequences for plant gene expression, growth, and defense. Alternatively, it is possible that the expression changes observed in tomato plants 3 weeks after psyllid infestation were a consequence of the accumulation of stress-related expression changes during psyllid infestation and sampling (with a razor blade). Continual stress can create negative feedback loops in stress-responsive genetic pathways [2]. This explanation is consistent with the overall deleterious impact of psyllid infestation observed in this study [9,22]. Future disease biology research should continue exploring the long-term effects that vectors have on their hosts independent of their associated pathogens. These results should also be taken into consideration for epidemiologic studies of diseases associated with Liberibacter and their psyllid vectors.

Insect source
B. cockerelli were field-collected from Weslaco, Texas in 2008 and used to establish laboratory colonies. Tomato psyllid colonies have since been maintained on tomato plants under a 16: 8-h (Light: Dark) photoperiod at room temperature (22 ± 2°C). The absence of Lso in these psyllid colonies was confirmed each month using the diagnostic PCR method previously described by Nachappa et al. [44]. Briefly, DNA from psyllids from the colony was extracted using the 10% CTAB method and subjected to PCR amplification of 'Candidatus Liberibacter solanacearum' 16S rDNA.

Plant material
Tomato plants, cultivar Moneymaker (Victory Seed Company; Molalla, OR), were grown from seed in Metro-Mix 900 (Sun Gro Horticulture, Agawam, MA) soil and individually transplanted to 10 × 10 cm square pots 4 weeks later. Plants were watered every other day and fertilized weekly according to the manufacturer's recommendation (Miracle-Gro® Water Soluble Tomato Plant Food; 18-18-21 NPK). All experiments were conducted at the same photoperiod (16: 8) and temperature (22 ± 2°C) used to rear psyllids.

Psyllid infestation and sample collection
Psyllid infestation were initiated when plants were 6 weeks old. Leaves branching below the apical meristem (i.e., leaves similar to the ones sampled for the transcriptome analysis) were caged with a small, white organza bag (amazon.com). Restricting psyllids to these leaves exposed them to systemic response of the plant to any prior infestation. Each bag either had no psyllids (control plants) or three adult male psyllids (psyllid-infested plants). Males were chosen to avoid the potentially confounding effect of oviposition on tomato gene expression. Seven days after infestation, caged tomato leaves were removed with a bleach-sterilized razor blade. Three weeks later, the top-most, fully developed leaf was sampled from each plant and immediately flash-frozen in liquid nitrogen. Samples were transferred to Eppendorf tubes and kept submerged under liquid nitrogen while ground with plastic, RNase-free pestles.

RNA purification, sequencing and bioinformatic analysis
Total RNA extraction was performed on leaf tissue harvested 3 weeks after psyllid infestation using the Plant RNeasy Mini Kit (Qiagen, Valencia, CA) following the manufacturer's protocol. Three biological replicates were sequenced per treatment (i.e., uninfested and psyllidinfested, six samples total). One fully-develop leaf and petiole were removed per biological replicate using sterilized razor blades. The top-most leaf was sampled to ensure that the gene expression changes observed were more likely to be associated with a plant systemic response. Samples were ground using sterilized plastic pestles. RNA samples were treated with RNase-Free DNase (Qiagen). Any remaining DNA was removed using the TURBO DNA-free™ Kit (Life Technologies, Carlsbad, CA). All remaining RNA was stored at − 80°C for downstream quantitative reverse transcription PCR (RT-qPCR) validation. The isolated RNA was submitted to the Texas A&M Genomics and Bioinformatic Service for quality analysis, library preparation, and sequencing.
For transcriptomic sequencing, cDNA libraries were developed using the TruSeq RNA Library Prep Kit v2 (Illumina®; San Diego, CA) following the manufacturer's protocol, generating 2 Х 150 bp read lengths. Libraries were multiplexed and sequenced on the Illumina PE HiSeq 2500 v4 platform. Sequence cluster identification, quality prefiltering, base calling, and uncertainty assessment were done in real time using Illumina's HCS 2.2.38 and RTA 1.18.61 software with default parameter settings. Library preparation, sequencing, and read processing were performed by the Texas A&M Genomics and Bioinformatic Service. The processed sequences were uploaded to the CyVerse Discovery Environment computational infrastructure [17] where bioinformatic analysis was performed using the HISAT2-StringTie-Ballgown RNA-Seq workflow [31]. Libraries reads were mapped to the S. lycopersicum genome (vSL3.0) using HISAT2. StringTie assembled hits to known transcripts based on the vITAG3.2 annotation and made non-redundant with StringTie-Merge. DEGs were identified using Ballgown. Genes were considered differentially expressed when comparative qvalues were below 0.01 [48]. DEG gene names were searched against the tomato genome database [15,25] as well as the PhytoMine search engine in Phytozome [18]. DEGs were assigned putative functions based on their homology with other plant genes with known function published in Ensembl Plants (version SL2.50) and the UniProt Knowledgebase [6]. Arabidopsis thaliana homologs of DEGs were uploaded to the NCBI Gene Expression Omnibus (GEO) functional genomics data repository in order to visual overrepresentation among molecular pathways using the g:Profiler functional profiler.

Transcriptome validation by RT-qPCR
To verify the results of the transcriptomic analysis, RT-qPCR analyses were performed on three genes differentially expressed in psyllid-infested plants: One putatively upregulated gene, an E3 ubiquitin-protein ligase that acts as a negative regulator of the response to water stress (Solyc06g084040.2 or DRIP2) [36] and two putatively downregulated genes, a peroxisomal protease potentially involved in drought stress response (Solyc04g080860.1 or LON2) and a chloroplastic Betacarotene isomerase D27 (Solyc08g008630.2 or D27) [38,65]. Since many of the regulatory genes differentially expressed in this study were involved in drought stress, an aquaporin (Solyc06g011350.2 or PIP2-4) that putatively underwent no regulatory change was selected as a control [29]. RT-qPCR experiments were conducted using RNA from the six sequenced tomato leaf samples (three per treatment) as well as six independently grown tomato plants (three per treatment), which were obtained by repeating the plant growth and infestation assays (three plants per treatment). This allowed for validation of the transcriptome results. An aliquot of 500 ng RNA was taken from each sample to develop cDNA libraries using the Verso™ cDNA Kit (Thermo Fisher Scientific, Waltham, MA), following the manufacturer's manual. The cDNA libraries were diluted to 1:5 prior to RT-qPCR. Each reaction consisted of 1.0 μL cDNA, 5.0 μL SensiFAST SYBR Hi-ROX mix (Bioline, Memphis, TN), 0.4 μL of each primer (400 nM), and 3.6 μL of molecular grade water. Primers were designed using Primer3 [56], which targeted exons within a DEG, had an optimal annealing temperature of 60.0-62.0°C, and generated 150 bp amplicons (Supplementary Table 1). RT-qPCR was performed in an Applied Biosystem QuantStudio 6 Flex system using the following parameters: 2 min at 95°C, followed by 40 cycles of 5 s at 95°C and 30 s at 60°C. The melting curve for each reaction was generated to assure amplicon specificity. All RT-qPCR reactions were performed in triplicate. Relative expression levels for each gene were analyzed using the 2 -ΔΔCT method [51] with glyceraldehyde 3-phosphate dehydrogenase (GADPH) as a reference gene [27]. Since expression levels did not assume normality, they were analyzed using the Mann-Whitney U ranked test in JMP® Version 13 (SAS Institute Inc., Cary, NC, 1989-2018).

Plant growth and psyllid development on previously infested and uninfested plants
Tomato plants were grown and treated using the same methods described above where 28 tomato plants were psyllid-infested and 27 plants were left uninfested. In order to minimize handling stress, plant growth was tracked using pictures taken 3 weeks after infestation to compare the total stem length of psyllid-infested plants to uninfested plants. Each picture included a 52 cm-long tray that served as a size standard. The total length (in pixels) of a tomato plant main stem was measured from the soil to the tip of the apical meristem using ImageJ1.X [58] and converted to centimeters using the length standard. This no-contact method of measurement was chosen to minimize plant wounding. Stem lengths were analyzed using a one-way student's t-test in JMP.
Three weeks after initial infestation, three female psyllids were transferred to a no-choice cage and allowed to oviposit on undamaged leaves of the tomato plants that had previously been psyllid-infested or uninfested. As before, psyllids were restricted to a single leaf inside an organza bag, using a different leaf than the one used during the initial infestation. This exposed them to plant systemic conditions. Three adult females were caged together in each bag; there was one bag per plant. After 48 h, psyllids were removed, and their eggs were counted. Eggs were left on their respective plants and allowed to hatch. Nymphs were counted every other day and left to develop into adults. Adults were collected as they emerged. Egg hatching and nymph survival rates were calculated for the psyllids reared on each plant. Additionally, initial egg number and nymphal survival rates were compared between psyllids reared on previously infested and uninfested plants. Since 100% of the nymphs that survived development also emerged, adult emergence rate was not compared. Egg number and nymph survival were analyzed using student's one-way ttests in JMP.