Evolved Transcriptional Responses and Their Trade-Offs after Long-Term Adaptation of Bemisia tabaci to a Marginally Suitable Host.

Although generalist insect herbivores can migrate and rapidly adapt to a broad range of host plants, they can face significant difficulties when accidentally migrating to novel and marginally suitable hosts. What happens, both in performance and gene expression regulation, if these marginally suitable hosts must be used for multiple generations before migration to a suitable host can take place, largely remains unknown. In this study, we established multigenerational colonies of the whitefly Bemisia tabaci, a generalist phloem-feeding species, adapted to a marginally suitable host (habanero pepper) or an optimal host (cotton). We used reciprocal host tests to estimate the differences in performance of the populations on both hosts under optimal (30°C) and mild-stressful (24°C) temperature conditions, and documented the associated transcriptomic changes. The habanero pepper-adapted population greatly improved its performance on habanero pepper but did not reach its performance level on cotton, the original host. It also showed reduced performance on cotton, relative to the nonadapted population, and an antagonistic effect of the lower-temperature stressor. The transcriptomic data revealed that most of the expression changes, associated with long-term adaptation to habanero pepper, can be categorized as "evolved" with no initial plastic response. Three molecular functions dominated: enhanced formation of cuticle structural constituents, enhanced activity of oxidation-reduction processes involved in neutralization of phytotoxins and reduced production of proteins from the cathepsin B family. Taken together, these findings indicate that generalist insects can adapt to novel host plants by modifying the expression of a relatively small set of specific molecular functions.

We target a relatively unstudied topic in generalist insects’ genome–environment interactions: the differences in performance and gene expression responses when experiencing a novel host for the first time or for long periods of time. The improved performance after multiple generations on the novel host was associated with expression changes in a relatively small set of genes that differed from the ones differentially expressed during the first exposure. Genes involved in cuticle formation and oxidation-reduction processes were upregulated, while genes coding for cathepsin B effectors were downregulated. Our study brings new insights into the way genomes of generalist insects function demonstrating that transcriptional changes in a limited set of molecular functions are sufficient for achieving high performance on novel hosts.

Introduction

The ability of organisms to respond to a changing environment is critical for their success, especially in novel or adverse habitats (Scheiner 1993; Via et al. 1995; Price et al. 2003; Fordyce 2006). It usually depends to some extent on the phenotypic plasticity of the organisms, the potential of their genotypes to produce different phenotypes in response to distinct environmental conditions (Gibert 2017). Phenotypic plasticity might be adaptive, producing individuals that are better suited to local conditions, especially if the plasticity was tested before in the organism’s “evolutionary history” and experienced natural selection (Zabinsky et al. 2019; Birnbaum & Abbot 2020). In other cases, phenotypic plasticity might be neutral or nonadaptive, especially if it shifts the phenotype further away from the local optimum (Schneider et al. 2014; Huang & Agrawal 2016). Moreover, in many habitats, stressful conditions might result from a combination of multiple biotic and/or abiotic stressors (Enders et al. 2015). When facing multiple stressors, the plastic responses to each of them may act synergistically and accelerate adaptation (Ragland & Kingsolver 2007; Snoeck et al. 2018), antagonistically, slowing or even abolishing the process (Alzate et al. 2017) or may not interact in a meaningful way to genetic or phenotypic variation (Enders et al. 2015; Jackson et al. 2016).

From the molecular perspective, phenotypic plasticity is believed to involve two main mechanisms, “gene regulation” and “allelic sensitivity”. In the first, different stimuli from the environment, such as changes in diet, temperature, illumination, or different forms of stress (Mukherjee et al. 2015), cause regulatory reprogramming, which leads to the activation and/or suppression of target genes (Gibert 2017). For example, in Drosophila melanogaster, female abdominal pigmentation is a plastic trait, as it is darker in females grown at 18°C than at 29°C (Gibert et al. 2000). This plasticity is thought to be adaptive as it increases the body temperature (Gibert et al. 2000). At the regulatory level, the temperature modulates the expression of the pigmentation gene tan, which encodes a hydrolase implicated in the production of melanin (True et al. 2005). “Allelic sensitivity”, on the other hand, involves gene variants that show different activity in specific habitats. For example, a Leu-to-Phe substitution present in a gene coding for aldehyde dehydrogenase in D. melanogaster was found to increase the turnover rate of acetaldehyde but to decrease the turnover rate of larger aldehydes (Chakraborty & Fry 2016). The Phe allele variant was found to be present at variable frequencies (5–30%) in temperate populations but absent or rare (<5%) in tropical populations (Fry et al. 2008).

It is clear that both aforementioned molecular mechanisms underlying phenotypic plasticity depend on the genetic pool of the population. Populations living in mild and nonstressful environments are likely to harbor multiple alleles in loci involved in their adaptation to these environments, as they produce selectively equivalent phenotypes (Schlichting & Wund 2014; Schneider & Meyer 2017). However, stressful environments can unmask vulnerable polymorphisms and initiate a selection process, which may leave in the population only alleles that provide the highest fitness values (Sollars et al. 2003; Jarosz & Lindquist 2010; Rohner et al. 2013). At the same time, an additional strong selection process might take place on the mode of expression of the genes that showed transcriptional plasticity upon the first exposure to the stressful environment (Nylin & Janz 2009; Schlichting & Wund 2014). In this case, plastic traits may “evolve” by acquiring quantitative genetic changes that can either increase or decrease their environmental responsiveness (Levis & Pfennig 2016). In some cases, the plastic traits may even lose completely their environmental responsiveness, resulting in constitutive expression (Schneider & Meyer 2017). These changes allow the refinement of the expression levels (via selection), for obtaining optimal performance that might outcompete that obtained by the ancestral plastic trait (Pfennig & Ehrenreich 2014; Levis & Pfennig 2016; Wang et al. 2017).

Populations of generalist insect herbivores provide an ideal model for studying the role of phenotypic plasticity in adaptation to changing environments (Birnbaum & Abbot 2020). These populations live in many cases in mild heterogeneous environments containing multiple suitable plant hosts. Occasionally, however, they can accidently migrate to more homogenous stressful environments containing a combination of abiotic and biotic stressors such as dryness, heat and the presence of a marginally suitable novel host (Kaunisto et al. 2016; Gutiérrez 2020). Yet, as these are generalist insect populations, this will likely not be the “end of the story”, meaning that we should also take into consideration the backward pathways in which the populations return to their mild non-stressful habitats after experiencing stressful episodes. Therefore, generalist insect herbivores offer dynamic study systems that facilitate the comparisons of closely-related populations or species experiencing recurrent exposures to complex environments (Birnbaum & Abbot 2020). This is expected to allow accurate examination of the broad impact of plastic traits on short-term (within a generation) performance and long-term (multiple generations) adaptive responses (Müller et al. 2017; Vandenhole et al. 2021).

For example, multiple studies have indicated that the immediate stress of switching to feed on a marginally suitable novel host plant can be mitigated in generalist insect herbivores by a plastic change in the expression of genes that code for detoxification enzymes (Matzkin 2012; de La Paz Celorio-Mancera et al. 2013; Yu et al. 2016). However, comparisons of transcriptional responses across species or populations that are subjected to different environments for longer periods, clearly indicate the presence of “evolved” constitutive and plastic expression differences that significantly differ from the original short-term plastic responses (de La Paz Celorio-Mancera et al. 2013; Vogel et al. 2014; Ragland et al. 2015). On the other hand, despite a long history of studies on plant–insect interactions, the interactive effects of host plant quality and an abiotic stressor such as temperature are as yet poorly understood (Kaunisto et al. 2016). In general, several studies have shown that host plant or diet quality can alter the thermal reaction norms for key life-history traits of generalist insect herbivores indicating a strong temperature × diet interaction (Lee & Roh 2010; Clissold & Simpson 2015; Jang et al. 2015). More importantly, differences in responses to mild thermal stressors are expected when feeding on low-quality or high-quality hosts. For example, over a temperature range of 20°–30°C, survival of Manduca sexta on tobacco (a high-quality host) was uniformly high at all rearing temperatures. In contrast, survival on the low-quality host Devil’s claw (Proboscidea louisianica), decreased with decreasing temperature, suggesting that lower temperatures were stressful for individuals reared on the low-quality host plant but not for those on the high-quality host plant. Moreover, at higher temperatures, the differences in growth and development of M. sexta on low- and high-quality host plants almost disappeared (Diamond & Kingsolver 2010). So far, however, the identification of the mechanisms underlying synergistic or antagonistic stressor interactions in insects have been most commonly reported for chemical–temperature and chemical–pathogen pairs and the effects of other stressor pairs (temperature × food-stress for example) have hardly been studied (Kaunisto et al. 2016).

We focused here on the generalist whitefly Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae), a phloem-feeding insect known as one of the most destructive pests of open-field, greenhouse crops and horticulture plants worldwide (Dinsdale et al. 2010; De Barro et al. 2011). Bemisia tabaci is a cryptic-species complex widely distributed throughout tropical and subtropical regions (Oliveira et al. 2001). Species in the complex are morphologically indistinguishable, but differ significantly in many biological traits including their host range (Malka et al. 2018). The most dominant species in the complex are Middle East-Asia Minor 1 (MEAM1) and Mediterranean (MED). These species are highly polyphagous and invasive as they successfully colonized large areas all over the globe (De Barro et al. 2011). Previous studies have suggested that MEAM1 and MED harbor sufficient phenotypic plasticity to be able to adapt to new marginally suitable host plants after generations of rearing on an unrelated host. For example, the ability of MEAM1 to adapt to cassava (a marginally suitable Euphorbiaceae host for this species) was investigated by passing MEAM1 through a series of intermediate hosts. Within a short period of 10 generations, MEAM1 was capable of producing eight generations/year on cassava plants (Carabali et al. 2005). Similarly, the MED species was capable of significantly improving its performance (survival and fecundity) on tobacco plant (a marginally suitable Solanaceae host) in a short period of 10 generations (Xia et al. 2017).

Our main goal in this study was to understand the molecular process that MEAM1 populations undergo during their long-term adaptation to a well-defended and marginally suitable host plant followed by their return to a benign suitable host. For this, we established multigenerational colonies that allowed us to explore the differences, both in the performance (survival and development rate) and transcriptomic levels, between MEAM1 populations that were subjected to cotton (a suitable host) and habanero pepper (a marginally suitable host). We used reciprocal host tests to document the differences in performance and to underline the “plastic” and “evolved” differences in transcriptomic-wide gene expression. Moreover, in order to test for possible synergistic or antagonistic interactions between the plant-host stressor and a common abiotic stressor, the comparative performance assays were conducted in two temperatures, 30°C, which is nearly-optimal, and 24°C which extracts a mild abiotic cold-stress (Butler et al. 1983; Nava-Camberos et al. 2001; Zidon et al. 2016).

Results

General Insights from the Performance Assays

Performance evaluation experiments were conducted 17 generations after the cotton and habanero pepper populations were established (see more details in the “Material and Method” section below). We used four reciprocal host treatments to estimate the differences in performance of the populations (fig. 1): suitable host to suitable host (cotton to cotton, C→C), marginally suitable host to marginally suitable host (habanero pepper to habanero pepper, hP→hP), suitable host to marginally suitable host (cotton to habanero pepper, C→hP) and marginally suitable host to suitable host plants (habanero pepper to cotton or hP→C). The performance assays focused on answering three main questions: 1) did adaptation to habanero pepper occur? 2) is the adaptation to habanero pepper complete? and 3) is there a price “paid” during the adaptation to habanero pepper that reduces the performance on cotton?

Fig. 1.

Experimental set-up. First, performance evaluation experiments were conducted using four treatments: insects from the suitable host were tested on suitable host plants (cotton to cotton or C→C), insects from the marginally suitable host were tested on marginally suitable host plants (habanero pepper to habanero pepper or hP→hP), insects from the suitable host were tested on marginally suitable host plants (cotton to habanero pepper or C→hP) and insects from the marginally suitable host were tested on suitable host plants (habanero pepper to cotton or hP→C). Then, females were collected and pooled by treatment (50 individuals in each biological replicate), RNA was extracted and sequencing data were analyzed using DESeq2 and topGO R packages.

Egg-to-adult survival (fig. 2) was significantly different between the four treatments (X2: 97.239, P < 0.0001). The temperature effect by itself was not significant (X2: 0.861, P = 0.3535) while the interaction (treatment × temperature) was significant (X2: 9.110, P = 0.0279). Prominent and significant higher survival on habanero pepper was observed in the hP→hP treatment relative to the C→hP treatment, at both 24°C and 30°C (X2: 10.663, P = 0.0022 and X2: 18.413, P = 0.0003, respectively), indicating that adaptation to habanero pepper occurred in the habanero pepper colonies. Significant lower survival on habanero pepper was observed in the hP→hP treatment when compared to the survival on cotton of the C→C population, at both 24°C and 30°C (X2: 59.044, P < 0.0001 and X2: 23.210, P < 0.0001, respectively), indicating that the adaptation to habanero pepper was not complete. Comparisons between the two cotton treatments (C→C and hP→C), indicated a small but significant reduction in the ability to survive on cotton of the habanero pepper-adapted population, at both 24°C and 30°C (X2: 5.964, P = 0.0209 and X2: 10.773, P = 0.0022, respectively). Comparisons between the same treatments across the two temperatures revealed significant lower survival only in the hP→hP treatment at 24°C (X2: 7.574, P < 0.0099).

Fig. 2.

MEAM1 “egg to adult” survival by treatment. All P-values were calculated using a generalized linear model with a logit link and a binomial distribution. FDR corrections were applied to a priori paired comparisons. The four treatments were: C→C, suitable to suitable (cotton to cotton) hosts; hP→hP, marginally suitable to marginally suitable (habanero pepper to habanero pepper) hosts; C→hP, suitable to marginally suitable (cotton to habanero pepper) hosts; hP→C, marginally suitable to suitable (habanero pepper to cotton) hosts.

Similarly, significant longer development time (reduced performance) was observed in the C→hP treatment compared to the hP→hP treatment both at 24°C and 30°C (X2: 11.9484, P = 0.0005 and X2: 11.9484, P = 0.0048, respectively), bringing an additional indication that adaptation to habanero pepper occurred in the habanero pepper colonies (fig. 3). Significant longer development time, at both 24°C and 30°C (X2: 43.0663, P ≤ 0.0001 and X2: 9.3846, P = 0.0022, respectively), was observed in the hP→hP treatment compared to the C→C treatment (fig. 3), indicating again an incomplete adaptation of the habanero pepper colonies to habanero pepper (hP→hP). In addition, significant longer development of the habanero pepper-adapted population on cotton (hP→C) was observed compared to the C→C treatment, at both 24°C and 30°C (X2: 13.3405, P = 0.0003 and X2: 19.9944, P ≤ 0.0001, respectively) (fig. 3), indicating that adaptation to habanero pepper is associated not only with reduced survival but also with prolonged development on cotton, although it was the original host plant of the habanero pepper population.

Fig. 3.

MEAM1 “egg to adult” development time by treatment. All P-values were obtained using the Kruskal-Wallis Test followed by Bonferroni correction for multiple comparisons. The four treatments were: C→C, suitable to suitable (cotton to cotton) hosts; hP→hP, marginally suitable to marginally suitable (habanero pepper to habanero pepper) hosts; C→hP, suitable to marginally suitable (cotton to habanero pepper) hosts; hP→C, marginally suitable to suitable (habanero pepper to cotton) hosts.

General Insights from the Transcriptomic Analyses

As our performance assays indicated that the mild-abiotic temperature stress of 24°C enhances the impact of the stressor of primary interest, a marginally suitable host, transcriptomic analyses were conducted only on samples collected from the 24°C-treated cages. First, we used ordination techniques to analyze the expression profile of MEAM1 adult females from each of the four plant treatments. Principal components analysis (PCA) revealed that 44.26 and 24.97% of total gene expression variation across treatments could be explained by the two first components, PC1 and PC2, respectively (fig. 4). The analysis indicated that the expression profiles of insects developing on cotton were much more variable in their PC1 scores than those of insects developing on habanero pepper (F4,3: 11.351, P = 0.03714 when comparing between the hP→C and hP→hP treatments and F4,2: 30.291, P = 0.03221 when comparing between the C→C and C→hP treatments). Overall, relatively small number of genes were differentially expressed between the treatments. In the C→hP versus hP→hP comparison, there were 135 up-regulated genes and 53 down-regulated genes. In the C→C versus hP→hP comparison, there were 63 up-regulated genes and 51 down-regulated genes. In the C→C versus hP→C comparison, there were 26 up-regulated genes and 57 down-regulated genes.

Fig. 4.

Principal Component Analysis of the female RNASeq samples (24°C). Together, PC1 and PC2 explain 69.23% of the variability. The C→C and hP→C treatments had five replicas, the hP→hP treatment had four replicas, and the C→hP treatment had three replicas due to low survival rate (∼5%) on habanero pepper.

To acquire first understanding, molecular function (MF) enrichment analyses, using Gene Ontology (GO) terms were conducted, separately for the detected up-regulated and down-regulated genes. “Structural constituent of cuticle” (GO:0042302), “chitin binding proteins” (GO:0008061) and “transferase activity of acyl groups, other than amino-acyl groups” (GO:0016747), underwent the most dramatic expression changes during the adaptation to habanero pepper. These functions were up-regulated and enriched in the hP→hP treatment when compared both to the C→hP and C→C treatments (supplementary tables S1and S1, Supplementary Material online, respectively). Additional enriched and up-regulated GO term was found in the hP→hP treatment: “G protein-coupled receptor activity” (GO:0004930), which was upregulated in hP→hP when compared both to the C→hP and C→C treatments (supplementary tables S1and S1, Supplementary material online, respectively). From the other end, the most noticeable GO term that was downregulated during the adaptation to habanero pepper was the “cysteine-type peptidase activity” (GO:0008234) with 12–16 genes depending on the comparison done. This GO term was down-regulated in the hP→hP treatment when compared both to C→hP and C→C treatments (supplementary tables S1and S1, Supplementary material online, respectively). The “cysteine-type peptidase activity” GO term also showed high plasticity (both upregulated and downregulated activity) when the hP→C and C→C treatments were compared (supplementary tables S1and S1, Supplementary material online, respectively).

Specific Genes and Their Expression Patterns during Adaptation to Habanero Pepper

Genes that showed transcriptional differences between the C→C, C→hP and hP→hP treatments were divided into four patterns of expression (Ragland et al. 2015: 1) a plastic response to habanero pepper of the cotton-maintained population that did not evolve during adaptation to habanero pepper (C→C ≠ [C→ hP = hP→hP]); 2) a plastic response to habanero pepper of the cotton-maintained population that was followed by an evolved enhancing response in the same direction during adaptation to habanero pepper (C→C < C→ hP < hP→hP or C→C > C→ hP > hP→hP); 3) a plastic response to habanero pepper of the cotton-maintained population leading to different expression levels only in the transferring population ([C→C = hP→hP] ≠ C→hP); and 4) evolved differences in gene expression during the adaptation to habanero pepper, without an initial plastic response in the cotton population ([C→C = C→hP] ≠ hP→hP).

Pattern 1) harbored one up-regulated phospholipase A2 gene. This gene family was previously suggested to be involved in cellular and humoral immune responses of insects to pathogenic bacteria (Stanley-Samuelson et al. 1991; Xu et al. 2021). It also harbored eight down-regulated genes involved in carbohydrate and protein metabolism (fig. 5; supplementary table S2, Supplementary material online). Pattern 2) harbored no up-regulated genes and three down-regulated cathepsin B genes (of 49 cathepsin B genes present in the B. tabaci genome) (fig. 5; supplementary table S2, Supplementary material online). Pattern 3) harbored changes in gene expression that likely reflect a short-term (within a generation) stress response to host switching that disappear with time, which do not contribute to habanero pepper adaptation. This gene group included five up-regulated genes putatively involved in detoxification and sugar metabolism/transport, and no down-regulated genes (fig. 5; supplementary table S2, Supplementary material online). Pattern 4) harbored the largest gene group among the four expression patterns and included 49 up-regulated genes and 16 down-regulated genes (fig. 5; supplementary table S2, Supplementary material online). Up-regulated genes, likely to be associated with adaptation to habanero pepper were found to be involved in the formation of the cuticle structure (10 out of 92 genes in the B. tabaci genome that are likely to contribute to the structural integrity of the insect’s cuticle) and oxidation-reduction processes (6 genes). Other genes showed more complex expression responses with some genes with similar functions being up-regulated or down-regulated in the habanero pepper-adapted population (detoxification, lipid and xenobiotic transport, signal transduction and proteolysis). Again, reduced activity of genes belonging to the cathepsin B family was a characteristic of the colonies adapted to habanero pepper (fig. 5; supplementary table S2, Supplementary material online).

Fig. 5.

Significant differentially expressed genes (FDR ≤ 0.05 and Log2Fold ≥ +1 or Log2Fold≤ −1), displaying four patterns of “plastic” and “evolved” responses (see text for more details). Comparisons were made between colonies from the suitable host that were tested on the suitable host (C→C), colonies from the suitable host that were tested on the marginally suitable host (C→hP) and colonies from the marginally suitable host that were tested on the marginally suitable host (hP→hP). (A) a “plastic” response with no “evolved” response across environments; (B) “plastic” and “evolved” responses in the same direction; (C) “plastic” short-term responses; (D) an “evolved” non-plastic response. Purple and green represent genes up- or down-regulated relative to the C→C treatment, respectively.

Specific Genes and Their Expression Patterns during the Return of the Habanero Pepper-Adapted Population to Cotton

Genes that showed expression differences between the hP→hP, hP→C and C→C treatments were divided into three patterns of expression (Ragland et al. 2015: 1) A plastic response to cotton of the population adapted to habanero pepper that did not evolve during the 17 generations spent on habanero pepper (hP→ hP ≠ [hP→C = C→C]); 2) a plastic response to cotton of the population adapted to habanero pepper leading to differential expression levels only in the returning population ([hP→hP = C→C] ≠ hP→C); and 3) evolved non-plastic differences in gene expression during the adaptation to habanero pepper ([hP→hP = hP→C] ≠ C→C).

Pattern 1) harbored genes that were plasticly up-regulated upon returning of the habanero pepper adapted population to cotton. These genes were mostly associated with the cathepsin B (6 genes) and L families (2 genes), other proteolysis processes (5 genes) and sugar metabolism (4 genes). The genes that were plasticly down-regulated upon the return to cotton were associated with oxidation-reduction and proteolysis processes, likely indicating a return to a less stressful environment (fig. 6; supplementary table S3, Supplementary material online). Pattern 2) harbored again changes in gene expression that likely reflect short-term (within a generation) stress responses to host switching that disappear with time. These included three up-regulated genes involved in protein digestion/degradation and six down-regulated genes mostly involved in oxidation-reduction processes and transport (fig. 6; supplementary table S3, Supplementary material online). Pattern 3) harbored genes that showed higher or lower expression in the cotton-maintained population compared to the hP→hP and hP→C treatments. Genes that showed higher expression were found to be involved in a variety of metabolic functions (dominated by protein degradation, 5 of 19 genes). Genes that showed lower expression were found to be involved mainly in oxidation-reduction and lipid transport processes (fig. 6; supplementary table S3, Supplementary material online). These gene groups might be associated with the observed reduced performance (survival and development rate) of the returning habanero pepper-adapted population to cotton plants relative to the cotton-maintained population (figs. 2 and 3).

Fig. 6.

Significant differentially expressed genes (FDR ≤ 0.05 and Log2Fold ≥ +1 or Log2Fold≤ −1), displaying three patterns of “plastic” and “evolved” responses (see text for more details). Comparisons were made between colonies from the marginally suitable host that were tested on the marginally suitable host (hP→hP), colonies from the marginally suitable host that were tested on the suitable host (hP→C) and colonies from the suitable host that were tested on the suitable host (C→C). (A) a “plastic” response with no “evolved” response across environments; (B) short-term “plastic” response; (C) an “evolved” non-plastic response. Purple and green represent genes up- or down-regulated relative to the hP→hP treatment, respectively.

Discussion

In this study, we aimed to provide insights on the performance- and genome-wide transcriptional changes that occur during long-term adaptations of generalist insect herbivorous to marginally suitable plant hosts, a topic that is still only scarcely investigated. Our findings suggest that extreme generalists like the whitefly B. tabaci (Malka et al. 2018) can successfully adapt to marginally suitable host plants on which their initial survival is very low, but likely require multiple generations for the process that also involves the payment of a general fitness “price” (see below). Quite surprisingly, the adaptation was found to be driven by a relatively small set of mostly “evolved” changes in gene expression. We highlight three complementing findings, obtained in our performance bioassays.

Major Findings from the Performance Assays

First, the adaptation to habanero pepper was not complete, as both the survival and development rate indicated significant lower performance of the habanero pepper-adapted population on habanero pepper (hP→hP) compared to the performance of the cotton-maintained population on cotton (C→C). Similar findings were obtained when the Fabaceae generalist species Acyrthosiphon pisum (“G” biotype) was grown for six months on the “non-neutral” hosts Medicago truncatula, Medicago sativa and Vicia villosa (Lu et al. 2016). Although significant increase in performance was observed, the biotype did not reach its performance level on the “universal” host, Vicia faba (Lu et al. 2016). On the other hand, Hu & Tsai (2020) showed complete restoration of B. tabaci performance, 10 generations after switching a population reared on Chinese kale for five years to cotton, cucumber, poinsettia, and tomato, all suitable hosts for B. tabaci MEAM1. This highlights the possibility that the phenomenon described here is different from other reports on B. tabaci and likely occurs only when the insect experiences an extremely challenging plant, which is outside its regular host range.

Second, there were significant quantitative (development rate) and qualitative (survival rate) reductions in the performance of the habanero pepper-adapted population on cotton (hP→C) when compared to the cotton-maintained population (C→C). This might indicate some loss of genetic variation required for optimal performance on cotton. Previous studies conducted on the generalist mite Tetranychus urticae, did not obtain these findings, as mites that went through a long-term adaptation period to novel/poor-quality hosts such as tomato or cucumber did not show reduced performance when returning to their original (ancestral) hosts bean and cotton (Agrawal 2000; Wybouw et al. 2015). In an additional study, Kühnle & Müller (2011) split a population of the Brassicaceae specialized mustard leaf beetle Phaedon cochleariae, reared for more than 40 generations on Brassica rapa, into three sub-populations that continuously fed on B. rapa, Nasturtium officinale or Sinapis alba for multiple generations. Eleven generations after the split, the three sub-populations showed no differences in performance on the original host B. rapa in nearly all tested performance parameters.

Taking all this in account makes our finding here a bit surprising, as cotton plants cannot be considered as a challenging host for B. tabaci (see in “Materials and Methods”). An alternative explanation, which cannot be ruled out at this point, is the existence of nonadaptive maternal effects (Mousseau & Dingle 1991), as the “mothers” of the analyzed hP→C generation developed and emerged on habanero pepper while their offspring developed and emerged on cotton. It is important to note, however, that detailed evidence for maternal effects in whiteflies are currently lacking and previous studies on the generalist and Fabaceae generalist aphid species Myzus persicae and A. pisum, respectively, found no or very small maternal effects on offspring performance (McLean et al. 2009; Nespolo et al. 2015). Future analysis of the habanero pepper-adapted population in advanced generations will indicate if the reduction in performance on the original host is further enhanced, which might indicate that a host specialization process is taking place, limiting the potential host range of the population.

Third, we observed a significant interaction between the temperature and plant treatments (significant lower survival only of the hP→hP treatment at 24°C compared to 30°C) likely indicating an antagonistic interaction between the insect’s response to the biotic (marginal-plant host) and abiotic (mild-cool temperatures) stressors. In general, temperatures in the range of 20–24°C are expected not to be stressful for B. tabaci being reared on high-quality food resources as they were previously associated with increased egg-to-adult survival (Tsueda & Tsuchida 2011; Zidon et al. 2016). This emphasizes again the possibility that successful colonization of habanero pepper should be considered as an extreme out-of-range challenge to the insect. The antagonistic interaction between cooler temperatures and low food quality may simply result from dividing the available energy resources when responding to the two stressors (Todgham & Stillman 2013). Alternatively, it may reflect contradictions in the regulation of common defense pathways such as the production of cuticular protein (Enders et al. 2015). However, the prediction of the outcome of these mechanism-based non-additive interactions is challenging as they can easily yield both synergistic and antagonistic results (Kaunisto et al. 2016).

Global Insights from the Transcriptomic Analyses

Next, the molecular data were analyzed globally. An integrated inspection of the transcriptomic changes by PCA indicated that the pepper samples (C→hP and hP→hP) were more aggregated among themselves, while the cotton samples (C→C and hP→C) were more dispersed. One possible explanation to this finding is that only a subset combination of alleles of various loci, allow survival on habanero pepper. On the other hand, many more combinations of alleles/loci can support successful development and survival on cotton. Therefore, the set of genes/alleles that are allowed to be differentially expressed is larger, as the selective forces are weaker (Flatt 2005; Schlichting & Wund 2014; Schneider & Meyer 2017). Although the continued rearing of the colonies on habanero pepper could have led to directional selection followed by stabilizing selection in order to achieve an optimal gene expression pattern (Wagner 2000; Siegal & Bergman 2002; Flatt 2005), we believe this not to be the case as the spreading of the hP→C samples (habanero pepper-adapted colonies returning to cotton) along the PC1 axis suggests that strong selection for specific genotypes did not take place during habanero pepper adaptation. The molecular mechanism/s that allow adaptive non-variable expression pattern on habanero pepper and at the same time large diversity of gene expression upon the return to cotton (of the habanero pepper-adapted colony) remain elusive at this stage. We can only speculate that they include one or more of the following activities: specific or combined activity of molecular chaperon hubs (Hsp90 as a prominent example), regulatory microRNAs, gene-regulatory networks with epistatic interactions and epigenetic effects on gene regulation. All these regulatory mechanisms are expected to be highly affected by the polygenic nature of the adaptive trait and the possible separation of segregating loci via meiotic recombination in the “non-stressful” cotton environment (Flatt 2005; Zabinsky et al. 2019). Finally, we cannot exclude the possibility that random gene expression differences between the relatively small number of samples in the hP→C treatment contributed to some extent to the observed diversity in gene expression.

Transcriptomic Insight for Long-Term Adaptation Process to Habanero Pepper

At the functional level, our transcriptomic analysis raised the possibility that changes in the cuticle, and its different components, such as the cuticular proteins and chitin, might play a leading factor in the adaptation to habanero pepper of B. tabaci (supplementary tables S1and S1, S2, Supplementary material online; fig. 5). Previous studies allow us to suggest two possible explanations for the enrichment of up-regulated cuticle-related genes in the habanero pepper-adapted population. The first explanation is related to the possibility that cuticle barrier enhancement is a general response to multiple abiotic stressors such as heat and UV radiation (Nguyen et al. 2009; Benoit et al. 2010). Following this line, the cuticle might also play an important role in protecting the insect from desiccation, if feeding on a marginal-suitable host interferes with acquiring sufficient water during feeding. The second explanation relates to the possible hardening and thickening of the stylets. Eight of the ten up-regulated cuticle proteins in the habanero pepper-adapted population (fig. 5; supplementary table S2, Supplementary material online) have a consensus RR-2 region, suggesting their presence in rigid and sclerotized cuticles (Willis 2010; Willis et al. 2012). Mathers et al. (2017) suggested that RR-2 cuticular protein-encoding genes, enable the aphid M. persicae (which shares a similar mode of feeding with B. tabaci) to adjust the stylet to different physical and chemical attributes of the plant cell wall and defense responses. Interestingly, studies in aphids have associated the RR-2 cuticular proteins with the “acrostyle”, the tip (last few microns) of the maxillary stylets of the aphids’ mouthparts, where the food canal and salivary canal are fused (Uzest et al. 2010). This region performs intracellular punctures during probing and phloem feeding (Uzest et al. 2010). Therefore, changes in its structure can allow adaptation to the physical resistance of the host-plant cell wall in the habanero pepper-adapted population.

Another overexpressed group that showed significant enrichment in the hP→hP treatment (fig. 5; supplementary table S2, Supplementary material online) was found to be involved in oxidation-reduction processes (GO:0055114, 6 genes). One of the main responses of plants to phloem-feeding involves the accumulation of reactive oxygen species (ROS) (Moloi & Van Der Westhuizen 2006). ROS are elicitors of plant defense-signaling pathways (Fürstenberg-hägg et al. 2013) and may also have direct adverse effects on midgut tissues (Smith & Boyko 2007). The ability to induce ROS scavenging genes and repress ROS production g­enes was associated in M. persicae with adaptation to an aphid-resistant pepper accession (Capsicum baccatum, PB2013071) (Sun et al. 2020). Four of the six overexpressed genes in the oxidation-reduction-enriched group were found to belong to gene families previously associated with herbivorous-insects host adaptation: 1) one cytochrome P450 gene, coding for enzymes known to be involved in oxidation and detoxification of plant secondary compounds (Birnbaum et al. 2017; Vandenhole et al. 2021), 2) one glucose dehydrogenase gene coding for enzymes previously found to be active in the saliva of several aphid species (Carolan et al. 2009; Cooper et al. 2010), likely playing a role in suppressing plant defensive responses, and in detoxification of plant defensive compounds (Nicholson et al. 2012), and 3) two aromatic peroxygenase genes, which function in detoxification by selectively hydroxylating the aromatic ring of toxic compounds. In fungi, these enzymes have been implicated in the degradation of complex plant biomolecules (Hammel & Cullen 2008; Chen et al. 2016).

Other enriched and over-expressed groups in the hP→hP treatment include xenobiotic and detoxification functions such as “transferase activity of acyl groups, other than amino-acyl groups” (GO:0016747) and “transferase activity, transferring hexosyl groups” (GO:0016758). Genes classified under the first GO term are considered to be related to the transport of fatty acids and xenobiotics in Caenorhabditis elegans (Choy & Thomas 1999; Choy et al. 2006; Watts & Browse 2006) and were shown to be essential for cuticular modification during the development in Drosophila (Dzitoyeva et al. 2003). The second GO term was mainly UDP-glucuronosyltransferases (UDPGTs). These genes code for enzymes that are involved in phase-two detoxification through the conjugation of a diverse range of plant secondary metabolites with sugars to produce glycosides, which are water-soluble and can be efficiently excreted (Mackenzie et al. 1997). High expression of UDPGTs is likely to contribute to the specific adaptation to tobacco of the M. persicae nicotianae race of M. persicae (Pan et al. 2019) and to be involved in the ability of Aphis nerii, a polyphagous specialist that feeds on more than 50 species of milkweed and oleander plants, to handle increased levels of cardenolide toxicity (Birnbaum et al. 2017). Similarly, feeding B. tabaci MEAM1 females for 48 h on dsRNA targeting three UDPGT genes, reduced female fecundity on cabbage by ∼30%, compared to dsEGFP fed females (L. Guo et al. 2020).

Transcriptomic Insight for Habanero Pepper-Adapted Population Returning to Cotton

Our transcriptomic results also identified one down-regulated MF in the habanero pepper-adapted population (when feeding on habanero pepper plants, hP→hP). The function involves cysteine-type peptidase activity (GO:0008234), mainly the cathepsin gene families, especially cathepsin B (49 of the 78 of the cathepsin genes in the B. tabaci genome are cathepsin B). Not much is known on the role of cathepsins in sap-feeders. Most of the knowledge comes from the aphid literature. The genes were shown to be expressed in the aphids guts, likely functioning in the degradation and uptake of phloem proteins (Rauf et al. 2019). Cathepsin B enzymes might also play a complementing defensive role during insect feeding (Kehr 2006), by contradicting plant defenses in at least three different ways. First, cathepsin B enzymes contain a novel insertion loop, which is responsible both for the dipeptidyl carboxypeptidase activity of the enzymes and their lower affinity to plant-derived cysteine peptidase inhibitors (Musil et al. 1991; Terra et al. 2019). Second, plant protease inhibitors are themselves proteins, which could be inactivated by insect cathepsin B enzymes (Rispe et al. 2008). Third, cathepsin B activities have been found in the salivary glands and saliva of the southern green stink bug, Nezara viridula (Lomate & Bonning 2016), raising the possibility that some cathepsin B–like forms might be secreted by phloem-feeders in order to manipulate the plant defense reactions (Rispe et al. 2008).

The cathepsin gene family, especially cathepsin B and L coding genes, is significantly expanded in the genome of B. tabaci, relative to 15 other arthropod species (Chen et al. 2016). Over-expression of cathepsin B and L coding genes was associated with resistance to xenobiotic compounds (Xie et al. 2014), insecticide resistance and efficient virus transmission (Chen et al. 2016). In our transcriptomic data, cathepsin B-coding genes where highly expressed in the C→C treatment compared to the hP→hP treatment (fig. 6; supplementary tables S1and S3, Supplementary material online). However, many of these genes did not differ in their expression level between the hP→C and C→C treatments (fig. 5; supplementary tables S1, S1and S3, Supplementary material online), suggesting that the genes are plastic and do not lose their plasticity during adaptation to habanero pepper. Taken together, these findings suggest that the cathepsin gene family, and especially the cathepsin B group, are likely to play a role in the ability of B. tabaci to feed on cotton. The mechanism involved is not clear at this stage but might suggest a combined effect of poor nutritional content and high abundance of defensive proteins in the phloem of cotton plants. On habanero pepper on the other hand, higher performance is achieved when cathepsin B genes are downregulated. A recent paper on the adaptation of M. persicae to tobacco plants sheds light on the mechanism involved (H. Guo et al. 2020). The authors identified a conserved cathepsin B protein that was constitutively upregulated in the salivary glands of a non-adapted to tobacco population compared to an adapted one. Moreover, the protein was shown to be secreted into the plant during feeding, and to interact with a plant defense protein, eliciting this way an effective host response. Our analysis indicated that seven of the 13 differentially expressed cathepsin B and L coding genes are highly expressed in published salivary gland transcriptomes (https://doi.org/10.3389/fevo.2018.00090; https://doi.org/10.1111/1744-7917.12856). Moreover, all seven were predicted to have a SignalP cleavage site indicating that they are likely to be secretory in nature (GitHub repository: https://github.com/KseniaJuravel/B.tabaci_salivary_glands). It can therefore be hypothesized that during long-term adaptation to marginally suitable host plants such as tobacco (M. persicae) and habanero pepper (B. tabaci), the adapting population can improve its fitness by reducing the level/s of specific cathepsin B proteins that elicit plant defenses.

Conclusions

This study adds important insights to the limited knowledge on the performance and genome-wide transcriptional changes that occur during long-term adaptation of phloem-feeding generalist species to marginally suitable hosts, likely to be outside their “comfort zone” of host plants. The performance assays indicated a non-straightforward process characterized by incomplete adaptation, unexpected reduction of fitness on the original host and the presence of antagonistic interactions with a mild-abiotic stressor. If continued for long time periods, this is expected to negatively affect the general (across host plants) fitness of the species and might initiate a specialization process that will further limit the host range of the species (Gould 1979; Wybouw et al. 2015). Although our transcriptomic data revealed some cryptic variation that could be uncovered by the short-term transfer from cotton to habanero pepper, the majority of transcriptional changes in the long-term habanero pepper-adapted population were categorized as “evolved” with no initial plastic response. Changes in a relatively small set of three MFs dominated the process: enhanced formation of cuticle structural constituents, enhanced activity of oxidation-reduction processes involved in neutralization of toxic plant compounds and reduced production of elicitors (cathepsin B proteins) that activate plant defenses. Further studies are required in order to identify the molecular regulatory mechanism/s that control this distinct transcriptional signature of adaptation and determine the possible involvement of complex gene-regulatory networks, molecular chaperon hubs, regulatory microRNAs and epigenetic inheritance in regulating the process.

Materials and Methods

Host Plants and B. tabaci Colonies

Multiple studies have indicated that cotton plants are a highly suitable host for the MEAM1 species of B. tabaci (Horowitz 1986; Zang et al. 2006; Naranjo et al. 2009). The MEAM1 population analyzed in this study was collected in Israel in 2003 (Ashalim, the Negev desert) and maintained since then on cotton plants (Gossypium hirsutum, Pima, cultivar Goliath). It has been previously reported that the MEAM1 species can be collected from commercially grown habanero pepper plants in the Yucatan Peninsula (Ballina-Gomez et al. 2013). However, in a preliminary experiment conducted at the beginning of 2017, the “Ashalim” population of MEAM1 presented an extremely low initial survival rate (∼5%) on habanero pepper (Capsicum chinense cultivar Pepper Magnum-Hot Habanero Orange, courtesy of Genesis seeds Ltd). As this phenomenon repeated itself with other MEAM1 populations, including ones recently-derived from the field (data not shown), we considered habanero pepper as a marginally suitable and stressful host for Israeli populations of MEAM1. In April 2017, 600 individuals were transferred from the original cotton-maintained populations to habanero pepper plants. Free mating was allowed in order to avoid inbreeding depression and genetic drift. Ten generations later, two habanero pepper stock populations of 600 hundred individuals were established. Since then, stock populations (colonies) were maintained on habanero pepper and cotton. Each generation, 600 individuals were collected from each colony and transferred to new cages harboring the same plant host from which the individuals were collected (those collected from cotton were returned to cotton plants and those collected from habanero pepper were returned to habanero pepper plants). The colonies’ cages were maintained in a greenhouse under standard conditions of 28 ± 2°C and a 14:10-hour light: dark cycle.

Performance Assay

Performance evaluation experiments were conducted at generation 17, using four distinct treatments (fig. 1): insects from the suitable host were tested on suitable host plants (cotton to cotton or C→C), insects from the marginally suitable host were tested on marginally suitable host plants (habanero pepper to habanero pepper or hP→hP), insects from the suitable host were tested on marginally suitable host plants (cotton to habanero pepper or C→hP) and insects from the marginally suitable host were tested on suitable host plants (habanero pepper to cotton or hP→C).

The performance experiments were conducted at two temperatures, 24°C or 30°C. These temperatures were selected based on previously developed thermal reproductive-performance curves which indicated that at 30°C, MEAM1 achieves 96.3% of its optimal reproductive performance while at 24°C only 79.9% (Zidon et al. 2016). Therefore, our working assumption was that the two temperatures mimic close-to-optimal (30°C), and mild-stressful (24°C), temperature conditions for reproduction of MEAM1 (Butler et al. 1983; Nava-Camberos et al. 2001).

Five pairs of newly emerged male and female adults were collected into leaf clip-cages, 1.0 cm in diameter, allowing to limit the adults to the abaxial surface of leaves. For each of the four treatments, 16 leaf clip-cages were placed on four leaves per plant, four plants per treatment. After an egg-laying period of 48 hours, the adults were removed. The colonies’ performance was assessed using two performance indices—1) survival, the proportion of emerging adults from the total number of eggs laid and 2) the developmental rate, the weighted mean of days for egg to adult development of surviving individuals. Differences in survival across treatments and temperatures (“treatment”, “temperature” and “treatment” × “temperature” as fixed effects) were analyzed using a generalized linear model with a logit link and a binomial distribution. False discovery rate (FDR) corrections were applied to a priori paired comparisons. Statistical significance was assumed at P ≤ 0.05. The weighted-mean development time data was inversely transformed as the values were not normally distributed or homoscedastic, and the non-parametric Kruskal-Wallis test was applied, followed by sequential Bonferroni comparisons (Sokal & Rohlf 1995). All experiments were conducted in an experimental room under standard conditions of 24 ± 1°C or 30 ± 1°C and a 12:12-hour light: dark cycle. Statistical analysis was conducted using the JMP statistical software program (JMP Pro, version 16.0.0; SAS, Cary, NC, USA). Plotting was performed using the ggplot2 (Wickham 2016), ggsignif (Ahlmann-Eltze 2019) and ggpubr (Kassambara 2020) R packages and subsequently modified by Inkscape (Inkscape Project 2020).

RNA Isolation and Illumina Sequencing

As we considered the 24°C temperature regime a more abiotic stressful environment for MEAM1 (Butler et al. 1983; Nava-Camberos et al. 2001; Tsueda & Tsuchida 2011; Zidon et al. 2016), newly emerged adults (1–7 days old) for transcriptomic analysis were collected only from the 24°C-treated cages, separated by sex (males and females) and pooled by treatment (50 individuals in each biological replicate). The C→C and hP→C treatments had five female replicas, the hP→hP treatment had four female replicas, and the C→hP treatment had three female replicas. The female samples were homogenized immediately upon collection with the first buffer (NucleoSpin RNA XS, Macherey-Nagel) using 1-mm sterilized zirconia beads and two pulses of 60 seconds at 5,000 rpm of a bead-beater (Minilys, Bertin Technologies). The male samples were saved for analyzing the epigenetic effects (DNA methylation) on gene expression (currently conducted). Total RNA was extracted from the homogenates according to the manufacturer’s instructions (NucleoSpin RNA XS, Macherey-Nagel). Library construction and sequencing were performed by Macrogen Inc. TruSeq RNA Library v2 libraries were constructed and sequenced on a NovaSeq 6,000 platform. On average, it produced approximately 30 million 100-bp pair-ended reads per sample (biological replicate) (supplementary tables S4 and S5, Supplementary material online).

Gene Expression Analysis

The raw reads obtained were quality screened using the FASTQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and cleaned with Trimmomatic v0.33 (Bolger et al. 2014) using the following parameters: LEADING:3 TRAILING:3 SLIDINGWINDOW:4:25 MINLEN:98. 6. Reads belonging to B. tabaci bacterial endosymbionts were filtered out by Kraken v2.0.6 (pair-ended mode with default settings), using a custom database (kraken2-build function) harboring Bacteria, Archaea, and Plasmid Refseq databases (downloaded September 2018) (Wood & Salzberg 2014). The filtered reads were used in all subsequent analyses.

The B. tabaci MEAM1 reference genome files (primary assembly and gff3 annotation file) were obtained from the NCBI RefSeq database (GCF_001854935.1 and ASM185493v1, downloaded January 2019). The reference transcripts database was obtained with the RSEM v1.2.28 (Li & Dewey 2011) built-in function (rsem-prepare-reference –hisat2-hca –gff3) using the aforementioned RefSeq files. Finally, reads mapping against the reference and expected counts were obtained using RSEM v1.2.28 (Li & Dewey 2011) and the rsem-calculate-expression built-in function (–hisat2-hca –estimate-rspd –paired-end). The percentage of mapped reads ranged from 74 to 77% (supplementary table S5, Supplementary material online).

Genes that did not have at least 10 reads in all samples were filtered out. The filtered count matrix was used as input for the DESeq2 R package v1.18.1 (Love et al. 2014). Differential expression analysis (one factor, four levels, “local” dispersion estimator) was performed by setting treatment contrasts using the design: ∼ 1 + treatment, with P-values adjusted to the false discovery ratio (FDR) of 0.05 and Log2Fold of ≥1 or ≤ −1 as cut-offs. PCA was performed on normalized count data that went through a variance stabilizing transformation (VST) step (Anders & Huber 2010) using the prcomp function from the stats R package (R Core Team 2020). Statistical analyses of PCA data (the equality of variances in PC1 between experimental treatments) were conducted using an F-test (the var.test function) in the dplyr R package (Wickham et al. 2021). Plotting was performed with the ggplot2 and ggpubr R packages (Wickham et al. 2021).

The two most informative salivary glands transcriptomes, ERS2502869 and SRS5714533, were obtained from the NCBI SRA database for testing if the differentially expressed cathepsin B and L coding genes (a total of 13) are expressed in the salivary glands and harbor secretory signals. Bowtie2 v2.4.2 (with default parameters) (Langmead and Salzberg 2012) was used for mapping the reads against the 13 cathepsin B and L coding genes. Only genes displaying high number of mapped reads and similar coverage across the coding region were considered as salivary gland-expressed genes. Analysis of secretion signal peptides was carried out using signalP (https://services.healthtech.dtu.dk/service.php?SignalP-5.0).

Gene Ontology Enrichment Analysis

Gene Ontology (GO) terms were assigned to B. tabaci reference genes using InterPro v5.33-72.0 (-iprlookup -goterms -pa -dp -appl Pfam, TIGRFAM, Hamap). Enrichment analysis was carried out using the topGO R package v2.34 (Rahnenfuhrer & Alexa 2020). The “weight01” algorithm and fisher’s exact test were applied with cut-off thresholds adjusted to the false discovery ratio (FDR) of 0.05 for P-values and Log2Fold of ≥1 for up-regulated genes or Log2Fold of ≤ −1 for down-regulated genes. Significantly enriched MF GO terms were identified separately for the two studied adaptation processes: 1) differences in adaptation to habanero pepper between the C→C, C→hP and hP→hP treatments; 2) differences in adaptation back to cotton between the hP→hP, hP→C and C→C treatments. Plotting was performed using the ggplot2 (Wickham 2016), ggpubr (Kassambara 2020) and gplots (Warnes et al. 2016) R packages.

Click here for additional data file.