Integrative transcriptomics reveals association of abscisic acid and lignin pathways with cassava whitefly resistance

Background Whiteflies are a global threat to crop yields, including the African subsistence crop cassava (Manihot esculenta). Outbreaks of superabundant whitefly populations throughout Eastern and Central Africa in recent years have dramatically increased the pressures of whitefly feeding and virus transmission on cassava. Whitefly-transmitted viral diseases threaten the food security of hundreds of millions of African farmers, highlighting the need for developing and deploying whitefly-resistant cassava. However, plant resistance to whiteflies remains largely poorly characterized at the genetic and molecular levels. Knowledge of cassava-defense programs also remains incomplete, limiting characterization of whitefly-resistance mechanisms. To better understand the genetic basis of whitefly resistance in cassava, we define the defense hormone- and Aleurotrachelus socialis (whitefly)-responsive transcriptome of whitefly-susceptible (COL2246) and whitefly-resistant (ECU72) cassava using RNA-seq. For broader comparison, hormone-responsive transcriptomes of Arabidopsis thaliana were also generated. Results Whitefly infestation, salicylic acid (SA), jasmonic acid (JA), ethylene (ET), and abscisic acid (ABA) transcriptome responses of ECU72 and COL2246 were defined and analyzed. Strikingly, SA responses were largely reciprocal between the two cassava genotypes and we suggest candidate regulators. While susceptibility was associated with SA in COL2246, resistance to whitefly in ECU72 was associated with ABA, with SA-ABA antagonism observed. This was evidenced by expression of genes within the SA and ABA pathways and hormone levels during A. socialis infestation. Gene-enrichment analyses of whitefly- and hormone-responsive genes suggest the importance of fast-acting cell wall defenses (e.g., elicitor recognition, lignin biosynthesis) during early infestation stages in whitefly-resistant ECU72. A surge of ineffective immune and SA responses characterized the whitefly-susceptible COL2246’s response to late-stage nymphs. Lastly, in comparison with the model plant Arabidopsis, cassava’s hormone-responsive genes showed striking divergence in expression. Conclusions This study provides the first characterization of cassava’s global transcriptome responses to whitefly infestation and defense hormone treatment. Our analyses of ECU72 and COL2246 uncovered possible whitefly resistance/susceptibility mechanisms in cassava. Comparative analysis of cassava and Arabidopsis demonstrated that defense programs in Arabidopsis may not always mirror those in crop species. More broadly, our hormone-responsive transcriptomes will also provide a baseline for the cassava community to better understand global responses to other yield-limiting pests/pathogens. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-023-04607-y.


Background
Cassava (Manihot esculenta) is a hardy tuber crop that feeds 800 million people in over 100 countries worldwide, including small shareholder farmers of sub-Saharan Africa, who rely on the crop for subsistence [1].However, since the 1990's, superabundant whitefly (Bemisia tabaci) populations have devastated cassava yields in Africa [2,3].As phloem-feeders, whiteflies deplete photosynthates, deposit mold growth-promoting honeydew and transmit viral diseases, together slowing cassava growth and root production [4,5].In South America, the whitefly Aleurotrachelus socialis also significantly reduces cassava yields (60-80%) [6,7].A potent whitefly resistance was discovered in the Ecuadorian genotype ECU72, manifesting as nymph mortality, which prevents adult emergence and blocks population expansion [8,9].ECU72's resistance also extends to four other whitefly species, providing a promising control strategy [10][11][12].A better understanding of the genetic basis of this resistance will inform breeding efforts to select for identified resistance traits in African cassava.
At the genetic/molecular level, plant defense against biotic stressors begins with recognition of the attacker [13].Generic molecular signatures derived from pathogens (pathogen-associated molecular patterns, PAMPs), insect herbivores (herbivore-associated molecular patterns, HAMPs), or from 'debris' from damaged host cells (damage-associated molecular patterns, DAMPs) are recognized by extracellular receptors to elicit a defense response via defense signals [14,15].Such signals are numerous and include the two major defense hormones SA (salicylic acid) and JA (jasmonic acid), as well as emerging defense signals ethylene (ET), abscisic acid (ABA) and reactive oxygen species (ROS) among others.Multiple signals are often necessary for basal immunity and resistance [16].
While whiteflies affect the yields of numerous crops [17], strong, fast-acting resistance specific to whiteflies has only been identified in cassava and alfalfa resulting in reduced egg deposition or mortality of early-stage nymphs [8,18,19].Broad spectrum to moderate resistance to whiteflies has also been found in Solanaceous crops and their wild relatives among others [20][21][22].Previous studies have provided some insight into the plant-defense signals elicited during responses to whitefly infestation.In whitefly-resistant lines of tomato, cotton, and cabbage, SA, JA/ET and ABA responses, respectively, are primarily elicited by whiteflies [23][24][25].
In whitefly-susceptible plant species, SA levels increase and SA or JA/ET signaling can increase [26][27][28].In addition, basal resistance to whiteflies can also be promoted by JA or ABA [29][30][31].To date, the molecular mechanisms underlying resistance remain largely uncharacterized.Our previous work showed that during infestation, whitefly-susceptible cassava genotypes elicit Pathogenesis-related (PR) genes involved in cell wall processes [32].A complementary metabolomics study, which assessed cell wall phenolics in cassava genotypes ECU72 and COL2246 suggested constitutive whitefly-resistance in ECU72 to involve cell wall reinforcement [33].
To obtain a more global understanding of resistant and susceptible responses of cassava to whiteflies at the transcript level, we define the transcriptomes of whiteflyresistant ECU72 and whitefly-susceptible COL2246 in response to A. socialis infestation, SA, JA, ET, and ABA.Hormone-responsive transcriptomes in cassava were analyzed at the genome-scale and as subsets of defensehormone-pathway genes.These data sets were compared to Arabidopsis and also integrated with whitefly-responsive transcriptomes to reveal the genotype-dependent biological processes occurring during infestation.Together, our integrative transcriptomics approach, supported by metabolomics analyses, identified an association of higher ABA levels and induction of ABA-and lignin-pathway and cell-wall-related genes with whiteflyresistance in ECU72.Additionally, higher SA levels and induction of SA-pathway genes were associated with whitefly-susceptibility in COL2246.

A fast transcript-level response to SA revealed in whitefly-resistant cassava
To define the defense hormone-responsive transcriptomes of cassava, we profiled the response of cassava genotypes ECU72 (whitefly-resistant) and COL2246 (whitefly-susceptible) to SA, JA, ET, and ABA at seven times post treatment (hpt) (Fig. 1).Differentially expressed genes (DEGs) were identified by temporal comparisons within a genotype ("temporal" DEGs, tDEGs) (Fig. 1a).The number of tDEGs varied substantially by hormone (Fig. 1a; Additional file 1: Tables S1-S8).While the responses to ET or ABA were similar, ECU72 and COL2246 had distinct temporal responses to SA and JA.Notably, as early as 0.5 to 1 h post SA or JA treatment, 5-to 65-fold more tDEGs were induced and repressed in ECU72 versus COL2246 (Fig. 1a; Additional file 2: Table S13).Comparison of ECU72 and COL2246 transcriptomes at each time point identified "genotype" DEGs (gDEGs) (Fig. 1b; Additional file 1: Tables S9-S12).SA elicited 1.5-2.8-foldmore gDEGs than any other hormone treatment and had a bimodal temporal program of gene expression with gDEGs peaking at 0.5 and 12 hpt (Fig. 1b; Additional file 2: Table S14).
Principal component analyses (PCAs) performed with reads detected in the SA, JA, ET, and ABA time courses revealed that while all responses had a temporal component, only SA and JA responses showed a clear distinction between genotypes (Fig. 1c; Additional file 3: Figure S1).PCAs resolved the two temporal phases, here called early and late.All COL2246 hormone responses followed the same early (0-2 hpt) and late (4-12 hpt) phases with a return to basal expression at 24 hpt.In contrast, ECU72's late phase of SA responses initiated earlier and lasted longer (2-24 hpt).Furthermore, PCA and heatmaps of temporal gene expression showed reciprocity in ECU72's and COL2246's responses to SA (Fig. 1c,d).Additionally ECU72's late phase clustered with the 0-h time point indicating a significant difference in constitutive expression of JA-responsive genes between the genotypes (Fig. 1c).
Pearson correlation analyses using genes detected during hormone treatments highlighted additional interactions between the hormone pathways in ECU72 and COL2246 (Additional file 4: Figure S2; Additional file 5: Table S15).Most striking, negative correlation values suggested antagonism between ECU72's late SA and early JA responses with ECU72's responses to ET and ABA.In contrast, most all other hormone responses were positively correlated between genotypes, suggesting gene coregulation (Additional file 4).

Reciprocal transcriptome responses to SA in cassava genotypes
The temporal patterns of gDEG expression in ECU72 and COL2246 following SA, JA, ET, and ABA treatments were visualized using heatmaps (Fig. 1d).Strikingly, most SA-or JA-responsive gDEGs were reciprocally regulated.In contrast, most ET-and ABA-responsive gDEGs had similar expression trends in ECU72 and COL2246 (Fig. 1d).To uncover the biological processes associated with hormone responses between the genotypes, k-means clustering and functional enrichment analyses of gDEGs were performed (Additional file 6: Figures S3-S6; Additional file 7: Tables S16-S19).In response to all hormone treatments, the GO term categories response to hormone, stimulus or biotic stimulus were enriched, confirming the effectiveness of hormone treatments in eliciting hormone responses (Additional file 6: Figures S3-S6).Clustering also revealed that differential response between the genotypes was attributed to reciprocal regulation for SA gDEGs, but was due to differing 0-h transcript levels for JA gDEGs (Additional file 6: Figures S3-S4).

Cassava and Arabidopsis hormone responses have diverged
To place cassava's hormone responses within a broader context, the SA-and JA-dependent transcriptomes of cassava were compared to those of Arabidopsis thaliana.Previous microarray or RNA-seq studies profiling SA and JA responses in Arabidopsis utilized various hormone concentrations and plants of different ages resulting in marked differences in DEGs identified in each study (Additional file 8: Figure S7a,b; Additional file 9: Tables S20 and S21).Therefore, RNA-seq analyses of SA-and JA-treatment time courses (0, 0.5, 1, 2, 4, 8, 12, and 24 hpt) in Arabidopsis were performed for comparison to cassava time courses (Additional file 10: Figure S8a; Additional file 11: Tables S22 and S23; Additional file 12: Table S24).Of the Arabidopsis SA and JA tDEGs we identified, over 20% were identified by previous transcriptome studies (Additional file 8: Figure S7c,d).
PCA analyses indicated that the Arabidopsis response to SA and JA was separated into two phases similar to cassava.However, Arabidopsis' response timing was more similar to COL2246, as ECU72's biphasic response to SA and JA was not observed (Fig. 1c; Additional file 10: Figure S8b).In comparing SA-, JA-, ET-, or ABA-responsive genes in Arabidopsis and cassava, strikingly different regulatory programs between the species were revealed (See figure on next page.)Fig. 1 Cassava transcriptome profiles in response to defense hormone treatments.a, b tDEG and gDEG counts in ECU72 and COL2246 during SA and JA treatments (0.5, 1, 2, 4, 8, 12, and 24 hpt).A visual definition of tDEGs and gDEGs is provided.Number of up-and down-regulated genes (red and blue, respectively) are displayed.Total number of DEGs for each treatment are provided.c PCAs of detected gene reads prior to and after SA and JA treatments (0, 0.5, 1, 2, 4, 8, 12, 24 hpt) in ECU72 and COL2246.Clustering of time points defined distinct early (E) and late (L) response phases in ECU72 versus COL2246 for which a visualization is provided.ECU72's early phase (0-1 hpt) samples were more similar to COL2246's late phase (4-12 hpt) and vice versa.Normalized read count values for three biological replicates are shown per time point.Time points and genotypes are labeled by color and shape, respectively.d Heatmaps displaying gDEGs in ECU72 and COL2246 in response to SA, JA, ET, and ABA treatments.Expression is displayed as log 2 FC values relative to 0 hpt.DEG expression values are provided in Additional file 1 and DEG counts in Additional file 2. DEGs had |log 2 FC|≥ 1 and FDR ≤ 5% Fig. 1 (See legend on previous page.)[34] (Fig. 2; Additional file 11; Additional file 13: Tables S25 and S26).Of the phytohormone-responsive genes in Arabidopsis, only a small core were responsive to the same phytohormone in both cassava genotypes.For these core genes, Arabidopsis' regulatory programs did not always align with cassava responses.For example, many hormone up-regulated genes in Arabidopsis were downregulated in ECU72 or COL2246 (Fig. 2).
Given the marked differences in the transcriptome responses of Arabidopsis and cassava to the defense hormones, we determined the expression programs of 162 Arabidopsis genes central to SA, JA, ET, and ABA biosynthesis, modification, transport, and signaling (encompassing perception, signaling and response) and their orthologs in cassava [35][36][37][38] S27-S30).Over 40% of Arabidopsis hormone-pathway genes had multiple orthologs in cassava, with many displaying diverged expression programs (Additional file 16: Table S31).
The species differences were most apparent in the SA pathway (Additional file 14: Figure S9).Many SApathway genes in Arabidopsis had regulatory trends similar to one cassava genotype and displayed reciprocal regulation between the cassava genotypes.Several of these genes (MeNPR1, MeWRKY70a-b, MeGRX480c, MeCBP60a1, and MeSARD1a) shift from up-to downregulation or vice versa from the early to late SA response phases in ECU72, suggesting they may act as transcriptional regulators for this phase shift (Additional file 14: Figure S9).Similarly, the early responses of ET-and ABA-pathway genes were substantially different between the two species.For instance, the ETinduced, single-copy gene AtPR3 had 15 orthologs in cassava, 13 of which were undetected or repressed after ET treatments [32] (Additional file 14: Figures S10 and  S11).In contrast, the JA-pathway genes of Arabidopsis, ECU72 and COL2246 had overall similar responses to JA.A notably lower or no induction of MeVSP and several MeLOX2 genes was however observed in COL2246 (Additional file 14: Figure S12).

Genotypic differences in cassava's responses to whitefly infestation
To define cassava's response to whiteflies, we analyzed the transcriptomes of whitefly-resistant ECU72 and whiteflysusceptible COL2246 at 0, 1 (adults/eggs), 7 (eggs), 14 (1 st instars), and 22 (2 nd -3 rd instars) days post infestation (dpi) with the Latin American whitefly A. socialis (Fig. 3; Additional file 17: Tables S32 and S33).The temporal expression of whitefly tDEGs in ECU72 and COL2246 were distinct.At 1 dpi, when adults are feeding and eggs are deposited, COL2246 had over fivefold more tDEGs than ECU72.However, by 14 to 22 dpi, when nymphs are feeding, the magnitude of ECU72's response exceeded that of COL2246; in particular, ECU72 repressed over twofold more genes than COL2246 at these later times (Fig. 3a,b; Additional file 18: Table S35).Analysis of gDEGs also indicated that ECU72 and COL2246 had different transcriptomes prior to and during adult feeding and oviposition (0-1 dpi), with more pronounced differences after initiation of nymph feeding 14-22 dpi (Fig. 3c; Additional file 17: Table S34; Additional file 18: Table S36).PCA analyses further revealed that while ECU72 displayed a distinct constitutive/early (0 and 1 dpi) and late (7-22 dpi) infestation response phase, COL2246 had no clear delineation of phases (Fig. 3d).
To identify biological processes associated with whitefly infestation gDEGs, k-means clustering and GO-term enrichment analyses were performed (Additional file 19: Figure S13; Additional file 20: Table S37).Of note, Cluster 2 gDEGs were enriched for GO-term categories related to cell wall processes, and displayed higher transcript levels in ECU72 versus COL2246 at 0-1 dpi followed by a steep decline by 7 dpi.All other clusters were enriched for GO-term categories related to defense (Additional file 19).

Whitefly infestation of cassava triggers phytohormone-dependent and -independent gene expression
To gain a global understanding of the association of hormone-responsive genes with whitefly infestation, we classified ECU72's and COL2246's whitefly-regulated DEGs by their responsiveness to a single or multiple defense hormones (Fig. 4a-c; Additional file 21: Tables S38 and S39).Most ECU72 and COL2246 tDEGs (29-47%) responded to signals independent of SA, JA, ET, and/or ABA.The tDEGs regulated by all four hormones (12-22%) made up the next largest hormone-response class.Many of the remaining whitefly-regulated tDEGs were responsive to SA, JA or both SA and JA (Fig. 4a,b).In addition, ABA appears to be an important regulator in ECU72's response to whitefly infestation at 14 to 22 dpi, as the number of ABA-responsive genes (in the "Other" category) was over twofold higher in ECU72 than in COL2246 (Fig. 4a,b; Additional file 21: Table S38).In contrast to tDEGs, most whitefly-responsive gDEGs were responsive to SA (5.2-14.0%)or all hormones (28.1-54.6%)at all infestation times (Fig. 4c; Additional file 21: Table S39).
The levels of SA, JA, ABA, and their derivatives during whitefly infestation (0, 0.5, 1, 7, 14, and 22 dpi) in ECU72 and COL2246 were extracted from an untargeted metabolomics data set [33] (Fig. 4d; Additional file 22: Tables S40 and S41).The JA precursor 12-OPDA (12-oxo-phytodienoic acid) was significantly higher in COL2246 at Fig. 2 Comparison of Arabidopsis and cassava hormone responses.a-d Venn diagrams and heatmaps comparing tDEGs in Arabidopsis and cassava (COL2246 and ECU72) identified during SA, JA, ET, and ABA treatments, respectively.Core genes shared between Arabidopsis and both cassava genotypes are visualized in heatmaps.Arabidopsis SA and JA tDEGs were identified in this study.Arabidopsis ET-and ABA-responsive tDEGs at 0.5, 1 and 3 hpt were retrieved from Goda et al. [34].Expression values are provided in Additional files 11 and 13 0 dpi.As OPDA regulates defense genes in a JA-independent and -dependent manner [39][40][41], 12-OPDA may mediate constitutive differences in gene expression between the genotypes.SA levels were higher at all times in COL2246 versus ECU72, reaching statistical significance at 22 dpi.Reciprocally, the levels of SAG (SA glucoside, SA's storage form) demonstrated a trend of higher levels in ECU72 from 0-7 dpi.In contrast, ABA levels were higher in ECU72 with statistical significance at 0.5, 1, and 22 dpi.The highest ABA levels were reached at 0.5 dpi, preceding transcripts changes by 12 h.In addition, the levels of the oxidation metabolite of ABA phaseic acid (PA) followed trends similar to ABA (Fig. 4d).Correlation analyses additionally showed that changes in hormone levels were associated with changes in transcript  Hormone-response categories included responses to single or multiple defense hormones.Hormone categories that contributed to more than 10% of the whitefly-infestation response are shown.The hormone-response category "All" reflects the ability of a DEG to independently respond to SA, JA, ET, and ABA.The "WF" category indicates that DEGs responded to whitefly infestation but none of the defense hormones tested.The hormone category "Other" includes all single or multiple hormone-response categories that constitute less than 10% of whitefly-responsive genes at a time point.Gene counts by hormone category are provided in Additional file 21. d Defense hormones detected during infestation in ECU72 and COL2246.JA and JA-Ile were not detected and no significant changes in MeJA levels were observed during whitefly infestation of ECU72 or COL2246.Asterisks (red) indicate significant difference in hormone level between genotypes as identified by Student's t-test (* = p ≤ 0.05).Complete lists of detected hormones and p-values is provided in Additional file 22 including SA, SAG (SA glucoside), MeJA (methyl jasmonate), 12-OPDA(12-oxo-phytodienoic acid), ABA, and PA (phaseic acid) levels in ECU72 and COL2246 during whitefly infestation (Additional file 23: Figure S14c,d).
To determine the possible involvement of hormonepathway genes in infestation, the hormone and whitefly responses of defense-hormone pathway genes that were whitefly gDEGs were visualized by heatmaps (Additional file 24: Figure S15; Additional file 15).Many genes that were whitefly-upregulated in both genotypes and more highly expressed in COL2246 at 14-22 dpi (Cluster 1) were involved in positive regulation of SA signaling (MeNPR1, MeSARD1a and MeWRKY70a-b) or SA responses (MePR1a and b).In contrast, genes that were induced in ECU72 but repressed in COL2246 during infestation (Cluster 3) were predominantly involved in ABA signaling, ET signaling or SA inactivation.These ABA signaling genes are known co-receptors involved in negative feedback loops in Arabidopsis [42], suggesting their whitefly induction may be involved in fine-tuning ECU72's ABA response.

Resistant cassava deploy multiple biochemical defenses in response to whitefly adult feeding and egg deposition
To assign biological functions to the gDEGs identified in both whitefly and hormone treatments, GO term-enrichment analyses were performed (Additional file 25: Tables S42 and S43; Additional file 26: Tables S44 and S45).gDEGs were grouped by infestation time point, hormone responses and up-or down-regulation in ECU72 versus COL2246 during whitefly infestation.Enrichment results are shown as categories in defense (Fig. 5) or other processes (Additional file 27: Figure S16).Among the defense-associated gDEGs that were up-regulated in ECU72 prior to (0 dpi) or early after infestation (1-7 dpi), 81% responded to all four hormones and were associated with GO terms including cell wall-related and immune Fig. 5 Functional enrichment of cassava gDEGs in response to whitefly and hormone treatments.GO-term enrichment was performed on whitefly gDEGs.Numbers of genes enriched for terms linked to defense are shown.gDEGs responsive to one-three hormones, all hormones, or that are hormone-nonresponsive (WF) are shown.Genes upregulated in ECU72 or COL2246 are displayed on the right and left sides of the x-axis, respectively.Counts and identities of genes within each GO term category are provided in Additional file 25: Table S42 and Additional file 26: Table S44 system processes, glucosinolate, phenylpropanoid and lignin metabolism, and hormone/stimulus responses (Fig. 5).
Visualizing the RPKM expression trends for 23 genes within these enriched GO-term categories showcased associations between infestation and hormone regulation (Additional file 28: Figure S17).Genes involved in cell wall remodeling (MeXTH23, MeGUX1 and MePME31), cell wall biosynthesis (MeCESA4), response to fungal cell wall elicitors (MeCAD8i) [43], and lignin biosynthesis (MeCOMTf, MeCCOAMTa, MeMYB63, MeLAC4, and MeLAC5) were more highly expressed in ECU72 versus COL2246 prior to infestation (0 dpi) and at most times during infestation and hormone treatments.Seven immunity genes displayed similar expression trends (Additional file 28).
As four glucosinolate metabolism genes (MeSOT17, MeCYP83B1a-c) [44] were more highly expressed in ECU72 in all treatments (Additional file 28), we speculated that glucosinolate levels may be altered in ECU72.Previous publications report the ability of exogenously expressed cassava enzymes to produce glucosinolates [45,46].Therefore, given that MeCYP83B1a-c were induced in whitefly-infested ECU72, and Arabidopsis AtCYP83A1 and AtCYP83B1 catalyze the conversion of valine, isoleucine and phenylalanine aldoximes into their corresponding glucosinolates, we measured their levels in cassava leaves [44].Desulfonated glucosinolates, the thioglucose moieties characteristic of glucosinolate structures, or molecular ions corresponding to glucosinolates derived from either valine, isoleucine or phenylalanine were not detected in cassava leaves (Additional file 29: Figure S18).Our results suggest that glucosinolates may be at undetectable levels or that cassava's MeCYP83B1 genes are acting in alterative pathways in infested leaves.
As several genes involved in lignin biosynthesis were identified as up-regulated gDEGs in the whitefly-resistant ECU72 (Additional file 28), we identified the set of cassava orthologs associated with this metabolic pathway that were whitefly gDEGs (Fig. 6; Additional file 30: Table S46).Among these ten genes, eight were more highly expressed in ECU72 during infestation (Fig. 6).

Whitefly nymph feeding induces a surge in SA and immune signaling in susceptible cassava
A major shift occurred at the time of nymph feeding (14-22 dpi) when many defense-associated GO terms were enriched among a surge in genes upregulated by infestation in COL2246 versus ECU72.Genes within these GOterm categories were mainly regulated by all hormones or were defense-hormone independent (Fig. 5).Visualization of 24 individual gene transcript levels exemplified such expression trends (Additional file 31: Figure S19).
For instance, within the immune system process category, MePR3a, MeFMO1, NLR MeRPP8, and MePERK1ac (receptors that perceive cell wall perturbations) were more highly expressed in COL2246 than in ECU72 during infestation and all hormone treatments.Four starch catabolism and three sesquiterpenoid biosynthetic genes also had higher expression in COL2246 versus ECU72 by late infestation.Although sesquiterpenoid production is known to be JA induced in plants, such as maize and rice [47,48], the cassava sesquiterpenoid biosynthetic geneswere surprisingly not regulated by SA, JA, ABA, or ET (Additional file 31).Interestingly, for genes enriched in processes other than defense (i.e., metabolism, carbohydrate metabolism, transport), a similar shift in enriched terms from 0-7 versus 14-22 dpi was observed (Additional file 27; Additional file 26: Table S45).
Given the reciprocal regulation of SA-responsive genes in ECU72 and COL2246 (Fig. 1d) and detection of many SA-pathway genes in the enriched GO-term categories (Fig. 5), we examined 23 whitefly gDEGs identified in our annotation of the cassava SA pathway and/ or identified in our enrichment analyses (Fig. 7; Additional file 15: Table S27).Of these genes, 15 were signaling genes, and all but one, were more highly expressed in COL2246 versus ECU72 at 14-22 dpi.Notably, four were SA-responsive gDEGs (MeEDS1a, MeWRKY41, MeNPR1, and MeWRKY70a) and four were gDEGs that responded to both SA and JA (MeSARD1a and b, MeNPR3, and MeWRKY70b).At one or more early infestation times (0-7dpi), MeSMTa, MeSMTb and MeSAMTa were more highly expressed in ECU72; these enzymes convert SA into the mobile signal of systemic resistance, MeSA (methyl salicylate).In contrast, two of the three MeUGT74F1 genes, which encode for enzymes that convert SA into its storage forms SAG (SA glucoside) or SGE (SA glucose ester), were more highly expressed in COL2246 at 0-7 dpi (Fig. 7).Higher expression of these SA-modifying genes may reflect an attempt of COL2246 to modulate its high levels of free SA that occur early after infestation.

Divergent hormone responses in cassava versus Arabidopsis
In comparing the hormone-responsive transcriptomes of cassava and Arabidopsis, we found significant differences in gene expression both globally and within hormone pathways (Fig. 2; Additional file 14).Additionally, the numbers of hormone-pathway genes present were often expanded or contracted between species (Additional file 16).As neofunctionalization is often observed in polyploid species [49], such expansions may suggest adapted functions in defense.For example, the PR3 chitinase gene is expanded in cassava, tomato and Eucalyptus relative to Arabidopsis [32,50,51], suggesting an adaptive response to insect/pathogen-derived chitin [52].Finally, we found that many canonical markers of SA, JA, ET, or ABA signaling in Arabidopsis are not sentinels of phytohormone pathways in cassava, similar to some tomato PR genes [53].This includes the classical Arabidopsis PR genes [32], as well as the JA-induced Arabidopsis marker AtVSP1.

Hormone responses of whitefly-resistant and -susceptible cassava
Our comparison of ECU72's and COL2246's hormone levels and hormone-responsive transcriptomes revealed several key findings.While COL2246 responded similarly to all hormones, ECU72 showed evidence of antagonism between its SA and JA versus ET and ABA responses (Additional file 4).Most striking, while ECU72 and COL2246 displayed similar temporal responses to JA, ET and ABA, they showed largely reciprocal regulation by SA.In cacao, a similar reciprocity of some SA-responsive genes was seen in fungus-resistant versus -susceptible genotypes [54].Also similar to our findings, van Leeuwen et al. [55] reported that Arabidopsis ecotypes showed variable responses to SA and Proietti et al. [56] showed variation in the degree of SA antagonism with JA in Arabidopsis accessions.These findings suggest hormone responses may vary considerably between and even within plant species.
Our finding that ECU72 displays a faster-acting and more prolonged SA response than COL2246 also suggests that a master transcriptional regulator may facilitate this phase shift.The SA-signaling genes MeNPR1, MeWRKY70a-b, MeGRX480c, MeCBP60a1, and MeSARD1a are possible candidates as their expression profiles change during this early to late phase shift in ECU72 (Additional file 14: Figure S9).In Arabidopsis, AtNPR1, AtWRKY70 and AtGRX480 promote SA   S27 responses important for resistance to biotrophic pathogens, while suppressing JA/ET responses important for resistance to insect herbivores (like whiteflies and caterpillars) and necrotrophic pathogens [28,[57][58][59].AtWRKY70 is also known to suppress ABA-induced stomatal closure [60].As MeNPR1, MeWRKY70 and MeSARD1a were more highly induced at 14 dpi in COL2246 (Fig. 7), it is tempting to speculate that one or more of these transcription factors enact an ineffectual SA response in response to whitefly infestation in COL2246.Conversely, their fine-tuned regulation in ECU72 may allow for an ABA-mediated resistance response to whitefly in ECU72.

Resistance to whitefly adults and eggs is associated with ABA and cell wall defenses
Our integration of whitefly and hormone transcriptome and metabolite data sets highlighted the association of ABA and SA with whitefly resistance and susceptibility in ECU72 and COL2246, respectively (Fig. 4; Additional file 24).The observed association of high ABA and low SA levels in ECU72 could possibly have been achieved through SA-ABA antagonism, something that has been previously described in Arabidopsis-pathogen interactions [16].Furthermore, ECU72's lower levels of free SA and higher levels of SAG (SA's vacuolar storage form of SA) at 0.5-7 dpi [61] also suggests SA modification mechanisms may be important for quelling ineffective SA responses in ECU72 during early infestation.
The role of ABA in ECU72's fast-acting resistance phenotype is intriguing as increases in ABA is an early response to PAMPs that promotes stomatal closure to interfere with pathogen access to a leaf 's interior spaces [62].With limited evidence for hemipteran-plant interactions, in whitefly-resistant ECU72, ABA-mediated stomatal closure could slow phloem feeding by decreasing transpiration rate, or, it could decrease the release of whitefly-attracting plant volatiles.Alternatively, ABA may promote stomatal opening to improve gas exchange for resumption of basal photosynthesis rates after initial responses to whitefly attack [63].Several studies have shown either ABA or osmotic stress responses to be important for plant resistance to whiteflies [24,29,30] and other insect pests [64], however the overall role of ABA in defense against insects remains unclear.
Enrichment analyses additionally showed that during early infestation (0-7 dpi), few defenses were mounted in COL2246.In contrast, cell wall defenses like lignin biosynthesis were already active in ECU72 (Figs. 5 and  6).The importance of cell-wall based defenses against whiteflies discovered here is corroborated by our findings of higher basal and whitefly-induced leaf lignin levels in ECU72 versus COL2246 by Perez-Fons et al. [33], as well as observed cell wall-based defenses in cotton responses to whiteflies [25,65,66].Our analyses additionally pointed to the importance of cell wall fortification and the sensing of cell wall damage during ECU72's early infestation response (Additional file 28).Strengthening the cell wall by enhancing lignin biosynthesis to deter the probing of whitefly stylets and the perception of damage through cell wall elicitors/DAMPs may be important components of ECU72's defense during early stages of infestation [67][68][69].It is also interesting to note that by late infestation (14-22 dpi) even though fewer nymphs are present on ECU72, it has a higher magnitude response than COL2246.This suggests that ECU72 is able to prompt prolonged whitefly defenses even with minimal insect signals.

Susceptibility to whitefly nymphs is associated with a rise in SA and immune signaling
A marked shift in the cassava transcriptome was observed at the times of voracious phloem feeding by first to third instar nymphs (14-22 dpi) when a surge in SA and immune responses in COL2246 was initiated (Fig. 5; Additional files 24 and 31).Per our direct observations, at this time, the susceptible plant has a heavy load of developing 2 nd -and 3 rd -instar nymphs.In contrast, on the resistant plant, the majority of nymphs have ceased development or perished thereby reducing the quantity of effectors/elicitors being delivered to the resistant plant; this may explain the absence of a burst in SA responses in ECU72, as was observed in COL2246.Indeed, it is known that an insect's developmental stage influences plant responses at the transcriptome and metabolome level [32,33,70].The role of SA in plant-whitefly interactions has been shown to vary based on the plant and pest species [71].However, like COL2246, SA responses mounted during whitefly infestation do not interfere with but they promote whitefly performance on whitefly-susceptible Arabidopsis, tobacco and lima bean [26,28,72].
Starch catabolism and sesquiterpenoid biosynthetic genes were additionally more highly expressed in COL2246 at 14-22 dpi (Additional file 31).As observed in other plants, the breakdown of starch may be a strategy to mobilize stored energy to compensate for photosynthate depletion due to insect infestation [73].Emission of terpenoid volatiles has also been previously found in response to aphids, whiteflies and other insects for varied purposes including attracting or deterring insects or their natural predators or signaling responses in nearby plants [63,[74][75][76].
Lastly, it is important to note that a large portion of whitefly-responsive genes were responsive to all or none of the hormones tested indicating that they may be controlled by other defense signals.One possibility is that reactive species (ROS) underlie the different SA levels and responses seen in ECU72 and COL2246, as ROS can regulate SA pathway crosstalk via redoxregulated transcription factors like GRX480 and NPR1 [77][78][79].

Conclusions
Here, we provide the first global analysis of cassava's response to whitefly infestation, as well as to the defense hormones SA, JA, ET, and ABA, in the whitefly-resistant and whitefly-susceptible cassava genotypes ECU72 and COL2246.Comparisons of hormone-dependent transcriptomes in cassava and the model plant Arabidopsis revealed a striking divergence in expression programs.We suggest that such interspecies divergence may be more common than currently understood, and caution the assumptive adoption of commonly used Arabidopsis defense sentinel genes in other plant species.
Comparisons of hormone-responsive transcriptomes and findings in the Arabidopsis literature suggested several possible SA-signaling genes that may facilitate the observed SA-ABA antagonism in ECU72, ECU72's fastacting SA response involving a phase shift, and the resulting largely reciprocal response of ECU72 and COL2246 to SA treatment.However, additional genetic testing in gene-edited or transgenic cassava is necessary to determine a possible role of these genes in defense against whiteflies.Such observed intraspecies variation in phytohormone responses, while understudied, may reflect fine tuning to the suites of pests, pathogens or conditions to which each genotype has adapted.
Our integrative analyses together suggest that COL2246's late infestation response to large numbers of actively feeding nymphs occurs due to the absence of effective early control strategies and ineffective SAmediated defenses.In contrast, ECU72's faster response to whitefly adults and eggs via ABA-mediated responses and lignin-based cell wall defenses may underlie its resistance to whitefly infestation.A possible link between ABA and stomatal responses to whitefly in ECU72 remains unclear and requires further investigation.Additionally, as many whitefly-regulated genes were unresponsive to tested hormones, further studies are required to identify remaining unknown signals important for regulating cassava's response to whitefly infestation.
Hormone-responsive transcriptomes generated in this study will serve as a valuable resource to the cassava defense community.This genetic material could also be used for the construction of introgression populations, elucidating QTLs with the potential to contribute to the reduction in yield-loss resulting from whitefly and other vector-borne disease states.

Plant growth, insect rearing, and infestation experiments
In vitro-grown cassava (M.esculenta) genotypes ECU72 and COL2246 [80,81] from the CIAT collection were grown for 3 months before use in whitefly-infestation and hormone-treatment experiments as described in Irigoyen et al. [32].The Aleurotrachelus socialis Bondar colony used for cassava infestation experiments were maintained at CIAT as described by Bellotti et al. [8].
Arabidopsis thaliana Col-0 seeds (sterilized with chlorine gas and cold-treated for 2 d) were sown on ½ MS 1% sucrose agar plates [82] and kept at room temperature under constant light.One week after plating, seedlings were moved to soil (autoclaved Sunshine Mix (Sun Gro Horticulture, Agawam, MA)) supplemented with 2% Osmocote (w/w) (The Scotts Company, Marysville, OH).Plants were grown in a growth chamber under incandescent and fluorescent lights (180 μE m −2 s −1 ) under a short-day light cycle (6-h light/18-h dark) at 24 °C for 27 days, then adjusted to a long-day light cycle (16-h light/8-h dark) for one day before use in hormone-treatment experiments.
All of the plant propagation and experiments performed with plants, hormones and/or insects complied with relevant institutional, national, and international guidelines and legislation.

RNA extraction, cDNA library preparation, sequencing, and data processing
Cassava RNA extraction was performed as described by Behnam et al. [83].For Arabidopsis leaves, finely-ground frozen tissue was used for RNA extraction, which followed a standard phenol-chloroform phase separation via centrifugation and LiCl precipitation of RNAs [53].RNA quality assessment, strand-specific cDNA library preparation and RNA-sequencing for cassava and Arabidopsis samples were performed as according to Irigoyen et al. [32].Libraries were prepared for the three biological replicates of each time point.For libraries from whiteflyinfestation experiments, the Illumina NextSeq500 and Illumina HiSeq2500 platforms were used to sequence single-end 75-bp reads (trimmed to 50-bp) and 50-bp reads, respectively.For libraries from Arabidopsis and cassava hormone-treatment experiments, the Illumina Next-Seq500 platform was used to sequence single-end 75-bp reads.Twelve to fifteen libraries were multiplexed per lane.An average of ~ 30-50, ~ 13-52 and ~ 20-38 million reads among the three biological replicates were obtained with reproducibility confirmed by Pearson correlation coefficient values of 0.85-1.00,0.73-0.99 and 0.69-0.98 for Arabidopsis hormone treatments, cassava hormone treatments and cassava whitefly infestations, respectively (Additional file 32: Figures S20-S23).Sequencing and read adaptor removal was carried out at the UCR Institute for Integrative Genome Biology Genomics Core.
Following read adaptor removal, quality filtering was performed.Reads with a Phred quality score below 20 for more than 10 bases were removed.For cassava samples, reads were aligned to the reference Manihot esculenta genome version 6.1 with associated annotation information available through Phytozome [84].For Arabidopsis samples, reads were aligned to the reference Arabidopsis thaliana genome version 10 with associated annotation information available through TAIR [85].Read alignment was performed using HISAT2 [86] and DEG calling using DESeq2 [87].Read filtering, alignment, and DEG calling were performed using systemPipeR [88].Detected genes were defined as genes with an average of ≥ 20 reads across a genotype's treatment time course.Temporal and genotype DEGs (tDEGs and gDEGs, respectively) were identified by comparisons of 0.5-24 hpt to 0 hpt and transcript levels in ECU72 versus COL2246 during treatment, respectively.DEGs were defined as having |log 2 fold change (FC)|> 1 and false discovery rate (FDR) ≤ 5%.Expression values obtained by RNA-seq were validated by qRT-PCR using seven tDEGs including PR and hormone biosynthetic genes (Additional file 33: Figure S24; Additional file 34: Table S47) as in Irigoyen et al. [32].

Ortholog identification
Hormone-and lignin-pathway genes (Additional files 15 and 30) in Arabidopsis thaliana (TAIR ver.10) [85] were assigned orthologs in the M. esculenta genome version 6.1 obtained from Phytozome (JGI) [84] using the online program eggNOG-mapper version 4.5 [89,90].Protein sequence queries were matched with the most similar sequence (termed a "seed ortholog", supported by an e-value ≤ 0.001 and a score ≥ 60) within the eggNOG database.eggNOG orthologous groups (OGs) based on precomputed phylogenies were also assigned.The locations of cassava seed orthologs relative to Arabidopsis genes within an OG phylogenetic tree was used to inform cassava ortholog nomenclature.In the case of Arabidopsis genes equally distant from one or more cassava orthologs, the Arabidopsis gene with the lowest numeral was used for naming (Additional files 15 and 30).Plantspecific OGs (virNOGs) were used in naming cassava orthologs over database-wide OGs (NOGs).All other cassava gene names were obtained from Arabidopsis orthologs annotated for the M. esculenta genome (Phytozome ver.6.1) [84].

Gene expression clustering and visualization
Ordering of gene expression values by k-means or hierarchical clustering was performed using the base R stats package v3.6.2 and the R package ComplexHeatmap [94], respectively.Heatmaps displaying gene expression values (mean RPKM or log 2 FC) were constructed using the R package ComplexHeatmap.Boxplot or line graphs were constructed using the R package ggplot2.Boxplot whiskers represent values within 1.5 × IQR, and box values represent the first quartile, median, and third quartile values.Arabidopsis ET and ABA microarray data sets (0.5, 1, and 3 hpt) from Goda et al. [34] were compiled for comparison of Arabidopsis and cassava hormone-responsive DEG counts.The R package VennDiagram [95] was used to construct Venn diagrams.Gene expression data is available in Additional files 1, 11 and 17.

Metabolite quantification and transcript correlation
Relative levels of hormones were obtained from an untargeted metabolomics data set previously prepared, analyzed and published by Perez-Fons et al. [33] and Drapal et [96].The same cassava samples were used in our study.Significant differences (p ≤ 0.05) between genotypes ECU72 and COL2246 at each time point for each phytohormone detected were calculated by Student's T-test (Additional file 22: Table S41).
The methods for detection and characterization of putative glucosinolates were established with a pool containing equal parts of infested ECU72 and COL2246 samples using the collection methods of Perez-Fons et al. [33] and the standard methodologies of Crocoll et al. [97] and Clark [98].Dried leaf powder (30 mg) was dissolved in 1 mL of 85% methanol and shaken for 4 min at room temperature, centrifuged at 20,000 g for 5 min and the pellet was discarded.An aliquot of 100 μL was used as crude extract and the remaining solution used for in-column sulfatase treatment.Sulfatase (Sigma-Aldrich) solution and DEAE-Sephadex A25 column (Sigma-Aldrich) were prepared as detailed in Crocoll et al. [97].Briefly, a 200-μL column solution in 20 mM potassium acetate (pH 5) was loaded into 1-mL glass pipette tip with glass wool and remaining crude extract (900 μL) was subsequently loaded onto the column.The column was washed twice with 70% ethanol and water (100 μL each).At this time, 20 μL of the sulfatase solution was added and the reaction incubated overnight at room temperature.Reaction products (desulfonated glucosinolates) were collected by eluting with 200 μL of water and stored at -20° C until LC-MS analysis.The procedure was repeated using a solution of sinigrin and glucotropaeolin (0.5 mg/mL, 200 μL) as references for aliphatic and phenyl derivatives of glucosinolates, respectively, as positive control, and a blank solution (70% ethanol) served as a negative control.In a separate experiment, the samples were prepared by extracting plant material at 100° C for 4 min to inactivate any residual myrosinase activity.Both experiments were compared and no difference in metabolite composition was observed.
For the analysis of the crude extracts and collected fractions, an Agilent's 1290 UPLC and a 6560 Ion mobility Q-TOF mass spectrometer equipped with an Agilent Jet Stream (AJS) electrospray source was used in negative mode as detailed in Drapal et al. [99].Compounds were separated in a Zorbax RRHD Eclipse Plus C18 2.1 × 50 mm, 1.8 μm using a two solvents gradient consisting of (A) 2.5% acetonitrile in water and (B) acetonitrile, both solvents containing 0.03% vol.formic acid.Gradient started at 2% B for 1 min, increase to 30% B over 5 min, stay isocratic for 1 min followed by an increase to 90% B in two minutes and stay isocratic for another two minutes.Initial conditions were restored and reequilibration lasted 3 min.The total runtime per sample was 15 min and flowrate was set at 0.3 mL/min.Nebulizer and sheath gas temperatures were 325° and 275° C, respectively; flowrate of drying and sheath gas (nitrogen) were 5 and 12 L/min respectively and nebulizer pressure 35 psi.Capillary VCap, nozzle and fragmentor voltages were set at 4000, 500 and 400 V respectively.A reference mass solution was continuously infused to ensure mass accuracy calibration at 24-25 K resolution.Injection volume was 1 μL.A mix solution of commercial standards of glucosinolates (sinigrin and glucotropaeolin) were used as reference material and for methodology validation purposes (Additional file 29: Figure S18).

GO term enrichment
The R package ClusterProfiler [103] was used to perform GO-term enrichment analyses of cassava whitefly-and/ or hormone-responsive gene sets.Results of enrichment analyses display only those terms within the "biological process (BP)" GO category.Terms were additionally grouped into categories based on shared ancestral GO terms.Significant GO terms are defined as those with adjusted p-values ≤ 0.05 and q-values ≤ 0.05.

Fig. 3
Fig. 3 Transcriptome profiles and PCA of cassava's whitefly infestation response.a-c Numbers of tDEGs in ECU72 and COL2246 and gDEGs during whitefly infestation.Number of up-and down-regulated genes (red and blue, respectively) are displayed and total number of DEGs for each treatment are provided.DEG expression values are provided in Additional file 17 and DEG counts in Additional file 18. DEGs had |log 2 FC|≥ 1 and FDR ≤ 5%.d PCA analyses of detected genes prior to and during whitefly infestation (1, 7, 14 and 22 dpi) in ECU72 and COL2246.More variation among biological replicates was observed in COL2246 than in ECU72.Normalized RNA-seq read count values for three biological replicates are shown per time point.Time points and genotypes are labeled by color and shape, respectively

Fig. 4
Fig. 4 Hormone regulation of whitefly-regulated DEGs and hormone levels in whitefly-infested cassava.a-c Stacked bar graphs displaying number of whitefly-regulated tDEGs in ECU72 (a), tDEGs in COL2246 (b) and gDEGs (c) categorized by their hormone-response class.Hormone-response categories included responses to single or multiple defense hormones.Hormone categories that contributed to more than 10% of the whitefly-infestation response are shown.The hormone-response category "All" reflects the ability of a DEG to independently respond to SA, JA, ET, and ABA.The "WF" category indicates that DEGs responded to whitefly infestation but none of the defense hormones tested.The hormone category "Other" includes all single or multiple hormone-response categories that constitute less than 10% of whitefly-responsive genes at a time point.Gene counts by hormone category are provided in Additional file 21. d Defense hormones detected during infestation in ECU72 and COL2246.JA and JA-Ile were not detected and no significant changes in MeJA levels were observed during whitefly infestation of ECU72 or COL2246.Asterisks (red) indicate significant difference in hormone level between genotypes as identified by Student's t-test (* = p ≤ 0.05).Complete lists of detected hormones and p-values is provided in Additional file 22 including SA, SAG (SA glucoside), MeJA (methyl jasmonate), 12-OPDA(12-oxo-phytodienoic acid), ABA, and PA (phaseic acid)

Fig. 6
Fig. 6 Cassava's lignin biosynthetic pathway.Expression of whitefly-infestation gDEGs during whitefly and SA, JA, ABA, and ET treatments are shown beside lignin-pathway genes.Of note, MeCOMTf and MeCCOAMTa showed strong upregulation in ECU72 versus COL2246 at most infestation and hormone-treatment time points.Expression is displayed as log 2 FC values in ECU72 as compared to COL2246 with ECU72 upregulated (red) and downregulated (blue) gDEGs displayed.Time points where gene expression was not significantly different between the genotypes are shown in white.Time points from left to right during infestation are 0, 1, 7, 14, and 22 dpi and during hormone treatments are 0.5, 1, 2, 4, 8, 12, and 24 hpt.Gene loci are provided in Additional file 30

Fig. 7
Fig. 7 Cassava's modification, signaling and transport pathway.Expression of whitefly-infestation gDEGs during whitefly and SA, JA, ABA and ET treatments are shown beside SA-pathway genes.Expression is displayed as log 2 FC values in ECU72 as compared to COL2246 with ECU72 upregulated (red) and downregulated (blue) gDEGs displayed.Time points where gene expression was not significantly different between the genotypes are shown in white.Time points from left to right during infestation are 0, 1, 7, 14, and 22 dpi and during hormone treatments are 0.5, 1, 2, 4, 8, 12, and 24 hpt.Gene loci are provided in Additional file 15: TableS27