Gene expression plasticity across hosts of an invasive scale insect species

Abstract

For plant-eating insects, we still have only a nascent understanding of the genetic basis of host-use promiscuity. Here, to improve that situation, we investigated host-induced gene expression plasticity in the invasive lobate lac scale insect, Paratachardina pseudolobata (Hemiptera: Keriidae). We were particularly interested in the differential expression of detoxification and effector genes, which are thought to be critical for overcoming a plant’s chemical defenses. We collected RNA samples from P. pseudolobata on three different host plant species, assembled transcriptomes de novo, and identified transcripts with significant host-induced gene expression changes. Gene expression plasticity was pervasive, but the expression of most detoxification and effector genes was insensitive to the host environment. Nevertheless, some types of detoxification genes were more differentially expressed than expected by chance. Moreover, we found evidence of a trade-off between expression of genes involved in primary and secondary metabolism; hosts that induced lower expression of genes for detoxification induced higher expression of genes for growth. Our findings are largely consonant with those of several recently published studies of other plant-eating insect species. Thus, across plant-eating insect species, there may be a common set of gene expression changes that enable host-use promiscuity.

About the Authors:

Nicholas Christodoulides

Contributed equally to this work with: Nicholas Christodoulides, Alex R. Van Dam, Nate B. Hardy

* E-mail: [email protected]

Affiliation: Department of Entomology and Plant Pathology, Auburn University, Auburn, Alabama, United States of America

ORCID http://orcid.org/0000-0002-2424-5753

Alex R. Van Dam

Contributed equally to this work with: Nicholas Christodoulides, Alex R. Van Dam, Nate B. Hardy

Current address: Department of Biology, University of Puerto Rico, Mayagüez, Puerto Rico

Affiliation: Biosynthetic Pathways Engineering, Department of Bioengineering, Denmark Technical University, Søltofts plads, Lyngby, Denmark

Daniel A. Peterson

¶‡ These authors also contributed equally to this work.

Affiliation: Department of Biology and Graduate Program in Organismic and Evolutionary Biology, University of Massachusetts, Amherst, Massachusetts, United States of America

Rasmus John Normand Frandsen

Affiliation: Biosynthetic Pathways Engineering, Department of Bioengineering, Denmark Technical University, Søltofts plads, Lyngby, Denmark

Uffe Hasbro Mortensen

Affiliation: Biosynthetic Pathways Engineering, Department of Bioengineering, Denmark Technical University, Søltofts plads, Lyngby, Denmark

Bent Petersen

Affiliation: Center for Biological Sequence Analysis, Department of Systems Biology, Technical University of Denmark, Kemitorvet, Lyngby, Denmark

Simon Rasmussen

Affiliation: Center for Biological Sequence Analysis, Department of Systems Biology, Technical University of Denmark, Kemitorvet, Lyngby, Denmark

Benjamin B. Normark

¶‡ These authors also contributed equally to this work.

Affiliation: Department of Biology and Graduate Program in Organismic and Evolutionary Biology, University of Massachusetts, Amherst, Massachusetts, United States of America

Nate B. Hardy

Contributed equally to this work with: Nicholas Christodoulides, Alex R. Van Dam, Nate B. Hardy

Affiliation: Department of Entomology and Plant Pathology, Auburn University, Auburn, Alabama, United States of America

Introduction

The overwhelming majority of plant-eating insect species are host-specialists that consume only one or a few closely-related plant species [1, 2]. But plant-eating insect communities predictably include a few species that are host-generalists, some of which can have extremely broad diets [2, 3]. These generalists are critical elements in terrestrial food webs; their high degree of ecological connectedness is a critical part of what makes ecosystems resilient and many specialist niches possible [4]. Generalists are also enriched in invasive and pest insect faunas, which gives them more economic weight than their species richness might suggest [5]. What genetic systems enable host-use promiscuity in plant eating insects? Some key host-use genes and proteins have been identified [6], and in the last few years there has been a number of studies examining how gene expression varies across hosts [7, 8]. Nevertheless, at this stage, few plant-eating insect species and host interactions have been studied, and we still know little with confidence about the genetic architecture of generalism in plant-eating insects.

We typically assume that the main limiters on plant-eating insect diets are plant defenses, and that a plant-eating insect’s diet breadth is determined by how many of those defenses it can overcome [9]. Researchers have paid particular attention to a plant-eating insect’s ability to detoxify plant defensive chemicals with a core suite of enzymes in a handful of large gene families: carboxylesterases, UDP-glycosyltransferases, glutathione-S-transferases, cytochrome oxygenase P450s, glutathione peroxidases, and ATP-binding cassette transporters (ABC transporters) [9, 10, 11]. However, this list may be too restrictive. For example, zinc-binding dehydrogenases have been implicated in plant-defense compensation in peach aphids [12]. And recently, researchers have come to recognize that the fitness of some plant-eating insects on their host plants depends on effectors, that is, secreted macromolecules that attenuate a plant’s inducible defenses [13, 14]. Disrupting a host plant’s inducible defenses may be especially important for phloem-sap-sucking insects—such as aphids, whiteflies and scale insects—as they have close and persistent contact with their hosts [15]. Less attention has been paid to how factors other than plant defense, such as nutrition, could limit a plant-eating insect’s diet, but non-defensive factors may be important.

The recent spate of comparative transcriptomics studies indicates that broad diets are, in fact, facilitated by differential expression of detoxification genes [8, 16, 17]. However, those studies also indicate that detoxification genes are only a small part of the picture. Most of the significant host-induced changes in gene expression do not appear to be related to detoxification. For example, researchers have found consistent host-dependent expression of many genes involved in primary metabolism [8, 18]. Furthermore, some studies have recovered evidence of a trade-off between expression of genes for detoxification and primary metabolism [7, 18]; for some plant-eating insects, increased investment in the expression of detoxification genes may come at the cost of reduced investment in growth and reproduction. But this is not the case for all insects. In the silverleaf whitefly, Bemisia tabaci, we find just the opposite: expression of genes for detoxification and primary metabolism rise and fall together [19]. Thus, the generality of this trade-off is not clear. Also unclear is the extent to which expression of effector genes varies across host plant species (but see [16] for evidence that it is important in host-use adaptation). Effector genes are involved in many biological processes and have many molecular functions. Consequently, there is no effector category in standard gene ontologies, and researchers would not detect enriched differential expression of effectors across hosts using standard enrichment analyses. Furthermore, with one exception [20], the published RNAseq-based differential expression analyses of insect diet variation have studied species with transient connections with their host plants (at least in comparison to sap-sucking insects), and might not be expected to strongly manipulate their host’s induced defenses.

The lobate lac scale, Paratachardina pseudolobata (Hemiptera: Coccoidea: Kerriidae), is a good example of host-generalism (polyphagy) in plant-eating insects, and why it is worth studying. Thought to be native to Southeast Asia [21], the lobate lac scale is invasive and widespread in nearly all areas of Florida. It has been recorded from more than 307 plant species, including several of ecological and economic importance. Females are parthenogenic and sessile phloem feeders. They form dense aggregations on branches and stems, and excrete sugar-rich waste that feeds harmful sooty molds [22]. Some native plant species seem to be particularly vulnerable to lobate lac scale infection, including wax myrtle (Myrica cerifera), a favorite nesting area for many water-wading birds [23]. Controlling lobate lac scale populations has proven difficult; natural enemies are lacking in Florida, and thick resinous wax protects females from many insecticides [21].

In this study we investigated the genetic architecture of host-generalism in Paratachardina pseudolobata. Specifically, we examined gene expression variation across field-collected samples taken from multiple locations and hosts in southern Florida. We were particularly interested in 1) determining how expression of effector genes depends on host species, 2) confirming the importance of the differential expression of detoxification genes across hosts, and 3) assessing if there are trade-offs between the expression of genes involved in primary metabolism and those involved in plant defense compensation.

Materials and methods

Sample collection

P. pseudolobata samples of whole, adult females were collected from four locations in southern Florida (Table 1). No specific permissions were required for these collections, and the insects are not protected or endangered species. Each sample was made up of about twenty specimens. A total of nine samples were taken from three host plants: two samples from Tetrazygia bicolor (Melastomataceae), three samples from Myrsine cubana (Myrsinaceae), and four samples from Psychotria nervosa (Rubiaceae). Specimens were plucked from their host plants and immediately homogenized in tubes filled with Trizol reagent. The ~20 specimens in each tube were pooled.

[Figure omitted. See PDF.]

Table 1. Samples of lobate lac scales.

https://doi.org/10.1371/journal.pone.0176956.t001

RNA extraction, library, preparation, and sequencing

Each RNA sample was homogenized using the Fast Prep FP120 Homogenizer Cell Disrupter (Thermo Fisher) for three pulses of 20 seconds using small steel beads. Then, 700μL of each sample was transferred to a fresh centrifuge tube with 200μL of chloroform and shaken vigorously for 15 seconds and then incubated on ice for 20 minutes, shaking every 10 minutes. Samples were then centrifuged for 15 minutes at 14,000rpm at 4°C. Approximately 400μL of each sample was added to an equal amount of 100% EtOH at 4°C. This mixture was then passed through the standard protocol of the TRIzol^® Plus RNA Purification Kit (Thermo Fisher). RNA yields and qualities were assessed using a Qubit fluorometer and an Agilent 2100 Bioanalyzer at the Denmark Technical University Multi-Assay Core (DMAC). CDNA Libraries were prepared at the University of California Davis Genome Center DNA Technologies Core. Samples were poly-A tailed, normalized, and pooled using Illumina TruSeq adapters. Samples were then loaded across three lanes of a single flow cell and sequenced in 100 base paired-end reads using an Illumina HiSeq3000.

Quality control and transcriptome assembly

We removed adapters and low-quality sequences (quality score cutoff = 25) from the raw sequence reads using Trimmomatic 0.35 [24] and discarded reads less than 36 bp long. We then merged and normalized the quality-trimmed reads by kmer coverage and length using Khmer [25]. We assembled normalized reads into transcriptomes de novo using both Velvet 1.2.08 and Trinity [26, 27]. We constructed the Trinity assembly using default parameters (kmer = 25). We constructed eight Velvet assemblies, each with scaffolding enabled, and a minimum transcript length of 200bp, but with eight different kmer lengths: 27, 29, 31, 33, 35, 43, 53, 63. We then combined all assemblies from both programs into a single Fasta file and passed it through the EvidentialGene pipeline [28] to cluster sequences based on similarity, and find the best consensus transcripts for a final assembly. We assessed the quality of the assembly with BUSCO and Transrate [29, 30]. The Transrate assembly score was 0.30, which is considered a passing score based on the number of input reads mapping to the assembly. BUSCO found that 89% of highly conserved arthropod sequences were present as single-copy or duplicated transcripts in the assembly, suggesting that it is mostly complete. For annotation and gene ontology assignments we passed the assembly through Annocript, BLASTing against the Uniref90, SwissProt and Conserved Domain databases with an e-value cut-off value of 1e-5 [31, 32]. We used these annotations for all subsequent analysis.

Differential expression and GO enrichment analysis

We quantified the number of each transcript represented in the quality-trimmed read libraries by mapping them to the final assembly using eXpress and Bowtie2 [33, 34]. To determine how gene expression in lobate lac scale populations varied across host plant species, we first created a matrix of transcript counts using Trinity scripts. Next, we used the Bioconductor package Limma to convert the transcript counts to log-counts per million and the voom function to model the mean-variance relationship with precision weights across the three host plant species with three pair-wise comparisons [35]. We excluded low-abundance transcripts from the analysis. We used a false discovery rate (FDR) cutoff value of 0.05 to classify transcripts as differentially expressed, and made an expression heat map via Trinity scripts. We also used TopGO with gene ontology (GO) annotations provided by Annocript to perform GO term enrichment analysis [36]. TopGO categorizes differentially expressed transcripts under GO terms relating to general biological processes and molecular functions. We summarized these results with REVIGO [37].

We then looked specifically at the expression of detoxification and effector genes across hosts. For detoxification genes, we examined transcripts which had been annotated as carboxylesterases, UDP-glycosyltransferases, glutathione-S-transferases, cytochrome oxygenase P450s, glutathione peroxidases, and ABC transporters. Investigation of insect effector genes was less straightforward, as they have only recently been characterized, fall into a variety of gene ontology classes, and are unlikely to be in curated protein databases. We began with a list of 67 candidate effector genes in the aphids Acrythosiphon pisum and Myzus persicae provided by [38]. We then extended this list to include 20 predicted effectors from aphid salivary gland transcriptomes [10]. We collated sequences for the genes in this list and identified lobate lac scale homologs in our assembly using the tblastx function in BLAST+ with an e-value cutoff of 0.00001.

Results

Assembly and annotation

Sequencing generated ~120 million paired-end reads per library. The final, consensus assembly contained 113,670 transcripts, with a contig N50 of 1055 base pairs and a longest contig length of 33,906 base pairs. Annotation yielded 82,024 transcripts with at least one database hit. In comparison to published transcriptomes, the lobate lac transcriptome is most similar to those of the two-spotted spider mite, Tetranychus urticae, and the pea aphid, Acrythosiphon pisum. The most diverse genes in the assembly included sugar transporters, protein kinases, Ras family proteins, ubiquitins, and cytochrome P450s. Transcripts annotated as detoxification genes in the assembly included 354 cytochrome P450 oxygenases, 12 carboxylesterases, 96 glutathione S-transferases, 78 UDP-glycosyltransferases, 10 glutathione peroxidases, and 646 ABC transporter transcripts. Additionally, we found 1015 transcripts closely matching candidate aphid effectors.

Differential expression

We identified 2,028 transcripts as differentially expressed across all comparisons. Most differentially expressed transcripts were upregulated on M. cubana relative to the other two plant species, but 280 and 58 genes were upregulated in T. bicolor and P. nervosa respectively (Fig 1). Only 23 transcripts matching detoxification genes were significantly differentially expressed: 3 UDP glycosyltransferases, 4 cytochrome P450s, 2 glutathione S-transferases, 1 glutathione peroxidase, and 13 subunits of ABC transporters. Each was upregulated in M. cubana, except for 2 ABC transporter subunit transcripts and 1 glutathione S-transferase, which were upregulated in T. bicolor compared to P. nervosa. Of the 1015 transcripts which closely resemble possible aphid effectors, 8 were significantly differentially expressed (Table 2). Seven effector homologs were upregulated on M. cubana and one on T. bicolor.

[Figure omitted. See PDF.]

Fig 1. Heatmap showing differential expresssion across host plants.

(Purple: Tetrazygia bicolor, Green: Psychotria nervosa, Red/: Myrsine cubana). Heatmap colors correspond to log2-transformed RPKM (Fragments Per Kilobase of transcript per Million mapped reads) values for each transcript. Yellow genes are upregulated in at least one pairwise comparison between the three host plants and purple genes are downregulated.

https://doi.org/10.1371/journal.pone.0176956.g001

[Figure omitted. See PDF.]

Table 2. Differentially expressed aphid effector homologs and what is known about their function.

https://doi.org/10.1371/journal.pone.0176956.t002

Functional enrichment analysis

The genes upregulated in lobate lac scale feeding on M. cubana were enriched for GO terms related to functions such as ion transport, ATP hydrolysis-coupled proton transport, regulation of gene expression, GTP binding, oxidoreductase activity on NADPH, phosphatase activity, ion binding, and metalloendopeptidase activity (Table 3). In contrast, on P. nervosa, upregulated expression was significantly enriched for macromolecule localization and protein catabolism. Enriched GO terms on T. bicolor included actomyosin structure, cellular component organization, and regulation of gene expression.

[Figure omitted. See PDF.]

Table 3. Functional enrichment analysis.

https://doi.org/10.1371/journal.pone.0176956.t003

Many of the GO terms enriched for differential expression on M. cubana are thought to be involved in detoxification (Table 3). Specific mechanisms are not clear, but increased expression of these genes has been associated with exposure to plant chemical defenses [44]. By contrast, these GO terms are not enriched on T. bicolor and P. nervosa. Instead, on those host plant species, we saw increased expression of genes involved in primary metabolism.

Discussion

The recent comparative genomics work is beginning to show us general features of the genetic architecture of host-use variation in plant-eating insects. In this study, we find additional support for the importance of host-dependent expression of detoxification genes. We also find further evidence that detoxification genes account for only a small part of all host-induced gene expression changes, and that there is a trade-off between expression of genes for detoxification and those for primary metabolism. We recover evidence of a striking diversity of effector proteins in lobate lac scale; more than a thousand transcripts, ~ 1% of all transcripts in the sample, are near matches to one of 87 putative aphid effectors, and several of these transcripts are differentially expressed in the lobate lac scale across host plant species. Nevertheless, effectors as a group were not enriched for differential expression; most of them are expressed at the same level across hosts.

Detoxification genes

Of the 1196 putative detoxification genes that were expressed, only 23 were significantly upregulated on any one host. Pervasive host-insensitivity in the expression of detoxification genes appears to be the rule across plant-eating insects [18, 45]. However, we did find that ABC transporters were more differentially expressed across hosts than expected by chance. The pool of differentially expressed genes was significantly enriched for genes involved in ATP hydrolysis-coupled proton transport (p-value = 1.6e-5) (Table 3), a process carried out by ATPase proteins, which are subunits of ABC transporter proteins. ABC transporter expression is associated with xenobiotic elimination and insecticide resistance in several insect species [44].

Only a few of the genes that we identified a priori as being involved in detoxification were found to be differentially expressed, but perhaps our a priori assignments were too exclusive. In fact, some of the enriched GO terms in lobate lacs scale from M. cubana could be indicative of host-induced differences in the expression of genes involved in detoxification pathways (Table 3). Most notably, oxidoreductase activity on NADPH is often linked to detoxification, as several detoxification proteins use NADPH as an electron acceptor [46]. These proteins include cytochrome P450 oxygenases, and glutathione peroxidases [47]. In other words, we found significantly enriched differential expression of some genes that are adjacent to detoxification genes in metabolic networks. Furthermore, oxidoreductase activity and zinc binding, other enriched GO terms, are functions of zinc-binding dehydrogenases, which appear to be important for plant defense compensation in peach aphids [12]. Thus, if we consider a more inclusive set of genes involved in plant-defense compensation, our expression data provides more support for the importance of host-induced differential expression.

Aphid effectors on M. cubana

Each of the differentially expressed effector transcripts is likely to have a similar function to its closest matching aphid transcript (Table 2). These effectors have several functions. M1 zinc metalloprotease is one of the most well-understood effector proteins. Sap-sucking insects use it to destroy plant signaling defense peptides and improve phloem sap quality by increasing free amino acid content [10]. A second effector that we know something about is calreticulin. Calcium ions transmit information in plant signaling pathways. Calreticulin binds calcium ions and disrupts these signals. This prevents plants from closing off compromised sieve tube elements, which would cut off the supply of phloem sap [10]. The over-expression of a calreticulin homolog may account for the enrichment of calcium binding processes we found in lobate lac scale on M. cubana. A third effector with characterized function is retinol dehydrogenase. It attracts retinols and fatty acids and can bind lipid jasmonate precursors to prevent jasmonate pathway defenses [43]. This may relate to the enriched lipid and phospholipid binding we observe in lobate lac scale on M. cubana.

Induction of effectors and detoxification genes in M. cubana

The great majority of differentially expressed genes were upregulated on M. cubana. This could indicate that M. cubana is simply better defended than the other plant species. Plants in the Myrsinaceae family are known to produce a variety of benzoquinone compounds, characterized by an aromatic ring and two carbonyl groups [48]. Benzoquinones and their derivatives have antibacterial, insecticidal, and anti-feeding properties [49]. Benzoquinones are toxic to several species of aphid as well as the extremely polyphagous red spider mite, which shows that they could be difficult to process by generalist insect species [50]. Many species of Myrsinaceae also accumulate saponins in their tissues. Saponins are cytotoxic chemicals composed of a hydrophobic aglycone (non-sugar) attached to glycoside (sugar residue) [48]. They play an important role in plant defense against insects and parasites, and they also elicit immune responses in people [51]. Myrsinaceous plants can also be induced to produce proteinase inhibitors known as phytostatins in response to wounding [52]. These proteinase inhibitors make their way to the insect gut and inhibit digestion. In comparison to M. cubana, P. nervosa and T. bicolor are thought to be less protected by defensive chemistry [53, 54]

DE of gene expression regulation

We found that some of the most conspicuous changes in gene expression across hosts occurs at genes involved in regulating gene expression. In lobate lac scale on T. bicolor, we found enrichment of genes for gene silencing by RNA. This was previously shown on some hosts of Myzus persicae (in which case miRNA activity is negatively correlated with the abundance of a P450 enzyme) [55]. In lobate lac scales on P. nervosa, we found enriched expression of genes for RNA binding and translation factor activity, both of which are thought to play a role in eukaryotic post-transcriptional regulation [56]. Previous RNA-seq studies of how gene expression differs across hosts of plant-eating insect species have also found differential expression of genes for ribosomal proteins and nucleic acid binding. This kind of differential expression may integrate across GO categories and be indicative of broad-scale expression changes that may nonetheless not amount to significant enrichment of differential expression within GO categories.

Caveats

In this study, some of the gene expression variation that we observed across hosts could be due to biotic and abiotic factors that covary with hosts across sites. For example, gene expression variation could correspond to variation in local light and soil conditions, or assemblages of other plant-eating insects and natural enemies. In future studies, we could minimize the noise caused by such factors by using common garden experimental designs, or sampling much more intensively across wild populations. The strictly bioinformatics approach that we have taken to functional analysis also has limitations. The function of candidate genes for host-generalism could be more accurately characterized with genome editing approaches such as RNA interference. And that may yield insights that could yield new tools to control populations of pests.

Conclusion

For decades, we assumed that host use in a plant-eating insect is limited by their capacity to cope with the chemical defenses of plant species, but we knew little about the details. In the last few years, comparative transcriptomics analyses of plant-eating insects have identified hundreds of genes that are differentially expressed across hosts, including many genes suspected to be involved in the detoxification of plant chemicals. They have also shown us that host-use adaptation is much more complex than activating or silencing the expression of a few detoxification genes, and may entail gene expression trade-offs between plant-defense compensation and primary metabolism. Here, by analyzing gene expression plasticity across three host species of the lobate lac scale insect, we find further support for each of these insights. Moreover, we find evidence that the genes that are directly involved plant-defense-compensation may vary across plant-eating insect species, and may be more diverse than commonly thought. We also find evidence that in some cases, host-defense compensation may depend on differential expression of genes that are adjacent to detoxification enzymes in metabolic pathways. These insights have the potential to pave the way for a more useful theory of host-use adaptation, and for new, more economical and sustainable tools for managing pest insect populations.

Acknowledgments

Thanks to Robin Giblin-Davis, William Kern, Ian Stocks, and Catherine Mannion for helping us find scale insects and learn the southern Florida flora. Thanks also to Mayrolin García-Morales for helping with the field work. NC thanks the Auburn University Department of Entomology and Plant Pathology for research assistantship support. We thank Marlene Dalgaard (DMAC) for use of their Qubit and Bioanalyzer. We thank Andrea Schreier for sample processing at UC Davis. We also thank Lutz Froenicke at UC Davis for advice on the sequencing scheme for this project. Funding: This work was supported by the National Science Foundation [DBI 1306489] awarded to AVD; and the Alabama Agricultural Experiment Station.

Author Contributions

1. Conceptualization: NC AVD NBH DAP BBN RJNF UHM BP SR.

2. Data curation: NC AVD NBH DAP BBN.

3. Formal analysis: NC NBH AVD.

4. Funding acquisition: AVD.

5. Investigation: NC AVD NBH DAP BBN.

6. Methodology: NC AVD NBH DAP BBN.

7. Project administration: NC AVD NBH DAP BBN.

8. Resources: NC AVD NBH DAP BBN RJNF UHM BP SR.

9. Software: NC AVD NBH DAP BBN.

10. Supervision: NC NBH AVD.

11. Validation: NC NBH AVD.

12. Visualization: NC NBH.

13. Writing – original draft: NC.

14. Writing – review & editing: NC AVD NBH DAP BBN.

Citation: Christodoulides N, Van Dam AR, Peterson DA, Frandsen RJN, Mortensen UH, Petersen B, et al. (2017) Gene expression plasticity across hosts of an invasive scale insect species. PLoS ONE 12(5): e0176956. https://doi.org/10.1371/journal.pone.0176956

References

1. Bernays E, Graham M. On the Evolution of Host Specificity in Phytophagous Arthropods. Ecology. 1988;69(4):886–92.

2. Forister ML, Novotny V, Panorska AK, Baje L, Basset Y, Butterill PT, et al. The global distribution of diet breadth in insect herbivores. Vol. 112, Proceedings of the National Academy of Sciences of the United States of America. 2014. 1–6 p.

3. Hardy NB, Peterson DA, Normark BB. Scale insect host ranges are broader in the tropics. Biol Lett. 2015;11:20150924. pmid:26701757

4. Bascompte J. Plant-animal mutualistic networks: the architecture of biodiversity. Appl Math. 2007;101(2):187–203.

5. Ross L, Hardy NB, Okusu A, Normark BB. Large Population Size Predicts the Distribution of Asexuality in Scale Insects. Evolution (N Y). 2012;

6. Czosnek H, Ghanim M. Management of insect pests to agriculture: Lessons learned from deciphering their genome, transcriptome and proteome. Management of Insect Pests to Agriculture: Lessons Learned from Deciphering their Genome, Transcriptome and Proteome. 2016. 1–290.

7. Govind G, Mittapalli O, Griebel T, Allmann S, Böcker S, Baldwin IT. Unbiased transcriptional comparisons of generalist and specialist herbivores feeding on progressively defenseless Nicotiana attenuata plants. PLoS One. 2010;5(1).

8. Ragland GJ, Almskaar K, Vertacnik KL, Gough HM, Feder JL, Hahn DA, et al. Differences in performance and transcriptome-wide gene expression associated with Rhagoletis (Diptera: Tephritidae) larvae feeding in alternate host fruit environments. Mol Ecol. 2015;24(11):2759–76. pmid:25851077

9. Li X, Schuler MA, Berenbaum MR. Molecular Mechanisms of Metabolic Resistance to Synthetic and Natural Xenobiotics. Annu Rev Entomol. 2007;52(1):231–53.

10. Carolan JC, Caragea D, Reardon KT, Mutti NS, Dittmer N, Pappan K, et al. Predicted effector molecules in the salivary secretome of the pea aphid (Acyrthosiphon pisum): A dual transcriptomic/proteomic approach. J Proteome Res. 2011;10(4):1505–18. pmid:21226539

11. You M, Yue Z, He W, Yang X, Yang G, Xie M, et al. A heterozygous moth genome provides insights into herbivory and detoxification. Nat Genet. 2013;45(2):220–5. pmid:23313953

12. Cooper WR, Dillwith JW, Puterka GJ. Salivary proteins of Russian wheat aphid (Hemiptera: Aphididae). Environ Entomol. 2010;39(1):223–31. pmid:20146860

13. Bellafiore S, Shen Z, Rosso M-N, Abad P, Shih P, Briggs SP. Direct Identification of the Meloidogyne incognitaSecretome Reveals Proteins with Host Cell Reprogramming Potential. PLoS. 2008;

14. Harris MO, Friesen TL, Xu SS, Chen MS, Giron D, Stuart JJ. Pivoting from Arabidopsis to wheat to understand how agricultural plants integrate responses to biotic stress. Vol. 66, Journal of Experimental Botany. 2015. p. 513–31.

15. Will T, Furch ACU, Zimmermann MR. How phloem-feeding insects face the challenge of phloem-located defenses. Front Plant Sci. 2013;4(August):336. pmid:24009620

16. Wybouw N, Zhurov V, Martel C, Bruinsma KA, Hendrickx F, Grbic V, et al. Adaptation of a polyphagous herbivore to a novel host plant extensively shapes the transcriptome of herbivore and host. Mol Ecol. 2015;24(18):4647–63. pmid:26211543

17. De La Paz Celorio-Mancera M, Wheat CW, Vogel H, Söderlind L, Janz N, Nylin S. Mechanisms of macroevolution: Polyphagous plasticity in butterfly larvae revealed by RNA-Seq. Mol Ecol. 2013;22(19):4884–95. pmid:23952264

18. Hoang K, Matzkin LM, Bono JM. Transcriptional variation associated with cactus host plant adaptation in Drosophila mettleri populations. Mol Ecol. 2015;24(20):5186–99. pmid:26384860

19. Alon M, Elbaz M, Ben-Zvi MM, Feldmesser E, Vainstein A, Morin S. Insights into the transcriptomics of polyphagy: Bemisia tabaci adaptability to phenylpropanoids involves coordinated expression of defense and metabolic genes. Insect Biochem Mol Biol. 2012;42(4):251–63. pmid:22212826

20. Eyres I, Jaquiry J, Sugio A, Duvaux L, Gharbi K, Zhou JJ, et al. Differential gene expression according to race and host plant in the pea aphid. Mol Ecol. 2016;25(17):4197–215. pmid:27474484

21. Pemberton RW. Potential for biological control of the lobate lac scale, Paratachardina lobata lobata (Hemiptera: Kerriidae). Florida Entomol. 2003;86:353–60.

22. Howard FW, Pemberton R, Hamon A, Hodges GS, Steinberg B, Mannion M, et al. (Chamberlin) (Hemiptera: Sternorrhyncha: Coccoidea: Kerriidae) 1. Univ Florida. 2014;1–14.

23. Pemberton RW. Invasion of Paratachardina lobata lobata (Hemiptera: Kerriidae) in South Florida: A snapshot sample of an infestation in a residental yard. In: Florida Entomologist. 2003. p. 373–7.

24. Bolger A, Lohse M U B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114. pmid:24695404

25. Brown CT, Howe A, Zhang Q, Pyrkosz AB, Brom TH. A Reference-Free Algorithm for Computational Normalization of Shotgun Sequencing Data. arXiv. 2012;1203.4802(v2):1–18.

26. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512. pmid:23845962

27. Zerbino DR, Birney E. Gibbs RA, and Buetow KH 2005. Velvet: algorithms for de novo short read assembly using de Brujin graphs. Genome Res. 2008 May;18(5):821–9. pmid:18349386

28. Gilbert D. EvidentialGene: mRNA Transcript Assembly Software [Internet]. Genome Informatics Lab of Indiana University Biology Department. 2013. http://arthropods.eugenes.org/genes2/about/EvidentialGene_trassembly_pipe.html

29. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva E V., Zdobnov EM. BUSCO: Assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2. pmid:26059717

30. Smith-Unna RD, Boursnell C, Patro R, Hibberd JM, Kelly S. TransRate: reference free quality assessment of de-novo transcriptome assemblies. bioRxiv. 2015;21626.

31. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST plus: architecture and applications. BMC Bioinformatics. BioMed Central; 2009 Dec 15;10(421):1.

32. Musacchia F, Basu S, Petrosino G, Salvemini M, Sanges R. Annocript: A flexible pipeline for the annotation of transcriptomes able to identify putative long noncoding RNAs. Bioinformatics. 2015;31(13):2199–201. pmid:25701574

33. Roberts A, Trapnell C, Donaghey J, Rinn JL, Pachter L. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. 2011;12(3):R22. pmid:21410973

34. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9. pmid:22388286

35. Law CW, Chen Y, Shi W, Smyth GK. voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15:R29. pmid:24485249

36. Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for gene ontology. Vol. 2, R package version. 2016.

37. Supek F, Bosnjak M, Skunca N, Smuc T, S . REVIGO summarizes and visualizes long lists of Gene Ontology terms. PLo e2.1800;6.

38. Bos JIB, Prince D, Pitino M, Maffei ME, Win J, Hogenhout SA. A functional genomics approach identifies candidate effectors from the aphid species Myzus persicae (green peach aphid). PLoS Genet. 2010;6(11).

39. Geldhof P, Vercauteren I, Knox D, Demaere V, Van Zeveren A, Berx G, et al. Protein disulphide isomerase of Ostertagia ostertagi: An excretory-secretory product of L4 and adult worms? Int J Parasitol. 2003;33(2):129–36. pmid:12633650

40. Wang W, Dai H, Zhang Y, Chandrasekar R, Luo L, Hiromasa Y, et al. Armet is an effector protein mediating aphid-plant interactions. FASEB J. 2015;29(5):2032–45. pmid:25678626

41. Carolan JC, Fitzroy CIJ, Ashton PD, Douglas AE, Wilkinsonl TL. The secreted salivary proteome of the pea aphid Acyrthosiphon pisum characterised by mass spectrometry. Proteomics. 2009;9(9):2457–67. pmid:19402045

42. Nicholson SJ, Hartson SD, Puterka GJ. Proteomic analysis of secreted saliva from Russian Wheat Aphid (Diuraphis noxia Kurd.) biotypes that differ in virulence to wheat. J Proteomics. 2012;75(7):2252–68. pmid:22348819

43. Iberkleid I, Vieira P, de Almeida Engler J, Firester K, Spiegel Y, Horowitz SB. Fatty Acid-and Retinol-Binding Protein, Mj-FAR-1 Induces Tomato Host Susceptibility to Root-Knot Nematodes. PLoS One. 2013;8(5).

44. Vogel H, Musser RO, de la Paz Celorio-Mancera M. Transcriptome Responses in Herbivorous Insects Towards Host Plant and Toxin Feeding. Annu Plant Rev Insect-Plant Interact. 2014;47:197–233.

45. Strode C, Wondji CS, David JP, Hawkes NJ, Lumjuan N, Nelson DR, et al. Genomic analysis of detoxification genes in the mosquito Aedes aegypti. Insect Biochem Mol Biol. 2008;38(1):113–23. pmid:18070670

46. Miles PW. O . The significance of antioxidants in the aphid plant interaction—the redox hypothesis.pdf. Entomol Exp Appl. 1993;275–83.

47. Mao W, Berhow M, Zangerl A, Mcgovern J, Berenbaum M. Cytochrome P450-Mediated Metabolism of Xanthotoxin by Papilio multicaudatus. Vol. 32, Journal of Chemical Ecology. 2006. p. 523–36.

48. Lukhele T, Krause, M C-SRW, Olivier MDK. Isolation, Characterisation and Biological Activity of Some Compounds From Rapanea Melanophloeos (L.) Mez. Mez. 2009.

49. Midiwo J, Mwangi R. Insect antifeedant, growth-inhibiting and larvicidal compounds from Rapanea melanphloes (Myrsinaceae). Int J. 1995;

50. Akhtar Y, Isman MB, Lee CH, Lee SG, Lee HS. Toxicity of quinones against two-spotted spider mite and three species of aphids in laboratory and greenhouse conditions. Ind Crops Prod. 2012;37(1):536–41.

51. Augustin JM, Kuzina V, Andersen SB, Bak S. Molecular activities, biosynthesis and evolution of triterpenoid saponins. Vol. 72, Phytochemistry. 2011. p. 435–57.

52. Volpicella M, Leoni C, Costanza A, De Leo F, Gallerani R, Ceci LR. Cystatins, serpins and other families of protease inhibitors in plants. Curr Protein Pept Sci. 2011;12(5):386–98. pmid:21418017

53. Adonizio A, Downum K, Bennett B, Mathee K. Anti-quorum sensing activity of medicinal plants in southern Florida. J Ethnopharmacol. 2006;105(3):427–35. pmid:16406418

54. Melendez PA, Capriles VA. Antibacterial properties of tropical plants from Puerto Rico. Phytomedicine. 2006;13(4):272–6. pmid:16492531

55. Peng T, Pan Y, Gao X, Xi J, Zhang L, Ma K, et al. Reduced abundance of the CYP6CY3-targeting let-7 and miR-100 miRNAs accounts for host adaptation of Myzus persicae nicotianae. Insect Biochem Mol Biol. 2016;75:89–97. pmid:27318250

56. Glisovic T, Bachorik JL, Yong J, Dreyfuss G. RNA-binding proteins and post-transcriptional gene regulation. FEBS Lett. 2008;582(14):1977–86. pmid:18342629