Transcriptome Analysis of Green Peach Aphid (Myzus persicae): Insight into Developmental Regulation and Inter-Species Divergence

Green peach aphid (Myzus persicae) and pea aphid (Acyrthosiphon pisum) are two phylogenetically closely related agricultural pests. While pea aphid is restricted to Fabaceae, green peach aphid feeds on hundreds of plant species from more than 40 families. Transcriptome comparison could shed light on the genetic factors underlying the difference in host range between the two species. Furthermore, a large scale study contrasting gene expression between immature nymphs and fully developed adult aphids would fill a previous knowledge gap. Here, we obtained transcriptomic sequences of green peach aphid nymphs and adults, respectively, using Illumina sequencing technology. A total of 2244 genes were found to be differentially expressed between the two developmental stages, many of which were associated with detoxification, hormone production, cuticle formation, metabolism, food digestion, and absorption. When searched against publically available pea aphid mRNA sequences, 13,752 unigenes were found to have no homologous counterparts. Interestingly, many of these unigenes that could be annotated in other databases were involved in the “xenobiotics biodegradation and metabolism” pathway, suggesting the two aphids differ in their adaptation to secondary metabolites of host plants. Conversely, 3989 orthologous gene pairs between the two species were subjected to calculations of synonymous and nonsynonymous substitutions, and 148 of the genes potentially evolved in response to positive selection. Some of these genes were predicted to be associated with insect-plant interactions. Our study has revealed certain molecular events related to aphid development, and provided some insight into biological variations in two aphid species, possibly as a result of host plant adaptation.


INTRODUCTION
Aphids (Insecta: Hemiptera), a group of economically important insect pests that consume plant phloem sap, cause substantial losses of crop yield by direct feeding on host plants and by vectoring plant viruses (Dixon, 1998). More than 450 species within Aphididae attack agricultural and horticultural plants, of which over 100 are categorized as significant and economically important pests (Blackman and Eastop, 1984). While some aphids are specific to plant species in a single taxonomic family, others have an exceptionally broad host range across many plant families. Green peach aphid (Myzus persicae) is a generalist with a host range comprising 40 different plant families including Brassicaceae, Solanaceae, and Fabaceae. Moreover, it is the most versatile viral vector, capable of transmitting more than 100 plant viruses (Ramsey et al., 2007). In contrast, pea aphid (Acyrthosiphon pisum) feeds specifically on legumes. Despite different feeding habits, they are both classified in the tribe Macrosiphini within the subfamily Aphidinae (von Dohlen et al., 2006). The close relationship between the two aphids is further supported by analysis of mitochondrial and nuclear sequences as well as transcriptomic sequence comparisons (Ramsey et al., 2007;Kim and Lee, 2008). Due to the difference in host range, green peach aphids most likely ingest toxic metabolites that pea aphids would not normally encounter, such as glucosinolates in Brassicaceae and alkaloids in Solanaceae, necessitating a more complex metabolic system (Ramsey et al., 2010).
Hemipteran immature nymphs and fully developed adults sometimes differ in their feeding behavior. Lygus hesperus nymphs prefer developing cotton squares, whereas adults prefer vegetative structures (Snodgrass, 1998). In three spittlebug species (Aeneolamia varia, A. reducta and Zulia carbonaria), foliage-feeding adults are more capable of feeding upon resistant hybrid crops than root-and stem-feeding nymphs (Cardona et al., 2010). Besides host and tissue preferences, quantity of food intake can vary (Banks and Macaulay, 1965). Profiling in nymphal and adult transcriptomes could reveal biological properties that are developmental stage-specific. In Asian citrus psyllid (Diaphorina citri) for instance, the transcriptome comparison revealed distinct patterns of protein and energy requirements between nymphs and adults (Vyas et al., 2015). This approach has also identified differentially expressed resistance/detoxification genes, e.g., cytochrome P450, glutathione S-transferase (GST), and ATP-binding cassette transporter genes from two developmental stages of a thiamethoxam-resistant strain of whitefly (Yang et al., 2013). Contrasting gene expression among different insect developmental stages on a large scale can not only shed light on development modulation, reproduction, and developmental stage-specific interaction with host plant, xenobiotics, and invading microbes, but can also facilitate the improvement of pest management strategies (Yang et al., 2013;Tian et al., 2015;Vyas et al., 2015). However, stage-specific gene expression in immature nymphs and fully developed adults has not yet been characterized in aphids.
While comparative genomic sequence analysis has furnished tremendous information regarding genetic factors underlying inter-species divergence (Chinwalla et al., 2002;Kaufman et al., 2002;Kirkness et al., 2003;Zdobnov and Bork, 2007;Arensburger et al., 2010;Bonasio et al., 2010;Werren et al., 2010), an increasing number of studies have applied RNA-seq for this purpose, particularly in species whose genome sequences are unavailable. For example, transcriptomic comparisons have been performed between different aphids, A. pisum vs. Sitobion avenae (Wang et al., 2014), whitefly (Bemisia tabaci) species complexes Middle East-Asia Minor 1 vs. Mediterranean (Wang et al., 2011), ranid frogs Rana chensinensis vs. Rana kukunoris , ornamental primrose species Primula poissonii vs. Primula wilsonii (Zhang L. et al., 2013), and fishes, Erythroculter ilishaeformis vs. Danio rerio (Ren et al., 2014). Comparisons among pea aphid, green peach aphid and grain aphid (S. avenae) have enabled investigation of the transcriptome evolution and understanding of the differences in host plant adaptation and insecticide resistance among them (Ollivier et al., 2010;Ramsey et al., 2010;Wang et al., 2014). Between grain aphid and pea aphid 340 gene orthologs are considered to be under positive selection based on the rates of nonsynonymous (Ka) and synonymous (Ks) substitutions (Wang et al., 2014). Such orthologs were also identified when Ollivier et al. (2010) compared coding sequences (CDSs) derived from the genome sequence of pea aphid and EST database derived from 5 tissues of green peach aphids reared on 5 host plants (Ramsey et al., 2007). Later, Ramsey et al. (2010) sequenced the transcriptome from mixed stages of green peach aphids using 454 pyrosequencing. Besides the reads mapped to the existing ESTs, they obtained 47,832 additional unigenes with a mean length of 160 bp, from which they identified more detoxification genes in green peach aphid than in pea aphid (Ramsey et al., 2010). However, limited transcriptomic information may not fully reflect the divergence between the two species.
In this study, we performed transcriptomic sequencing of green peach aphid nymphs and adults using Illumina RNAseq technology, de novo assembled sequencing reads, and annotated the resulting unigenes. Gene expression profiling between nymphs and adults identified genes potentially involved in development modulation. Furthermore, comparative transcriptomic analyses identified genes unique to green peach aphid (relative to pea aphid) and orthologous gene pairs under positive selection. Data analysis has helped expose certain genetic factors underlying host plant adaptation by the two destructive aphid species.

Plant Growth and Insect Rearing
Arabidopsis ecotype Col-0 plants were grown in LP5 potting medium (Sun Gro Horticulture, Agawam, WA, USA) in an environmental chamber at 23 • C (day)/21 • C (night), 65% relative humidity (RH), and a photosynthetic photon flux density of 88 µmol m −2 s −1 with a 12-h light/12-h dark photoperiod. The green peach aphid (a tobacco-adapted red lineage from Dr. Georg Jander, Boyce Thompson Institute for Plant Research, Cornell University) had been maintained on Col-0 for over 40 generations. Age-synchronized nymphs and adults were subjected to RNA extraction as described below.

RNA Isolation and Transcriptome Sequencing
Neonate nymphs (within 16 h) were placed on 4-week-old Col-0 plants for 4 or 8 days respectively. Sixty 4-day-old nymphs and 60 8-day-old adults were collected, immediately frozen in liquid nitrogen, and stored at −80 • C for RNA extraction. Three independent biological replicates were performed for transcriptome sequencing analysis.
Total RNA was extracted with TRIzol Reagent (Invitrogen, Carlsbad, CA, USA). RNase-Free DNase (Qiagen, Valencia, CA, USA) was added to remove residual DNA. Samples were then further purified using RNeasy Mini Kit (Qiagen) according to the manufacturer's instructions. Purified total RNA samples were quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and qualified by Agilent Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Transcriptome sequencing was performed on an Illumina HiSeq 2500 platform with 125-nucleotide (nt) paired-end reads at Texas A&M AgriLife Genomics and Bioinformatics Services (College Station, TX, USA).

Sequence Assembly and Annotation
After trimming the adaptor sequences and removing short or low-quality reads (>5% unknown nucleotides or more than 20% nts with >10% error rate), the processed reads were assembled using Trinity software (Trinity Software, Inc., Plymouth, NH, USA) and clustered with TGICL Clustering tools (The Institute for Genomic Research, Rockville, MD, USA) (Pertea et al., 2003;Grabherr et al., 2011). The publically available databases, NCBI non-redundant (Nr), NCBI non-redundant nucleotide (Nt), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Cluster of Orthologous Groups of proteins (COG) were used to perform BLAST analyses to annotate the functions of these assembled unigenes (E-value cutoff of 10 −5 ). Blast2GO software (http://www.geneontology.org) was used for gene ontology (GO) annotations (Conesa et al., 2005).

Differential Gene Expression and RT-qPCR Confirmation
Genes differentially expressed between nymphs and adults were identified based on Fragments Per Kilobase per Million mapped reads (FPKM) values, which adjusts the number of fragments mapped to a transcript by the total number of fragments mapped to all unigenes and the length of the transcript (Mortazavi et al., 2008;Ji et al., 2013). The false discovery rate (FDR) was used for the P-values in multiple tests and analyses. A FDR ≤ 0.001 and an absolute value of the log 2 ratio ≥ 1 provided significance threshold for gene expression differences.
To validate the FPKM analysis, expression of 20 selected genes were measured in nymphs and adults by RT-qPCR. For each total RNA sample, 2 µg RNA was used to synthesize cDNAs with random hexamer primers (Invitrogen) and M-MuLV reverse transcriptase (New England Biolabs, Beverly, MA, USA). qPCR reactions were performed using Power SYBR Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's protocol and run on the CFX384 TM Real Time System (BioRad, Hercules, CA, USA). Dissociation curve analyses were performed to ensure amplification specificity. Mean fold change in gene expression was calculated as described previously (Chi et al., 2011). Primer sequences are provided in Table S1. The 18S rRNA gene of green peach aphid (Acc. No. AF487712.1) was amplified as the internal control.

Functional Analysis of Differentially Expressed Unigenes
GO enrichment analysis was performed to recognize the main biological functions of differentially expressed unigenes. The hypergeometric test was performed to find significantly enriched GO terms in differentially expressed unigenes compared to the whole reference transcriptome background (Su et al., 2012;Ji et al., 2013). The P-value was calculated with the formula: where N and n are defined as the number of genes in the transcriptome and differentially expressed genes with GO annotations, respectively. The variables M and m represent the gene number in the transcriptome annotated to a certain GO term and differentially expressed genes within the group (M-m ≥ 0), respectively. The calculated P-value was subjected to Bonferroni correction. GO terms with corrected P-value, i.e., Q < 0.05 were considered significantly enriched. KEGG analyses were performed to identify significantly enriched pathways represented by differentially expressed unigenes. The hypergeometric test was used in a similar way to that for GO enrichment analysis and the terms with Q < 0.05 were determined as enriched pathways.

Ka and Ks Analyses
To predict CDS regions, unigenes were first aligned by BLAST analyses with E-value cutoff of 10 −5 to public databases in the priority order of Nr, Swiss-Prot, KEGG, and COG. Coding regions with the best match in BLAST were considered to be the CDS. Unigenes unable to be aligned to any databases were scanned by ESTScan, which may predict some coding regions. The CDSs of pea aphid were predicted from the mRNA sequence data (https://www.aphidbase.com/aphidbase/content/download/ 3250/33670/file/aphidbase_2.1b_mRNA.fasta.bz2).
After filtering the redundant CDSs that may result from alternative splicing, predicted CDSs of the two aphid species were used to identify orthologous genes using OrthoMCL (Li et al., 2003). Only single-copy ortholog pairs longer than 150 bp were considered as putative orthologous gene pairs. Ka, Ks, and Ka/Ks-values were computed using the YN method implemented in the software KaKs Calculator Version 1.2 (Yang and Nielsen, 2000;Wang et al., 2011Wang et al., , 2014. As the sequencing errors were distributed among synonymous and non-synonymous sites at equal frequencies, they were not expected to strongly influence the results of analyses (Tiffin and Hahn, 2002;Wang et al., 2011Wang et al., , 2014.

Illumina Sequencing Analysis and De novo Assembly
High-throughput RNA-seq generated the most extensive current transcriptome for the green peach aphid. After quality checks, about 74.1, 74.0, and 74.5 million reads were obtained from the three replicates of nymphs and 74.6, 76.0, and 74.3 million reads from adults ( Table 1). All reads were deposited in the NCBI Short Read Archive (SRA, the accession number SRP073458). The reads were assembled into 89,944, 85,416, and 82,810 contigs with mean lengths of 474, 502, and 460 nt for nymphs and 81,641, 78,710, and 87,354 contigs with mean lengths of 472, 484, and 464 nt for adults ( Table 1). Using paired-end joining and gapfilling, these contigs were finally assembled into a total of 62,627 consensus sequences with a mean length of 1460 nt. GC contents were 39.00% for nymphs and 39.63% for adults, comparable to that of the pea aphid (38.80%) (Wang et al., 2014).
GO assignments were used to classify the functions of the predicted unigenes; 14,260 sequences were categorized into 46 GO terms consisting of three domains: biological process, cellular component and molecular function (Figure 1). The most abundantly expressed genes in "biological process" were involved in cellular process (9028), single-organism process (7075), and metabolic process (6557). In "molecular function, " genes involved in catalytic (6894), binding (6678), and transporter (1137) activities were most abundantly expressed (Figure 1).
To better understand the biological pathways that are active in the green peach aphid, we mapped all sequences to the canonical reference pathways in the KEGG database. As a result, 23,695 sequences were assigned to 187 insect-related KEGG pathways (Table S3), with 3286 unigenes (15.47%) being involved in metabolic pathways. These annotations could be useful for further investigation of specific processes, functions and pathways.

Comparison of Gene Expression Profiles between Nymphs and Adults
When different developmental stages were compared, 1639 genes showed higher expression in nymphs and 605 higher in adults (Figure 2, Table S4). We performed RT-qPCR on selected genes to validate these gene expression data. Of the 20 selected genes, 18 were in agreement with RNA-seq results, suggesting good quality of transcriptomic analysis (Table S5). To gain insight into the major biological pathways represented by the differentially expressed genes, 21 enriched insect-related pathways (Q < 0.05) were identified using the hypergeometric test ( Table 2); 14 were associated with "metabolism" and 3 with "digestive system, " suggesting differential metabolic and digestive activities between nymphs and adults (Banks and Macaulay, 1965;Randolph et al., 1975). The most enriched pathway being "metabolism of xenobiotics by cytochrome P450" is intriguing because it may reflect developmental stage-specific interaction with the host plant. Presumably, nymphal, and adult aphids ingest different amounts of allelochemicals, given that more detoxification genes, e.g., 16 of the 23 differential P450 genes, and all differential esterase (6) and GST (1) genes, were expressed in higher abundance in adults ( Table 3). Developmental stage-dependent variations in expression patterns have often been observed in the detoxification genes (Harrison et al., 2001;Strode et al., 2006;Yang et al., 2013). High expression of CYP321B1 is detected in the late larval stage of tobacco cutworm (Spodoptera litura) (Wang et al., 2016). In B. tabaci, relatively high expression of CYP6CM1 is found in adults, correlating with the observation that specific resistance to neonicotinoid imidacloprid is largely restricted to adults (Nauen et al., 2008;Jones et al., 2011). Similarly, high expression of CYP6P9 in adults of Anopheles funestus, but not in larvae, explains the adult resistance (Amenya et al., 2008). In a pyrethroid resistant strain of Anopheles gambiae, CYP6Z1 is expressed in adults but undetectable in larvae or pupae (Nikou et al., 2003). Direct correlation between expression levels of detoxification genes at different developmental stages and resistance to pesticides is also exemplified by the beet webworm (Pyrausta sticticalis) (Leonova and Slynko, 2004) and citrus red mite (Panonychus citri) (Liao et al., 2013;Zhang K. et al., 2013). Banks and Macaulay (1965) reported that adult aphids have higher food consumption than nymphs. Ahmad (1982) stated that increased amounts of dietary allelochemicals due to increased food consumption may explain elevated P450mediated metabolic activity. In parallel, green peach aphid adults likely ingest more plant materials, thus more allelochemicals from host plants, necessitating higher detoxification capacity.
The differentially expressed genes were also assigned to 20 GO enriched functional groups; ontology distributions are shown in Figure 3. Enriched in the "biological process" and "molecular function" include cuticle formation-related groups such as "structural constituent of cuticle, " "chitin-based cuticle attachment to epithelium" and "molting cycle, chitin-based cuticle." The insect cuticle, composed of chitin and cuticle proteins, not only supports and maintains the physical structure, but also serves as a natural barrier against adverse external impacts (Andersen et al., 1995). Cuticle protein comparisons among insects at different developmental stages show that, rather than being an inert structure, the insect cuticle is developmentally modified (Chihara et al., 1982;Dombrovsky et al., 2003). Consistent with these findings, among the 81 differentially expressed transcripts of cuticular proteins and their precursors we detected, 79 were highly expressed in  nymphs ( Table 3). Insects of this developmental stage repeatedly shed their cuticles and replace them with new layers, thus their cuticle biosynthesis is likely more active. No doubt, hormones play an essential role in insect ecdysis. Enrichment of the "steroid hormone biosynthesis" pathway among the differential genes ( Table 2) supports this notion. The major steroid hormone ecdysone plays an essential role in larval ecdysis, a process mediated by hormones, such as ecdysone and ecdysis triggering hormone (ETH) (Robbins et al., 1968;Ewer et al., 1997). Interestingly, the ETH-encoding gene Unigene5077 was highly expressed in green peach aphid nymphs ( Table 3).

Transcriptomic Divergences between Green Peach Aphid and Pea Aphid
Transcriptome comparisons of different aphid species could provide useful information in understanding transcriptome evolution and the genetic factors underlying the biological divergence of these species. To identify genes specific to green peach aphid (relative to pea aphid), we compared the transcriptome we obtained in this study with publically available mRNA sequence data of pea aphid. tBLASTx identified homologous pea aphid mRNAs for 41,912 of our unigenes, leaving 20,595 having no hits. After removing sequences shorter than 250 bp (too FIGURE 2 | Fold change distribution of green peach aphid unigenes between nymphs and adults. The x-axis shows the fold change (log 2 ratio) of gene expression in nymphs compared to adults. |Log 2 | values of 2244 unigenes are higher than 1, indicating potential importance during developmental transition. short to be translated into polypeptides meaningful for comparisons) and BLASTn hits from pea aphid mRNAs and Nt databases, the remaining 13,752 were considered green peach aphid-specific unigenes under the rearing conditions described (Table S6).   Arabidopsis was selected as our host plant because it is readily consumed by green peach aphid. Its short life cycle, abundant genetic resources and well developed RNAi technique (Ramsey et al., 2007;Pitino et al., 2011;Bhatia et al., 2012;Elzinga et al., 2014) can greatly facilitate our more in-depth studies of candidate genes derived from the current study. One caveat however, is that choice of hosts may impact aphid gene expression. Few studies have been conducted to compare transcriptome profiles of the same insect species feeding on different host plants, but some information is available on differential gene expression of insect populations reared on different varieties/lines of the same host plant species (Ji et al., 2013;Bansal et al., 2014;Yu et al., 2014). It appears that the vast majority of genes are present (but likely varied in expression level) among different insect populations, and that genes solely expressed in one population are rare. Whether this observation can be extended to insect populations feeding on different host plant species is yet to be determined.
The Nr, Nt, Swiss-Prot, COG, KEGG, and GO annotations of green peach aphid-specific unigenes were then performed (Table S6). Only 4.52% were predicted to have defined functions (Table 4), and functions of the remaining sequences need further study in the future. Likewise, KEGG classification identified only 30 unigenes, the most predominant group being "xenobiotics biodegradation and metabolism" (13.33%) (Figure 4). This finding correlates well with the fact that green peach aphids feed on a wider variety of plant species, and may have to FIGURE 3 | Significantly enriched GO categories among the differentially expressed genes between nymphs and adults. GO categories with Q < 0.05 were considered significantly enriched. Classification consists of three domains: biological process, cellular component and molecular function. The y-axis shows the value of −log 10 Q of the category. The GO term with highest −log 10 Q was determined the most significant enrichment.

Ka and Ks Analysis between Green Peach Aphid and Pea Aphid
Contrasting with the above analysis where the focus was on genes unique to green peach aphid, here we concentrated on singlecopy orthologous genes between the two aphids. From the 33,963 green peach aphid CDSs (mean length, 1275 bp) derived from our RNA-seq, 3989 that had one-to-one orthologs in pea aphid CDSs were identified, and 3824 contained both substitution types, from which Ka/Ks ratios were calculated (Table S7). The Ka/Ks ratio provides information about the evolutionary forces operating on a particular coding gene and has been widely used to measure the intensity and mode of selection; Ka/Ks = 1 indicates a neutral evolution; Ka/Ks < 1 suggests that nonsynonymous mutations are deleterious and purged from the population; Ka/Ks>1 indicates that nonsynonymous amino acid substitutions offer fitness advantages and are fixed in the population at a higher rate than synonymous substitutions (Hurst, 2002). However, this cutoff value for positive selection has recently been adjusted to 0.5 by Swanson et al. (2004). They found that 15 of 16 genes with 0.5 < Ka/Ks < 1 showed statistical evidence for adaptive evolution (Swanson et al., 2004). Since then, this new value has been adopted for "positive selection" determination in many studies (Kelleher et al., 2007;Elmer et al., 2010;Yang et al., 2012;Zhang L. et al., 2013;Ren et al., 2014;Cheng et al., 2015;Mu et al., 2015;He et al., 2016;Pereira et al., 2016). In our study, a total of 24 pairs of orthologs had a Ka/Ks ratio greater than 1, and 124 had a Ka/Ks ratio between 0.5 and 1 (Table S8).
Relative to the earlier study by Ollivier et al. (2010), our CDS construction is more complete than that of EST-based (33,963 CDSs, 1275 bp mean length vs. 6652 CDSs, 667 bp mean length), due to improvements in sequencing technology. Nevertheless, some putative orthologs under positive selection were identified by both studies, such as C002 (Table S8). Other genes related to insect-plant interactions include those encoding mucins (Ka/Ks = 1.09 and 0.94), the essential components of peritrophic matrix. Fast-evolving mucin proteins presumably contribute to aphid adaptation to different dietary pro-oxidants, phenolic, and lipophilic xenobiotics associated with their respective host plants (Hiraishi et al., 1991;Felton and Summers, 1995;Barbehenn, 1999Barbehenn, , 2001Barbehenn and Stannard, 2004;Hegedus et al., 2009). Likewise, the homolog of salivary protein gene Me17 (Ka/Ks = 0.69), identified in multiple aphid species with dissimilar plant host ranges, is thought to play important roles as the effector in suppressing defense responses in different host plants and in promoting aphid colonization (Atamian et al., 2013;Pitino and Hogenhout, 2013;Elzinga et al., 2014). Nicotinic acetylcholine receptors (nAChR) in insects are often the target sites for naturally occurring and synthetic insecticides (Millar and Denholm, 2007;Bass et al., 2011). A high mutation rate in the nAChR β-2 subunit (Ka/Ks = 0.55) could help green peach aphids adapt to tobacco and become resistant to nicotine, as is the strain used in this study (Devine et al., 1996;Nauen et al., 1996). Another interesting ortholog pair encode odorant-binding protein 10 (OBP10) (Ka/Ks = 0.52). Nucleotide and amino acid sequence comparisons between the two species indicated that all substitutions occurred in the predicted mature protein region, and 19 of the 62 substitutions resulted in 12 hydrophilic and hydrophobic amino acid conversions ( Figure S1). Sun et al. (2012) observed that the two aphid species showed similar as well as dissimilar behavioral responses to certain tested odors. Presumably, fast evolution in OBPs could contribute to the change in their binding activity, which in turn could facilitate host shift or impact host range (Matsuo et al., 2007;Sun et al., 2012).

CONCLUSIONS
Our RNA-seq data have increased molecular resources available for the green peach aphid, a major agricultural pest as well as a biological model for insect-plant interaction studies. The transcriptomic analyses have deepened our understanding of aphid development and aphid-plant interactions. Our results have also provided useful insight into the molecular mechanisms underlying the biological variations in aphids, especially in adaptation to different host plants.