Reduction in fecundity and shifts in cellular processes by a native virus on an invasive insect.

Pathogens and their vectors have coevolutionary histories that are intricately intertwined with their ecologies, environments, and genetic interactions. The soybean aphid, Aphis glycines, is native to East Asia but has quickly become one of the most important aphid pests in soybean-growing regions of North America. In this study, we used bioassays to examine the effects of feeding on soybean infected with a virus it vectors (Soybean mosaic virus [SMV]) and a virus it does not vector (Bean pod mottle virus [BPMV]) have on A. glycines survival and fecundity. The genetic underpinnings of the observed changes in fitness phenotype were explored using RNA-Seq. Aphids fed on SMV-infected soybean had transcriptome and fitness profiles that were similar to that of aphids fed on healthy control plants. Strikingly, a significant reduction in fecundity was seen in aphids fed on BPMV-infected soybean, concurrent with a large and persistent downregulation of A. glycines transcripts involved in regular cellular activities. Although molecular signatures suggested a small regulatory RNA pathway defense response was repressed in aphids feeding on infected plants, BPMV did not appear to be replicating in the vector. These results suggest that incompatibilities with BPMV or the effects of BPMV infection on soybean caused A. glycines to allot available energy resources to survival rather than reproduction and other core cellular processes. Ultimately, the detrimental impacts to A. glycines may reflect the short tritrophic evolutionary histories between the insect, plant, and virus.

Introduction

Ecological interactions among viruses, the plants they infect, and the insects they are transmitted by are evolutionarily intertwined and complex (Froissart et al. 2010; Hodge and Powell 2010; Mauck et al. 2010). Viruses can affect insect vector life history directly or indirectly, with positive or negative consequences to fecundity, longevity, and other fitness traits (Rubinstein and Czosnek 1997; Jimenez-Martinez et al. 2004; Maris et al. 2004; Belliure et al. 2005; Colvin et al. 2006; Stout et al. 2006). Moreover, viruses can even alter the life history of insects that do not transmit them but which feed on infected hosts (Tu et al. 2013; Xu et al. 2013). To date, little is known about the underlying molecular mechanisms associated with virus-induced changes to insect fitness. A detailed understanding of the ecological and genetic bases of insect–virus interactions may uncover novel and specific molecular targets for disease control.

Soybean mosaic virus (SMV) and Bean pod mottle virus (BPMV) are two of the most destructive viruses of soybean (Glycine max (L.) Merr.), resulting in significant yield losses and reduced seed quality. The historical range of the comovirus, BPMV, is not well studied, but it was first found in Eastern and Midwestern North America (Zaumeyer and Thomas 1948; Ghabrial et al. 1977; Mabry et al. 2003) and may have originated on beans (Phaseolus spp.). Ancestrally, this plant host, as well as the primary BPMV vector, the bean leaf beetle (Cerotoma trifurcata), were distributed within the Americas. Only recently has BPMV been found in other continents (Shahraeen et al. 2005), and it has yet to be documented in East Asia. In contrast, the potyvirus, SMV, is present in most countries where soybean is grown (Clinton 1916; Gardner and Kendrick 1921; Cui et al. 2011). SMV has been reported in East Asia, where soybean is native, since the early 1980s, but was likely present much longer (reviewed in Li et al. 2010). Consequently, SMV shares longer coevolutionary histories with soybean and its native pests, including the soybean aphid, Aphis glycines Matsumura (Wang and Ghabrial 2002; Wu et al. 2004), than does BPMV.

Native to East Asia, A. glycines was first detected in Wisconsin in 2000 and has since spread throughout much of the North Central US and Eastern Canada (Venette and Ragsdale 2004; Wu et al. 2004). The aphid’s population growth and vector capability have quickly made it one of the most important arthropod pests of soybean in North America (Ragsdale et al. 2007). As with all aphid-borne potyviruses, A. glycines transmits SMV in a nonpersistent manner (Hogenhout et al. 2008; Cui et al. 2011). Nonpersistently transmitted viruses do not breach the gut barrier or infect the insect; rather they are retained in the insect stylet or foregut prior to transmission. Aphids can acquire the virus in a single probe of an infected leaf and subsequently transmit the virus for only a few minutes to a few hours (Hooks and Fereres 2006). BPMV is transmitted by beetles in the family Chrysomelidae and is not reported to be experimentally or naturally transmitted by A. glycines or any other aphid species.

Virus-associated molecular, cellular, and physiological changes are well documented in infected host plants (Maule 2007). However, the responses of insect vectors to infected plants have been examined far less and have mostly focused on persistently transmitted viruses. A handful of studies have investigated transcriptome changes in aphids and other hemipteran vectors in response to feeding on virus-infected host plants (Brault et al. 2010; Luan et al. 2011; Gotz et al. 2012; Xu et al. 2012; Cassone et al. 2014). Other studies have examined changes in insect vector fitness in response to virus exposure (Rubinstein and Czosnek 1997; McKenzie 2002; Jiu et al. 2007; Mann et al. 2008; Sidhu et al. 2009). To date, no study has attempted to characterize any associations between insect vector fitness phenotypes and genetic responses to plant virus exposure. Because BPMV is not transmitted by the soybean aphid and the nonpersistently transmitted SMV interacts only briefly with the insect and does not cross vectormembranes, we hypothesized that the fitness and molecular responses of the vector to these viruses will be limited.

Most plant virus research has focused on aphid vector transmission dynamics and molecular biology, including the documentation of over 150 virus–aphid associations (Hogenhout et al. 2008), and the sequencing of the pea aphid (Acyrthosiphon pisum) genome (Richards et al. 2010). In this study, we carried out fitness bioassays to measure soybean aphid fecundity and survival responses (if any) to feeding on SMV- and BPMV-infected soybean. Next-generation sequencing and quantitative real-time polymerase chain reaction (PCR) were carried out to characterize the transcriptional profiles for the different treatments. SMV and BPMV were chosen to examine divergent responses of the aphid to viruses for which it is a vector and nonvector, and early (4 h) and late (7 days) time points were selected to characterize both the immediate and residual impacts to the insects feeding on virus-infected plants. Unexpectedly, aphids feeding on BPMV-infected soybean showed significant reductions in fecundity coupled with substantial downregulation of transcripts involved in regular cellular activities.

Materials and Methods

Aphid Colony and Virus Maintenance

Experiments were carried out using a laboratory colony of A. glycines biotype 3 established from field collections in Indiana and maintained in a growth chamber at the Ohio Agricultural Research and Development Center, Wooster, OH (Hill et al. 2010). Aphids were maintained on “Sloan” seedlings placed in 15 cm × 7.5 cm × 15 cm cages in growth chambers under controlled conditions of 25 °C and 75% RH with a 16 h:8 h light–dark cycle.

The BPMV isolate was maintained in Sloan soybean through serial leaf-rub inoculation with inoculum made from leaves of infected plants (Louie et al. 2000). SMV was maintained in soybean through serial transmission by A. glycines. To generate experimental plants, inoculum made by grinding infected leaf tissue into 10 mM KHPO4, pH 7 (1:4 w/v), was mechanically inoculated onto 1.5-week-old soybean plants. After visual assessment of symptom development, virus presence was confirmed and titer estimated in using enzyme-linked immunosorbent assays (Todd et al. 2010).

Fitness Studies

Experimental Design

Thirty apterous 2-day-old adult A. glycines were randomly chosen from maintenance cages and placed on the youngest trifoliate leaf of a V2-stage SMV-infected, BPMV-infected, or healthy soybean plant (i.e., experimental treatments). Aphids were restricted to the trifoliate leaf using 4.5 cm × 1.5 cm × 5.5 cm dacron cages secured with padded edges to avoid damaging the plant. After 7 days, the survival rate and gross fecundity were determined for each treatment/replicate combination. Survival rate was defined as the number of living adult aphids at the end of the experiment/the number of adult aphids at start of the experiment and gross fecundity was determined as the total number of offspring present. Independent sample t-tests (two tailed, P < 0.05) were carried out to examine differences in survival rate and fecundity between treatments. Each treatment was replicated ten times using adults from different cohorts and multiple maintenance cages.

Transcriptome Analysis

A total of 120 apterous aphids (<24 h) were randomly selected in equal numbers from two maintenance cages and starved for 3 h. For each treatment/time combination, 20 aphids were caged on the first trifoliate leaf of 3–3.5-week-old SMV-infected, BPMV-infected, or healthy (control) soybean plants as described above. After 4 h or 7 days, the aphids from each cage were collected, frozen in liquid nitrogen, and stored at −80 °C until RNA isolation. Each treatment/time point was replicated four times using aphids from different cohorts (24 samples total) (table 1).

Table 1

Experimental Design for RNA-Seq Transcriptional Profiling of Aphis glycines Biotype 3 Using the Illumina HiSeq Platform

Species	Treatment	Time Points
		4 h	7 days
A. glycines♀	SMV	4a	4
	BPMV	4	4
	Healthy control	4	4

aNumber of replicates for each treatment/time point combination.

RNA Isolation and cDNA Library Synthesis

Total RNA was extracted from the pools of 20 aphids using the RNeasy Mini Kit (QIAGEN, Germantown, MD). RNA quantity and quality was assessed using the Experion Automated Electrophoresis System (Bio-Rad, Hercules, CA). All samples selected for cDNA synthesis had RNA Quality Indicator values of 10.

RNA (1 µg/sample) was used to generate adaptor-ligated double-stranded cDNA libraries for RNA-Seq using the TruSeq Sample Prep Kit V1 and V2 (Illumina, San Diego, CA) following the manufacturer’s protocol. Quantification of ds-cDNA was done using the Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA) and quality assessed using the Experion Automated Electrophoresis System (Bio-Rad, Hercules, CA). Samples were diluted to 17.5 nM and pooled to generate the multiplexed cDNA library (24 adaptor-tagged pooled samples total: 3 treatments × 2 time points × 4 replicates).

Illumina Sequencing and Transcriptome Assembly

The cDNA library (3.6 × 10−5 nmol) was sequenced on one flow cell lane using the Illumina HiSeq 2000 platform at the Ohio State University Comprehensive Cancer Center. Four fluorescently labeled nucleotides and a specialized polymerase were used to determine the clusters base by base in parallel. The mean library insert sequence length was 274 bp, and both ends of the library were sequenced to generate 100 nt raw paired-end reads. Illumina Analysis Package CASAVA 1.8.2 was used to perform bcl conversion and demultiplexing. Image deconvolution and quality value calculations were carried out using the Illumina GA pipeline v1.6.

Raw reads were imported into CLC Genomics Workbench (v6.0.1, CLC Bio) and trimmed for quality, adapter indexes, and poly (A) tails using the default settings (ambiguous limit = 2, quality limit = 0.05). Redundant reads were removed using the duplicate removal plugin in CLC. Processed reads were de novo assembled into scaffolded contigs of at least 250 nt using an algorithm based on de Bruijn graphs and the optimized parameters defined in supplementary table S1, Supplementary Material online. Using custom R scripts (available on request from BJC), the data set was filtered to contain only contigs that had a minimum of five mapped reads for any three replicates in at least one treatment/time combination. The set of transcripts was obtained by collapsing reads with ≥90% sequence similarity into clusters and retrieving the longest contig using CD-HIT-EST (Li and Godzik 2006). The raw sequence reads can be retrieved from the National Center for Biotechnology Information (NCBI) short sequence read archive under the accession number SRP031835.

Functional Annotation and Ortholog Comparison

De novo assembled transcripts were assigned hierarchical gene ontologies (GOs) terms on the basis of biological processes, molecular functions, and cellular components using the platform-independent Java 6 implementation of the BLAST2GO software (Gotz et al. 2008). The top ten BLASTx hits with a cut-off E value of 10−3, and similarity cutoff of 55% were considered for GO annotation.

Pair-wise comparisons of transcripts to cDNA databases of four insect species were carried out using desktop downloaded tBLASTx software and an E value < 10−10. The ortholog databases were retrieved on June 18, 2012, from the ftp files of NCBI or Ensembl: Ac. pisum (pea aphid, order Hemiptera, 37,994 sequences), Apis mellifera (honey bee, order Hymenoptera, 18,542 sequences), Tribolium castaneum (red flour beetle, order Coleoptera, 14,366 sequences), and Drosophila melanogaster (fruit fly, order Diptera, 19,233 sequences).

Transcriptome Analysis

Using the CLC Bio Transcriptomics Analysis tool, de novo assembled transcript counts for paired reads were normalized by calculating the number of unique (i.e., unambiguous) reads per kilobase of exon model per million mapped reads (RPKM) (Mortazavi et al. 2008). Quality assessment of the normalized data file was conducted using the quality control function in CLC Bio. Heat maps for transcriptional profiles were generated using heatmaps.2 in the gplots package in R.

Bayes-moderated t-tests were carried out on the normalized data set to detect virus-responsive transcripts at 4 h and 7 days. One-way analysis of variance (ANOVA) was used to identify transcripts that were differentially expressed among treatments at each time point. Significance was defined at a P value < 0.05 (false discovery rate < 0.20) with a minimum five RPKM-change difference for all analyses unless otherwise specified. For comparison of statistical methods, differentially expressed transcripts were also identified independently using the R module edgeR (Robinson and Oshlack 2010), which uses an overdispersed Poisson model to moderate the dispersion.

The DAVID v6.7 annotation clustering module (Huang et al. 2009a, 2009b) was used to classify differentially expressed transcripts into functional groups. Transcripts were first converted to D. melanogaster transcript IDs, then enrichment of GO and other annotation terms in candidate sublists were explored using the functional annotation clustering tool. The enrichment score ranks the biological significance of gene groups based on overall EASE scores (modified Fisher’s exact test) of all enriched annotation terms, thereby accounting for the relative importance of the groups as part of a probing rather than strictly statistical analysis. Significant clusters were defined using the following parameters: similarity term overlap = 6; similarity threshold = 0.6; initial group membership = 5; final group membership = 10; multiple linkage threshold = 0.15; EASE = 1.8 (equivalent to P ≤ 0.01).

Transcript Accumulation Profiling Using Quantitative Reverse Transcription PCR (RT-qPCR)

Each targeted transcript was validated using two treatments (BPMV-infected and healthy soybean) and time points (4 h and 7 days). To better gauge the likelihood that the downregulation of transcripts was a laboratory colony-specific phenomenon, transcripts were validated from two independent colonies of A. glycines from different geographic locales: the experimental biotype 3 colony and a biotype 1 colony established from multiple field collections in Ohio. Aphids were collected and stored identically to the RNA-Seq experiment. Each treatment/time point was replicated three times using adults from different cohorts. Total RNA was extracted from pools of 20 aphids, quantified, and quality assessed using the protocols outlined earlier. RNA was treated with DNase I (Invitrogen, Carlsbad, CA), and cDNA was transcribed using the iScript synthesis kit (Bio-Rad, Hercules, CA). Primer pairs for targeted mRNAs were designed using Primer3 (Rozen and Skaletsky 2000) and are shown in supplementary table S2, Supplementary Material online. In preliminary experiments, primer concentrations of 50 nM, 300 nM, and 900 nM were tested to determine optimal qPCR conditions for each gene. A TATA-box binding protein gene was used as the endogenous qPCR control with forward and reverse primers pairs F-GGGTTTTGAATAGTTTGTA and R-GCTACTCCACATAGTATG (Bansal et al. 2012).

The qPCR reactions (15 µl) were performed in duplicate using IQ 2X SYBR Green Supermix (Bio-Rad, Hercules, CA) and 300 nM of each primer on a CFX96 Real-Time PCR Detection System according to the manufacturer’s recommendations. Cycling conditions were 50 °C for 2 min, 95 °C for 10 min, 41 cycles of denaturation at 95 °C for 15 s and annealing/extension at 60 °C for 1 min, and 95 °C for 15 s. PCR efficiency (E) was evaluated by performing a dilution series experiment using a target assay and the equation E = 10(−1/slope) (Pfaffl 2001).

Expression levels were measured separately for reference and target genes, using three biological replicates for each treatment/time point combination. Relative transcript abundance was calculated using the 2−ΔΔT method (Livak and Schmittgen 2001). Threshold cycle (CT) values reported by the CFX96 Real-Time PCR Detection System were normalized to the reference gene and converted to relative log2-fold differences between treatments for each time point. One-tailed independent t-tests (P < 0.05) were used for statistical determination of differential expression between each treatment/time point combination.

Detection of BPMV in Infected Plants

Our results (below) suggested that BPMV may be replicating in A. glycines. To test this, aphids were given a 7-day acquisition access period (AAP) on BPMV-infected plants followed by a 7 h starvation and 1.5 h feeding period on healthy soybean to clear virus from the insect gut. The 7 days AAP was chosen to test for BPMV infection because it represented the longest exposure of aphids to BPMV in our time series. Total RNA was harvested from pools of 15 aphids, quantified, and quality assessed using the protocols outlined earlier. RNA from aphids fed on healthy and BPMV-infected plants without the clearance period served as the negative and positive controls, respectively. RNA was converted to cDNA using the iScript synthesis kit (Bio-Rad, Hercules, CA). Primer pairs targeted the second segment of the BPMV genome were designed using Primer3 (Rozen and Skaletsky 2000) as follows: BPMV2: F-GAAGAAGGGCCCAACAACATGTCTATCTCTCAGCAGACCGTTTGGA and R-GAAGTATGCGGCCGCTGCAGAAGATTCTGCATTTTCAGAAGA. PCR reactions (15 µl) were performed using GoTaq Green Master Mix (Promega, Madison, WI) and the following cycling conditions: 94 °C for 2 min; 35 cycles of 94 °C for 30 s, 58 °C for 30 s, and 72 °C for 30 s; then, 72 °C for 10 min. BPMV presence/absence was assessed by agarose gel electrophoresis.

Results

Fecundity and Survival of Aphids on Virus-Infected and Healthy Soybean

Differences in survival rates and fecundity were assessed among A. glycines fed on SMV-infected, BPMV-infected, and healthy soybean for 7 days. A total of 300 aphids were bioassayed per treatment across ten biological replicates harvested from different cohorts (fig. 1). Fecundity was significantly lower in aphids fed on BPMV-infected soybean (P < 0.001) compared with aphids in the other treatments. No dead offspring were found on the enclosed trifoliates after 7 days for any treatment. Reduced fecundity in the BPMV treatment did not appear due to a smaller pool of aphids available to produce offspring, as adult survival was not significantly different between treatments (P > 0.05). Across all replicates, aphids feeding on BPMV-infected soybean produced about half as many viable offspring after 7 days as those feeding on SMV-infected or healthy plants.

F

Survival (A) and fecundity (B) of Aphis glycines adult after feeding on BPMV-infected (white), SMV-infected (dark gray), or healthy control (light gray) soybean for 7 days. Error bars indicate standard deviation variation between ten replicates per treatment/fitness parameter.

Transcriptome Assembly

Although transcript sequence data were becoming increasingly available for A. glycines (Bai et al. 2010; Liu et al. 2012), a well-characterized transcriptome is not currently available. A cDNA library composed of 24 pooled samples (table 1) was sequenced, which generated 343,192,082 paired end reads of 100 nt. After trimming (quality, adapters, and poly (A) sequences) and duplicate removal, 7.2–18.9 million reads per sample were obtained. The A. glycines transcriptome was assembled de novo into 28,838 cDNA transcripts of ≥234 nt with a mean length of 872 nt. The assembly was similar to that of a previous A. glycines de novo assembly using the Rnnotator automated pipeline (Martin et al. 2010) in Galaxy (Goecks et al. 2010) at the OARDC, MCIC (Bansal R, unpublished results). A total of 364 transcripts were removed from the data set because they were of prokaryote or virus origin (40% from Buchnera aphidicola), including transcripts of 3.7 kb and 6 kb that matched sequences for both molecules of the BPMV genome (Macfarlane et al. 1991; Di et al. 1999) and a 9.6 kb transcript that matched the SMV genome sequence (Hill and Benner 1980). The abundance of the SMV and BPMV transcripts was assessed in aphids at both time points in their respective treatments (table 2). In both cases, the number of transcripts increased between 4 h to 7 days, likely due to the presence of residual virion in aphids continuously exposed to infected plants. An additional 399 transcripts were eliminated from the data set because they most closely matched soybean transcripts. The remaining 28,075 transcripts formed the data set for subsequent analyses.

Table 2

Average RPKM Expression Values of the Monopartite SMV and Bipartite BPMV (RNA 1 and BPMV RNA 2) Genome Segments in Aphis glycines for Each Treatment/Time Point Combination

Genome Seg.	SMV BPMV Control
	4 h	7 days	4 h	7 days	4 h	7 days
SMV	7.771	340.1	0.253	0.126	0.123	0.061
BPMV RNA 1	0.661	1.307	374.4	1152	0.307	1.339
BPMV RNA 2	1.886	4.27	668.8	3663	1.653	3.656

Functional Annotation and Pairwise Comparisons with Other Insect Species

BLAST2GO was used to explore the hierarchical associations between A. glycines transcripts by assigning GO and enzyme ontology (EC) during the BLAST, mapping, and annotation steps. A total of 12,909 transcripts had a significant BLASTx hit to the nr database (E value <10−3), of which 91% (n = 11,757) had the highest identity to pea aphid (Ac. pisum) transcripts. Of the transcripts with a BLAST match, 6,657 mapped to GO and/or EC terms in the nr database, and 4,205 could be annotated by biological process, molecular function, and/or cellular component (supplementary table S3, Supplementary Material online).

Pairwise sequence comparisons were carried out between the A. glycines transcriptome and the transcriptomes of four other insect species (fig. 2). Nearly 30% (n = 8,187) of transcripts had a significant match to all four insect databases, and roughly half (n = 14,210) had a match in at least one of the four databases (E value < 10−10). As expected, the majority (>98%) of transcripts that matched only one of the species corresponded with the related aphid, Ac. pisum.

F

Comparison of transcriptomes for Aphis glycines and four insect species. The Venn diagram shows pair-wise ortholog matches with tBLASTx (E values <10−10) for Apis mellifera (18,542 transcripts), Drosophila melanogaster (19,233 transcripts), Acyrthosiphon pisum (37,994 transcripts), and Tribolium castaneum (14,366 transcripts).

Transcriptional Changes in A. glycines Fed on Virus-Infected Soybean

RPKM-normalized transcript expression levels were examined to determine whether feeding on virus-infected soybean elicited changes in A. glycines that differed from those fed on healthy soybean. Quality control parameters indicated that the normalized sample values were normally distributed and homogenous. Bayes-moderated t-tests (P < 0.05) were implemented for each time point to compare the transcriptomes of aphids fed on SMV- and BPMV-infected plants to the healthy control. At 4 h, only 57 differentially expressed transcripts were shared between the SMV and BPMV fed aphids, with approximately half downregulated in the aphids fed on virus-infected soybean relative to the control. Considerably more differentially expressed transcripts were identified at 7 days (n = 1,075), and the majority of these were downregulated relative to insects fed on healthy plants (74%). Most of the differentially expressed transcripts at both time points (4 h: 47 [83%]; 7 days: 1,002 [93%]) were also identified using the edgeR package, indicating strong correlation between the statistical methods.

Functional annotation of the A. glycines transcriptome is in its preliminary stages. For this reason, we did not carry out any formal quantitative analyses of differentially expressed transcripts by function. Rather, we opted for a probing approach that relied on ortholog prediction using the well-characterized D. melanogaster transcriptome (BDGP5) assembly. Drosophila prediction was particularly informative for assigning function to differentially expressed transcripts. Although less than one-third of the A. glycines transcriptome matched a Drosophila ID, the majority of differentially expressed transcripts (below) could be assigned a putative function (74% of SMV and 63% for BPMV).

The DAVID annotation clustering module (Huang et al. 2009a, 2009b) was used to classify the differentially expressed A. glycines transcripts into functional clusters. Transcripts were first assigned their ortholog Drosophila ID (where possible) then divided into four lists: up or downregulated in aphids fed virus-infected soybean for 4 h or 7 days. No clusters with significant enrichment scores (i.e., >1.8 in minus logscale) were identified in the list containing transcripts downregulated relative to the healthy control at 4 h, whereas two clusters were identified in the transcript list comprised of upregulated transcripts at this time point (supplementary fig. S1, Supplementary Material online). The upregulated clusters contained transcripts involved in translation and cytoskeleton organization. A total of eight clusters were identified among the sets of differentially expressed transcripts at 7 days. The upregulated clusters comprised transcripts that function in protein translation, cytoskeleton organization, cell–cell signaling, and ion channel activity. Downregulated clusters are composed of transcripts associated with protein localization, cell cycle, and helicase activity. A full list of the annotation clusters with significant enrichment scores are displayed in supplementary figure S1, Supplementary Material online.

Differential Responses of A. glycines Fed on SMV- and BPMV-Infected Soybean

Although feeding on virus-infected soybean for 7 days brought about changes in transcript accumulation relative to those fed on healthy soybean, SMV- and BPMV-specific changes in transcript abundance also occurred. One-way ANOVA was used to detect transcripts with a significant treatment effect at 4 h and 7 days in the subset of A. glycines transcripts that did not elicit a common virus response described earlier. A total of 4,143 and 2,095 transcripts were differentially expressed in a virus-specific manner at 4 h and 7 days, respectively. A heat map of the 4,143 differentially expressed transcripts at 4 h indicated that nearly all of the differential accumulation was derived from the BPMV treatment (n = 4,116) and that 87% of these were downregulated relative to aphids fed on healthy or SMV-infected soybean plants (fig. 3). Similar in trend though not as dramatic, the majority of differential transcript expression at 7 days was from the aphids fed on BPMV-infected plants (n = 1,666), in which 70% (n = 1,171) transcripts were downregulated. For a subset of the differentially expressed transcripts at each time point, the response was opposite in the BPMV and SMV treatments. This was the case for 6% (n = 256) and 11% (n = 225) of differentially expressed transcripts at 4 h and 7 days, with the majority upregulated in the BPMV treatment and downregulated in the SMV treatment (4 h: n = 153 [60%]; 7 days: n = 176 [78%]).

F

Heat map of transcriptional profiles for Aphis glycines transcripts with a significant treatment effect at 4 h. Analyses of 4,143 transcripts that were differentially expressed based on one-way ANOVA with P-level significance at 0.05 and RPKM change >5. Each row represents an individual transcript; each column labeled 1–4 represents a replicate sample for insects fed on BPMV-infected, SMV-infected, or healthy (control) soybean. Within a row, the relative expression level is represented by a color that reflects the row z score (shown in the red–green key), calculated by subtracting the mean expression value for the row from individual sample values and dividing by the standard deviation of the row.

To examine the biological responses associated with virus-specific transcript accumulation, the differentially expressed transcripts were partitioned into four lists for each time point (up/downregulated in A. glycines fed on SMV-/BPMV-infected soybean relative to the alternative treatment and healthy control), and functional annotation was carried out using DAVID. Significantly enriched annotation clusters were identified in three of four lists for 4 h and 7 days, with between 1 and 11 clusters per treatment/time point combination (fig. 4). For aphids fed on BPMV-infected soybean, the clusters were most notably characterized by a downregulation of transcripts functioning in core cellular processes (e.g., transcription, intracellular transport, and helicase activity) and upregulation of transcripts encoding translation and energy production proteins. A single cluster of downregulated cellular respiration transcripts was identified for aphids fed on SMV-infected soybean.

F

Functional annotation clusters of Aphis glycines virus-specific transcripts in response to feeding on BPMV-infected, SMV-infected, or healthy (control) soybean. Drosophila orthologs were identified, and clusters of functionally related transcripts were identified using DAVID (Huang et al. 2009a, 2009b). The numbers of significant transcripts in each cluster are indicated above the columns, and DAVID enrichment scores are displayed inside each column.

Persistent Downregulation of Transcripts in A. glycines Fed on BPMV-Infected Soybean

Two-thirds (n = 1,815) of the 2,741 differentially expressed transcripts (virus-specific and virus-shared combined) in aphids fed on BPMV-infected soybean for 7 days were also differentially expressed and in the same direction at 4 h (supplementary fig. S2, Supplementary Material online). Most of these transcripts were downregulated (89%). Consistent with the 4 h functional annotation, the constitutively downregulated transcripts were mostly implicated in regular cellular activities, such as transport, transcription, and cell cycle (nine clusters in total). Considerably fewer transcripts (n = 83) were persistently differentially expressed in the aphids fed on SMV-infected soybean. Only one cluster comprising cellular respiration transcripts was downregulated in the SMV treatment.

Immune-Related Responses in A. glycines to Virus Exposure

Most previous studies have found that a major component of insect vector transcriptional responses to plant virus exposure is immune related (Luan et al. 2011; Whitfield et al. 2011; Xu et al. 2012; Cassone et al. 2014). This prompted us to manually inspect the candidate lists for transcripts assigned as immune related in the de novo assembly (Waterhouse et al. 2007; Bartholomay et al. 2010; Waterhouse et al. 2010) whether or not they were assigned to clusters of enriched annotation terms by DAVID. Of the 653 immune-linked genes in the D. melanogaster (BDGP5) assembly, only 114 were detected in the A. glycines transcriptome assembly (tBLASTx, E value < 10−10). Comparing between treatments at 4 h, 34 were differentially expressed, of which 24 (71%) were downregulated in aphids fed on BPMV-infected plants (table 3). However, the overrepresentation of downregulated immune genes in the BPMV treatment at 4 h was not significant (χ2 = 0.225 [1], P = 0.636). Even fewer immune transcripts (n = 22) were differentially expressed at 7 days, and ten of these were up/downregulated at both time points in the BPMV and/or SMV treatments. The persistently downregulated transcripts mostly belonged to peroxidase and autophagy gene families (data not shown).

Table 3

Immune-Related Transcripts Up/Downregulated in Aphis glycines After Feeding on SMV-Infected or BPMV-Infected Soybean for 4 h

Immune Gene Familya	BPMV/Upb	BPMV/Downb	SMV/Upb	SMV/Downb
Antimicrobial peptides	–	–	–	–
Autophagy genes	Atg6, Atg18c	Atg5c, Atg6, Atg7, Atg9, Atg12c, APHAG1	–	Atg5c
Caspases	–	Nc	–	Icec
Clip domain serine proteases	CLIP37	–	–	–
IMD pathway	–	Casp	–	–
JAKSTAT signal transduction	–	Hop	–	–
Inhibitors of apoptosis	–	dIAP1	–	dIAP1
Relish-Like proteins	–	dl-RC	–	–
Peroxidases	Prx2540, Pdx	Jafrac2c, PHGPx	–	PRDX16c, Pdxc, PRDX19c
Scavenger receptors	emp	CRQ, SCR19	–	CRQ
Serine protease inhibitors	–	–	–	Spn3
Superoxide dismutases	Sod3c	–	–	–
SRRP	–	LOQS, Dcr1, SPNE, PASHA, Dcr2c, AUB, AGO3	–	–
Thioester-containing proteins	–	TEP3	–	–
Toll-like receptors	–	–	–	Cact
Toll pathway members	–	Myd88	–	–

Note.–Functional annotation was based on Drosophila prediction using tBLASTX (E value < 10−10).

aGene assignments were based on Waterhouse et al. (2007, 2010) and Bartholomay et al. (2010).

bDifferentially expressed transcripts identified by one-way ANOVA as up/downregulated in one treatment relative to the other virus treatment and healthy control.

cGenes that were also persistently up/dowregulated 7 days after feeding on SMV-infected or BPMV-infected soybean.

Although the A. glycines transcriptome appears to lack transcripts annotated as immune related, one exciting pattern emerged: there was a significant overrepresentation of downregulated transcripts in the small RNA regulatory pathways (SRRPs) in aphids fed on BPMV-infected soybean for 4 h (χ2 = 9.451 [1], P = 0.002). SRRPs include the microRNA, small interfering RNA (siRNA), and PIWI small RNA pathways (Campbell et al. 2008), and components of all three were represented in the downregulated group: microRNA group transcripts were Dcr1, PASHA, and LOQS; siRNA group transcripts Dcr2; and PIWI RNA transcripts SPNE, AGO3, and AUB. Additional siRNA (AGO2, R2D2) and PIWI (PIWI, ARM1, and Rm62F) transcripts were identified as downregulated using ortholog prediction from the Anopheles gambiae PEST (AgamP3.7) assembly.

RT-qPCR Validation of Differentially Expressed A. glycines Transcripts

Differential expression of six transcripts was verified by RT-qPCR using two laboratory colonies of A. glycines. The transcripts were randomly selected from the subset of transcripts that were identified by RNA-Seq as constitutively downregulated in A. glycines fed on BPMV infected soybean for 4 h and 7 days. In addition, one transcript not differentially expressed between treatments at either time point served as a control. Primer pairs for the selected transcripts were developed with primer pair PCR efficiencies between 1.88 and 2.1 (table 4). All six transcripts were confirmed to be significantly downregulated in the BPMV treatment relative to the SMV and control treatments, with similar magnitudes of expression difference between the different quantification methods.

Table 4

Comparison Gene Expression for Four Transcripts Significantly Downregulated in Aphis glycines Fed on BPMV-Infected Soybean Using RNA-Seq and RT-qPCR

Transcript IDa	Functionb	Orthologc	RNA-Seq 4 hd	RT-qPCR 4 he	RT-qPCR 4 hf	RNA-Seq 7 daysd	RT-qPCR 7 dayse	RT-qPCR 4 hf
895g	Translation	Acyrthosiphon pisum	−0.09	−0.02	0.11	0	−0.09	0.07
1551	RNA processing	Ac. pisum	−3.04	−2.92	−2.43	−2.08	−2.54	−2.38
2413	NA	NA	−2.64	−3.58	−3.01	−2.78	−3.49	−2.64
2441	NA	NA	−2.61	−2.78	−2.75	−1.66	−1.89	−1.99
2559	NA	NA	−2.97	−3.90	−3.04	−1.97	−2.40	−1.82
4576	lipase maturation factor 2-like	Ac. pisum	−2.61	−2.88	−2.91	−2.31	−2.39	−2.22
14910	NA	NA	−3.96	−3.41	−3.51	−3.43	−2.74	−1.98

Note.—NA, not available.

aTranscript IDs derived from de novo assembly (supplementary table S1, Supplementary Material online).

bFunctional description derived from orthologs identified by BLASTx against nr database

cTop ortholog match from nr database.

dFold-change expression in BPMV treatment relative to healthy control calculated from RPKM from RNA-Seq.

eFold-change expression in BPMV treatment relative to healthy control of biotype 3 aphids calculated from 2−ΔΔT.

fFold-change expression in BPMV treatment relative to healthy control of biotype 1 aphids calculated from 2−ΔΔT.

gNegative control not differentially expressed by RNA-Seq or RT-qPCR.

BPMV Does Not Replicate in A. glycines

BPMV is not thought to have strong interaction with nor replicate in either vector or nonvector insects. However, the substantial phenotypic, genetic, and molecular responses of A. glycines to BPMV-infected soybean were surprising. To rule out the possibility that BPMV was strongly associated with or replicating in the aphid, RT-PCR assays were carried out on aphids fed on BPMV-infected plants for 7 days to allow virus access to the insect, followed by a 7 h starvation and 1.5 h feeding period on healthy soybean to clear the insect gut. BPMV was not detected in any of the five replicates examined (supplementary fig. S3, Supplementary Material online). These data indicate that, as expected, there was no strong association of the virus with the aphid even after 7 days of exposure.

Discussion

Aphis glycines can produce up to 15 generations on a soybean plant, all of which are primarily composed of apterous females produced through viviparous parthenogenesis (Venette and Ragsdale 2004). Thus, apterous females can spend their entire life on a single soybean plant, suggesting chronic exposure to virus-infected plants is not uncommon under natural conditions. In these experiments, feeding for 7 days on BPMV- or SMV-infected soybean did not affect adult aphid survival. However, aphids feeding on BPMV-infected soybean produced considerably fewer viable offspring than aphids feeding on healthy or SMV-infected plants over the same period. Little is known about the effects of nonpersistently transmitted viruses on plant–aphid interactions. Some evidence indicates that aphid population growth (Blua and Perring 1992; Mauck et al. 2010) and survival (Hodge and Powell 2008) is reduced on plants infected by the potyviruses Zucchini yellow mosaic virus and Bean yellow mosaic virus. In contrast, aphids fed on potato infected with the nonpersistently transmitted Potato virus Y showed an increase in population growth (Castle and Berger 1993). Our results indicate that feeding on SMV-infected plants for 7 days did not impact aphid fecundity and survival.

Negative impacts to fitness were previously detected for A. glycines feeding on soybean infected with BPMV (Donaldson and Gratton 2007). Changes in insect vector fecundity have been found for other plant–virus–vector pathosystems, but no clear patterns have emerged. In some instances, insect vector fecundity is reduced by feeding on virus infected plants (Rubinstein and Czosnek 1997; Mann et al. 2008; Sidhu et al. 2009), yet in others, the opposite is observed (McKenzie 2002; Jiu et al. 2007). Recently, reductions in fitness of nonvector species feeding on virus-infected plants were observed, but fecundity was unchanged (Tu et al. 2013; Xu et al. 2013). Additional studies are required to infer whether the observed decreases in A. glycines fecundity resulted from BPMV-induced changes in host plant phloem components (e.g., nutritional substances, toxins, and miRNAs) (Tu et al. 2013) or were due to other interactions occurring inside the aphid vector.

To explore the genetic and molecular underpinnings associated with the observed differences in fecundity between A. glycines fed on BPMV-infected, SMV-infected, and healthy soybean, transcriptome profiles were generated for each treatment, and changes in accumulation were assayed after 4 h and 7 days of feeding. There was a striking downregulation of A. glycines transcripts within hours of feeding on BPMV-infected soybean, and this downregulation persisted at 7 days. The downregulated transcripts encode proteins that predominately function in regular cellular processes, such as intracellular transport, transcription, and helicase activity (fig. 4). In contrast to this downregulation was increased expression of transcripts required for cellular respiration and glycolysis, suggesting that aphids may be diverting resources from parthenogenetic reproduction and other cellular activities to energy production to survive on BPMV-infected soybean.

Aphis glycines fed on SMV-infected soybean had transcriptome profiles that were more similar to the healthy control treatment, which was consistent with our expectations for a nonpersistently transmitted virus. These potyviruses interact only briefly with the vector, usually in the stylet or foregut, and do not cross vector membranes or other insect barriers, so the genetic and molecular responses of the vector to these viruses were not expected to be extensive. One notable difference between aphids feeding on healthy and SMV-infected plants was the downregulation of transcripts associated with regular cellular activities at 7 days, which was similar to, but less extensive than, that seen in the BPMV treatment. This may reflect aphid responses to SMV-induced deterioration of host plant quality, as has been found for other nonpersistently transmitted viruses (Hodge and Powell 2008; Mauck et al. 2010).

From an evolutionary perspective, the detrimental impacts of feeding on BPMV-infected soybean to A. glycines may reflect the short tritrophic histories among the insect, plant, and virus. Both the aphid and soybean are native to East Asia, whereas BPMV is largely confined to the Americas, indicating that their ecological, environmental, and genetic interactions span only a few years (Zaumeyer and Thomas 1948; Ghabrial et al. 1977; Mabry et al. 2003). SMV was established throughout Asia prior to the A. glycines North American invasion (Li et al. 2010) and thus has interacted with the plant and vector for considerably longer than has BPMV.

The planthopper and leafhopper responses to feeding on virus-infected hosts include upregulation of immune-related genes (Luan et al. 2011; Xu et al. 2012; Cassone et al. 2014). However, similar to the pea aphid, Ac. pisum (Gerardo et al. 2010), A. glycines appears to have a reduced and simplified humoral immune system. Notably absent from the transcriptome were transcripts belonging to several immune-related families, including antimicrobial peptides, caspase activators, and lysozymes. There were single putative members for several other families (e.g., peptidoglycan receptor proteins, fibrinogen-related proteins, galactoside-binding lectins, MD2-like receptors, and 1,3-beta-d glucan-binding proteins). Most aphids require symbiosis with the g-proteobacterium B. aphidicola to induce an immune response that is deficient in the aphid host (Gerardo et al. 2010; Douglas et al. 2011). Accordingly, most prokaryote transcripts in A. glycines had similarity to B. aphidicola, and this symbiont likely plays an integral role in its pathogen defense response.

Most insect-vectored plant viruses infect only the plants that the vector feeds upon. Only a handful of plant viruses are known to also infect aphids (Hogenhout et al. 2008), all of which are enveloped viruses in the family Rhabdoviridae. A recent study (Li et al. 2014) indicated that the pollen-borne plant virus, Tobacco ringspot virus (TRSV), replicates in honeybees and reduces fecundity in these insects. Interestingly, TRSV and BPMV belong to the same subfamily, Comoviridae, within the family Secoviridae. Two major trends in our results initially suggested that BPMV replicates in A. glycines. First, the substantial transcriptome downregulation was coupled with elevated expression levels of translation and cytoskeleton organization transcripts that persisted at 7 days. BPMV has a positive-sense single-stranded RNA genome (Giesler et al. 2002). As such the viral genomic RNA is directly translated after entry into a cell (Wang and Li 2012). At least some insects have the cellular machinery required for replication of positive strand RNA viruses (Scotti et al. 1983; Ball and Johnson 1998; Dasgupta et al. 2003; Aliota and Kramer 2012). Moreover in plants, virus-induced cytoskeleton organization is a common strategy for intra- and intercellular virus movement (Ruthel et al. 2005; Wei et al. 2006; Boutte et al. 2007; Harries et al. 2010; Schoelz et al. 2011). Second, most SSRP transcripts were downregulated at 4 h, including all components of the PIWI RNA pathway. The siRNA and miRNA pathways are important in aspects of the antiviral defense in animal-infecting arthropod vectors (Fire et al. 1998; Campbell et al. 2008; Schnettler et al. 2013). Less is known about the role of the PIWI RNA pathway in antiviral responses, but recent studies implicated the pathway in inhibiting replication of positive-strand RNA viruses in mosquitoes (Hess et al. 2011; Schnettler et al. 2013). However, following a period of feeding on healthy plants to clear the insect gut of virus, BPMV could not be detected in A. glycines fed on BPMV-infected plants for 7 days using RT-PCR.

Conclusions

To our knowledge, this is the first study to explore the transcriptional basis of changes in fitness of an insect vector after exposure to a transmitted virus (SMV) or a nontransmitted virus (BPMV). Aphis glycines fed on SMV-infected soybean had transcriptome and fitness profiles similar to the healthy control treatment, which was expected for a nonpersistently transmitted virus. However, the sizeable and persistent downregulation of the transcriptome and associated significant reduction in fecundity of aphids fed on BPMV-infected soybean were unexpected given that A. glycines does not transmit BPMV. Suppressing metabolically intensive activities such as reproduction may be essential for aphid survival on soybean negatively impacted by BPMV infection. BPMV does not apparently replicate in A. glycines, despite the transcriptional changes that suggest a defense response. The detrimental impacts to A. glycines fed on BPMV-infected soybean may be attributed to their short tritrophic evolutionary histories, which span only a few years. Future fitness bioassays, behavioral bioassays, and targeted metabolomic and proteomic approaches are needed to characterize any direct and indirect effects of BPMV exposure, as well as provide a better understanding of the genetic and molecular basis of these virus-associated differences. More broadly, our results highlight the need for studies aimed at characterizing the virus transmission dynamics associated with native and invasive species pathosystems.

Supplementary Material

Supplementary figures S1–S3 and tables S1–S3 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).