Original Research ARTICLE
Transcriptome Analysis of Green Peach Aphid (Myzus persicae): Insight into Developmental Regulation and Inter-Species Divergence
- 1Institute of Plant Protection, Jiangsu Academy of Agricultural Sciences, Nanjing, China
- 2Department of Entomology, Texas A&M University, College Station, TX, USA
- 3Ministry of Agriculture Key Laboratory of Agricultural Entomology, Institute of Insect Sciences, Zhejiang University, Hangzhou, China
- 4Department of Plant Pathology and Microbiology, Texas A&M University, College Station, TX, USA
- 5Department of Soil and Crop Sciences, Texas A&M University, College Station, TX, USA
- 6Biotechnology Research Institute, Chinese Academy of Agricultural Sciences, Beijing, China
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.
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 (Yang et al., 2012), 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 RNA-seq 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.
Materials and Methods
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 log2 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 CFX384TM 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., 2011, 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., 2011, 2014).
Results and Discussion
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 gap-filling, 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).
Functional Annotation and Classification of the Assembled Unigenes
Of the 62,627 unigenes, 33,543 were annotated by referencing to the Nr database (Table S2); 66.66% of the annotated sequences had very strong homology (E < 10−60), 12.02% showed strong homology (10−60 < E < 10−30) and the rest 21.32% showed homology (10−30 < E < 10−5) to known sequences. With respect to species, 92.30% of the unique annotated sequences matched to pea aphid, 1.45% to Tribolium castaneum, 0.49% to Bombus impatiens, and 0.41% to Camponotus floridana.
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).
Figure 1. Distribution of green peach aphid sequences by GO category. GO classification includes three domains: biological process, cellular component, and molecular function. The y-axis shows the number of matching unigenes in a category.
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 P450-mediated metabolic activity. In parallel, green peach aphid adults likely ingest more plant materials, thus more allelochemicals from host plants, necessitating higher detoxification capacity.
Figure 2. Fold change distribution of green peach aphid unigenes between nymphs and adults. The x-axis shows the fold change (log2 ratio) of gene expression in nymphs compared to adults. |Log2| values of 2244 unigenes are higher than 1, indicating potential importance during developmental transition.
Table 2. Significantly enriched insect-related KEGG pathways represented by the genes differentially expressed between nymphs and adults.
Table 3. Differentially expressed detoxification and cuticle formation-related genes in adult and nymph.
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).
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 −log10Q of the category. The GO term with highest −log10Q was determined the most significant enrichment.
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 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 encounter more types of toxic plant metabolites than pea aphids.
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 single-copy 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, 1999, 2001; Barbehenn 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).
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.
KZS, JF, and RJ conceived and designed the experiments. RJ performed the experiments. RJ, YW, YC, MZ, HZ, LZ, and KZS contributed to the transcriptome data analysis. RJ, JF, and KZS wrote the manuscript. All authors contributed to the discussion and approved the final manuscript.
This research was supported by the National Natural Science Foundation of China (31501636), the China Postdoctoral Science Foundation (2014M551529), Jiangsu Agricultural Science and Technology Independent Innovation Fund (CX(15)1055), Texas A&M Genomics Seed Grant, Texas Invasive Ant Research and Management Seed Grant, and the USDA-AFRI grant (2014-67013-21781).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We would like to thank Dr. Ron Salzman for his critical review of the manuscript. We appreciate technical assistance in sample preparation for sequencing from Dr. Charles Johnson and Dr. Richard Metz.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fpls.2016.01562
Amenya, D. A., Naguran, R., Lo, T.-C., Ranson, H., Spillings, B. L., Wood, O. R., et al. (2008). Over expression of a cytochrome P450 (CYP6P9) in a major African malaria vector, Anopheles funestus, resistant to pyrethroids. Insect Mol. Biol. 17, 19–25. doi: 10.1111/j.1365-2583.2008.00776.x
Arensburger, P., Megy, K., Waterhouse, R. M., Abrudan, J., Amedeo, P., Antelo, B., et al. (2010). Sequencing of Culex quinquefasciatus establishes a platform for mosquito comparative genomics. Science 330, 86–88. doi: 10.1126/science.1191864
Atamian, H. S., Chaudhary, R., Cin, V. D., Bao, E., Girke, T., and Kaloshian, I. (2013). In planta expression or delivery of potato aphid Macrosiphum euphorbiae effectors Me10 and Me23 enhances aphid fecundity. Mol. Plant Microbe Interact. 26, 67–74. doi: 10.1094/MPMI-06-12-0144-FI
Bansal, R., Mian, M., Mittapalli, O., and Michel, A. P. (2014). RNA-Seq reveals a xenobiotic stress response in the soybean aphid, Aphis glycines, when fed aphid-resistant soybean. BMC Genomics 15:972. doi: 10.1186/1471-2164-15-972
Barbehenn, R. V. (1999). Non-absorption of ingested lipophilic and amphiphilic allelochemicals by generalist grasshoppers: the role of extractive ultrafiltration by the peritrophic envelope. Arch. Insect Biochem. Physiol. 42, 130–137. doi: 10.1002/(SICI)1520-6327(199910)42:2<130::AID-ARCH3>3.0.CO;2-C
Barbehenn, R. V., and Stannard, J. (2004). Antioxidant defense of the midgut epithelium by the peritrophic envelope in caterpillars. J. Insect Physiol. 50, 783–790. doi: 10.1016/j.jinsphys.2004.05.012
Bass, C., Puinean, A. M., Andrews, M., Cutler, P., Daniels, M., Elias, J., et al. (2011). Mutation of a nicotinic acetylcholine receptor β subunit is associated with resistance to neonicotinoid insecticides in the aphid Myzus persicae. BMC Neurosci. 12:51. doi: 10.1186/1471-2202-12-51
Bhatia, V., Bhattacharya, R., Uniyal, P. L., Singh, R., and Niranjan, R. S. (2012). Host generated siRNAs attenuate expression of serine protease gene in Myzus persicae. PLoS ONE 7:e46343. doi: 10.1371/journal.pone.0046343
Bonasio, R., Zhang, G., Ye, C., Mutti, N. S., Fang, X., Qin, N., et al. (2010). Genomic comparison of the ants Camponotus floridanus and Harpegnathos saltator. Science 329, 1068–1071. doi: 10.1126/science.1192428
Cardona, C., Miles, J. W., Zu-iga, E., and Sotelo, G. (2010). Independence of resistance in Brachiaria spp. to nymphs or to adult spittlebugs (Hemiptera: Cercopidae): implications for breeding for resistance. J. Econ. Entomol. 103, 1860–1865. doi: 10.1603/EC10004
Cheng, T., Fu, B., Wu, Y., Long, R., Liu, C., and Xia, Q. (2015). Transcriptome sequencing and positive selected genes analysis of Bombyx mandarina. PLoS ONE 10:e0122837. doi: 10.1371/journal.pone.0122837
Chi, Y. H., Ahn, J.-E., Yun, D.-J., Lee, S. Y., Liu, T.-X., and Zhu-Salzman, K. (2011). Changes in oxygen and carbon dioxide environment alter gene expression of cowpea bruchids. J. Insect Physiol. 57, 220–230. doi: 10.1016/j.jinsphys.2010.11.011
Chinwalla, A. T., Cook, L. L., Delehaunty, K. D., Fewell, G. A., Fulton, L. A., Fulton, R. S., et al. (2002). Initial sequencing and comparative analysis of the mouse genome. Nature 420, 520–562. doi: 10.1038/nature01262
Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610
Devine, G. J., Harling, Z. K., Scarr, A. W., and Devonshire, A. L. (1996). Lethal and sublethal effects of imidacloprid on nicotine-tolerant Myzus nicotianae and Myzus persicae. Pestic. Sci. 48, 57–62. doi: 10.1002/(SICI)1096-9063(199609)48:1<57::AID-PS435>3.0.CO;2-9
Dombrovsky, A., Huet, H., Zhang, H., Chejanovsky, N., and Raccah, B. (2003). Comparison of newly isolated cuticular protein genes from six aphid species. Insect Biochem. Mol. Biol. 33, 709–715. doi: 10.1016/S0965-1748(03)00065-1
Elmer, K. R., Fan, S., Gunter, H. M., Jones, J. C., Boekhoff, S., Kuraku, S., et al. (2010). Rapid evolution and selection inferred from the transcriptomes of sympatric crater lake cichlid fishes. Mol. Ecol. 19, 197–211. doi: 10.1111/j.1365-294X.2009.04488.x
Elzinga, D. A., De Vos, M., and Jander, G. (2014). Suppression of plant defenses by a Myzus persicae (green peach aphid) salivary effector protein. Mol. Plant Microbe Interact. 27, 747–756. doi: 10.1094/MPMI-01-14-0018-R
Ewer, J., Gammie, S. C., and Truman, J. W. (1997). Control of insect ecdysis by a positive-feedback endocrine system: roles of eclosion hormone and ecdysis triggering hormone. J. Exp. Biol. 200, 869–881.
Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883
Harrison, T. L., Zangerl, A. R., Schuler, M. A., and Berenbaum, M. R. (2001). Developmental variation in cytochrome P450 expression in Papilio polyxenes in response to xanthotoxin, a hostplant allelochemical. Arch. Insect Biochem. Physiol. 48, 179–189. doi: 10.1002/arch.1070
He, Q., Jones, D. C., Li, W., Xie, F., Ma, J., Sun, R., et al. (2016). Genome-wide identification of R2R3-MYB genes and expression analyses during abiotic stress in Gossypium raimondii. Sci. Rep. 6:22980. doi: 10.1038/srep22980
Hegedus, D., Erlandson, M., Gillott, C., and Toprak, U. (2009). New insights into peritrophic matrix synthesis, architecture, and function. Annu. Rev. Entomol. 54, 285–302. doi: 10.1146/annurev.ento.54.110807.090559
Hiraishi, H., Terano, A., Ota, S., Mutoh, H., Sugimoto, T., Razandi, M., et al. (1991). Oxygen metabolites stimulate mucous glycoprotein secretion from cultured rat gastric mucous cells. Am. J. Physiol. Gastrointest. Liver Physiol. 261, 662–668.
Ji, R., Yu, H., Fu, Q., Chen, H., Ye, W., Li, S., et al. (2013). Comparative transcriptome analysis of salivary glands of two populations of rice brown planthopper, Nilaparvata lugens, that differ in virulence. PLoS ONE 8:e79612. doi: 10.1371/journal.pone.0079612
Jones, C. M., Daniels, M., Andrews, M., Slater, R., Lind, R. J., Gorman, K., et al. (2011). Age-specific expression of a P450 monooxygenase (CYP6CM1) correlates with neonicotinoid resistance in Bemisiatabaci. Pestic. Biochem. Physiol. 101, 53–58. doi: 10.1016/j.pestbp.2011.07.004
Kelleher, E. S., Swanson, W. J., and Markow, T. A. (2007). Gene duplication and adaptive evolution of digestive proteases in Drosophila arizonae female reproductive tracts. PLoS Genet. 3:e148. doi: 10.1371/journal.pgen.0030148
Kim, H., and Lee, S. (2008). A molecular phylogeny of the tribe Aphidini (Insecta: Hemiptera: Aphididae) based on the mitochondrial tRNA/COII, 12S/16S and the nuclear EF1α genes. Syst. Entomol. 33, 711–721. doi: 10.1111/j.1365-3113.2008.00440.x
Kirkness, E. F., Bafna, V., Halpern, A. L., Levy, S., Remington, K., Rusch, D. B., et al. (2003). The dog genome: survey sequencing and comparative analysis. Science 301, 1898–1903. doi: 10.1126/science.1086432
Leonova, I., and Slynko, N. (2004). Life stage variations in insecticidal susceptibility and detoxification capacity of the beet webworm, Pyrausta sticticalis L. (Lep., Pyralidae). J. Appl. Entomol. 128, 419–425. doi: 10.1111/j.1439-0418.2004.00866.x
Liao, C.-Y., Zhang, K., Niu, J.-Z., Ding, T.-B., Zhong, R., Xia, W.-K., et al. (2013). Identification and characterization of seven glutathione S-transferase genes from citrus red mite, Panonychus citri (McGregor). Int. J. Mol. Sci. 14, 24255–24270. doi: 10.3390/ijms141224255
Matsuo, T., Sugaya, S., Yasukawa, J., Aigaki, T., and Fuyama, Y. (2007). Odorant-binding proteins OBP57d and OBP57e affect taste perception and host-plant preference in Drosophila sechellia. PLoS Biol. 5:e118. doi: 10.1371/journal.pbio.0050118
Mu, H., Sun, J., Fang, L., Luan, T., Williams, G. A., Cheung, S. G., et al. (2015). Genetic basis of differential heat resistance between two species of congeneric freshwater snails: insights from quantitative proteomics and base substitution rate analysis. J. Proteome Res. 14, 4296–4308. doi: 10.1021/acs.jproteome.5b00462
Nauen, R., Bielza, P., Denholm, I., and Gorman, K. (2008). Age-specific expression of resistance to a neonicotinoid insecticide in the whitefly Bemisia tabaci. Pest Manag. Sci. 64, 1106–1110. doi: 10.1002/ps.1654
Nauen, R., Strobel, J., Tietjen, K., Otsu, Y., Erdelen, C., and Elbert, A. (1996). Aphicidal activity of imidacloprid against a tobacco feeding strain of Myzus persicae (Homoptera: Aphididae) from Japan closely related to Myzus nicotianae and highly resistant to carbamates and organophosphates. Bull. Entomol. Res. 86, 165–171. doi: 10.1017/S0007485300052408
Nikou, D., Ranson, H., and Hemingway, J. (2003). An adult-specific CYP6 P450 gene is overexpressed in a pyrethroid-resistant strain of the malaria vector, Anopheles gambiae. Gene 318, 91–102. doi: 10.1016/S0378-1119(03)00763-7
Ollivier, M., Legeai, F., and Rispe, C. (2010). Comparative analysis of the Acyrthosiphon pisum genome and expressed sequence tag-based gene sets from other aphid species. Insect Mol. Biol. 19, 33–45. doi: 10.1111/j.1365-2583.2009.00976.x
Pereira, R. J., Barreto, F. S., Pierce, N. T., Carneiro, M., and Burton, R. S. (2016). Transcriptome-wide patterns of divergence during allopatric evolution. Mol. Ecol. 25, 1478–1493. doi: 10.1111/mec.13579
Pertea, G., Huang, X., Liang, F., Antonescu, V., Sultana, R., Karamycheva, S., et al. (2003). TIGR Gene Indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics 19, 651–652. doi: 10.1093/bioinformatics/btg034
Pitino, M., and Hogenhout, S. A. (2013). Aphid protein effectors promote aphid colonization in a plant species-specific manner. Mol. Plant Microbe Interact. 26, 130–139. doi: 10.1094/MPMI-07-12-0172-FI
Ramsey, J. S., Rider, D. S., Walsh, T. K., De Vos, M., Gordon, K. H. J., Ponnala, L., et al. (2010). Comparative analysis of detoxification enzymes in Acyrthosiphon pisum and Myzus persicae. Insect Mol. Biol. 19, 155–164. doi: 10.1111/j.1365-2583.2009.00973.x
Ramsey, J. S., Wilson, A. C., De Vos, M., Sun, Q., Tamborindeguy, C., Winfield, A., et al. (2007). Genomic resources for Myzus persicae: EST sequencing, SNP identification, and microarray design. BMC Genomics 8:423. doi: 10.1186/1471-2164-8-423
Ren, L., Tan, X.-J., Xiong, Y.-F., Xu, K., Zhou, Y., Zhong, H., et al. (2014). Transcriptome analysis reveals positive selection on the divergent between topmouth culter and zebrafish. Gene 552, 265–271. doi: 10.1016/j.gene.2014.09.053
Robbins, W., Kaplanis, J., Thompson, M., Shortino, T., Cohen, C., and Joyner, S. (1968). Ecdysones and analogs: effects on development and reproduction of insects. Science 161, 1158–1160. doi: 10.1126/science.161.3846.1158
Strode, C., Steen, K., Ortelli, F., and Ranson, H. (2006). Differential expression of the detoxification genes in the different life stages of the malaria vector Anopheles gambiae. Insect Mol. Biol. 15, 523–530. doi: 10.1111/j.1365-2583.2006.00667.x
Su, Y.-L., Li, J.-M., Li, M., Luan, J.-B., Ye, X.-D., Wang, X.-W., et al. (2012). Transcriptomic analysis of the salivary glands of an invasive whitefly. PLoS ONE 7:e39303. doi: 10.1371/journal.pone.0039303
Sun, Y. F., De Biasio, F., Qiao, H. L., Iovinella, I., Yang, S. X., Ling, Y., et al. (2012). Two odorant-binding proteins mediate the behavioural response of aphids to the alarm pheromone (E)-ß-farnesene and structural analogues. PLoS ONE 7:e32759. doi: 10.1371/journal.pone.0032759
Swanson, W. J., Wong, A., Wolfner, M. F., and Aquadro, C. F. (2004). Evolutionary expressed sequence tag analysis of Drosophila female reproductive tracts identifies genes subjected to positive selection. Genetics 168, 1457–1465. doi: 10.1534/genetics.104.030478
Tian, C., Tay, W. T., Feng, H., Wang, Y., Hu, Y., and Li, G. (2015). Characterization of Adelphocoris suturalis (Hemiptera: Miridae) transcriptome from different developmental stages. Sci. Rep. 5:11042. doi: 10.1038/srep11042
Tiffin, P., and Hahn, M. W. (2002). Coding sequence divergence between two closely related plant species: Arabidopsis thaliana and Brassica rapa ssp. pekinensis. J. Mol. Evol. 54, 746–753. doi: 10.1007/s00239-001-0074-1
von Dohlen, C. D., Rowe, C. A., and Heie, O. E. (2006). A test of morphological hypotheses for tribal and subtribal relationships of Aphidinae (Insecta: Hemiptera: Aphididae) using DNA sequences. Mol. Phylogenet. Evol. 38, 316–329. doi: 10.1016/j.ympev.2005.04.035
Vyas, M., Fisher, T. W., He, R., Nelson, W., Yin, G., Cicero, J. M., et al. (2015). Asian citrus psyllid expression profiles suggest candidatus liberibacter asiaticus-mediated alteration of adult nutrition and metabolism, and of nymphal development and immunity. PLoS ONE 10:e0130328. doi: 10.1371/journal.pone.0130328
Wang, D., Liu, Q., Jones, H. D., Bruce, T., and Xia, L. (2014). Comparative transcriptomic analyses revealed divergences of two agriculturally important aphid species. BMC Genomics 15:1023. doi: 10.1186/1471-2164-15-1023
Wang, R.-L., Zhu-Salzman, K., Baerson, S. R., Xin, X.-W., Li, J., Su, Y.-J., et al. (2016). Identification of a novel cytochrome P450 CYP321B1 gene from tobacco cutworm (Spodoptera litura) and RNA interference to evaluate its role in commonly used insecticides. Insect Sci. doi: 10.1111/1744-7917.12315. [Epub ahead of print].
Wang, X.-W., Luan, J.-B., Li, J.-M., Su, Y.-L., Xia, J., and Liu, S.-S. (2011). Transcriptome analysis and comparison reveal divergence between two invasive whitefly cryptic species. BMC genomics 12:458. doi: 10.1186/1471-2164-12-458
Werren, J. H., Richards, S., Desjardins, C. A., Niehuis, O., Gadau, J., Colbourne, J. K., et al. (2010). Functional and evolutionary insights from the genomes of three parasitoid Nasonia species. Science 327, 343–348. doi: 10.1126/science.1178028
Yang, N., Xie, W., Jones, C. M., Bass, C., Jiao, X., Yang, X., et al. (2013). Transcriptome profiling of the whitefly Bemisia tabaci reveals stage-specific gene expression signatures for thiamethoxam resistance. Insect Mol. Biol. 22, 485–496. doi: 10.1111/imb.12038
Yang, W., Qi, Y., Bi, K., and Fu, J. (2012). Toward understanding the genetic basis of adaptation to high-elevation life in poikilothermic species: a comparative transcriptomic analysis of two ranid frogs, Rana chensinensis and R. kukunoris. BMC genomics 13:588. doi: 10.1186/1471-2164-13-588
Yang, Z., and Nielsen, R. (2000). Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol. Biol. Evol. 17, 32–43. doi: 10.1093/oxfordjournals.molbev.a026236
Yu, H., Ji, R., Ye, W., Chen, H., Lai, W., Fu, Q., et al. (2014). Transcriptome analysis of fat bodies from two brown planthopper (Nilaparvata lugens) populations with different virulence levels in rice. PLoS ONE 9:e88528. doi: 10.1371/journal.pone.0088528
Zhang, K., Niu, J. Z., Ding, T. B., Dou, W., and Wang, J. J. (2013). Molecular characterization of two carboxylesterase genes of the citrus red mite, Panonychus citri (Acari: Tetranychidae). Arch. Insect Biochem. Physiol. 82, 213–226. doi: 10.1002/arch.21087
Zhang, L., Yan, H.-F., Wu, W., Yu, H., and Ge, X.-J. (2013). Comparative transcriptome analysis and marker development of two closely related Primrose species (Primula poissonii and Primula wilsonii). BMC Genomics 14:329. doi: 10.1186/1471-2164-14-329
Keywords: Myzus persicae, Acyrthosiphon pisum, nymph and adult, transcriptome, developmental regulation, synonymous and nonsynonymous substitutions, host plant adaptation
Citation: Ji R, Wang Y, Cheng Y, Zhang M, Zhang H-B, Zhu L, Fang J and Zhu-Salzman K (2016) Transcriptome Analysis of Green Peach Aphid (Myzus persicae): Insight into Developmental Regulation and Inter-Species Divergence. Front. Plant Sci. 7:1562. doi: 10.3389/fpls.2016.01562
Received: 28 April 2016; Accepted: 04 October 2016;
Published: 21 October 2016.
Edited by:Jyoti Shah, University of North Texas, USA
Reviewed by:Muthu Venkateshwaran, University of Wisconsin–Platteville, USA
Vamsi J. Nalam, Indiana University–Purdue University Fort Wayne, USA
Copyright © 2016 Ji, Wang, Cheng, Zhang, Zhang, Zhu, Fang and Zhu-Salzman. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.