# Expression differences in Aphidius ervi (Hymenoptera: Braconidae) females reared on different aphid host species

Gabriel I Ballesteros, Jürgen Gadau, Fabrice Legeai, Angelica

Gonzalez-Gonzalez, Blas Lavandero, Jean-Christophe Simon, Christian C
Figueroa

<figure><img src="image_1.png" /></figure>

To cite this version:

Gabriel I Ballesteros, Jürgen Gadau, Fabrice Legeai, Angelica Gonzalez-Gonzalez, Blas Lavandero, et al.. Expression differences in Aphidius ervi (Hymenoptera: Braconidae) females reared on different aphid host species. PeerJ, 2017, 5, 10.7717/peerj.3640. [https://doi.org/10.7717/peerj.3640](https://doi.org/10.7717/peerj.3640). [arXiv:2005.04601](https://arxiv.org/abs/2005.04601)

HAL Id: hal-02628911

https://hal.inrae.fr/hal-02628911

Submitted on 27 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Expression differences in *Aphidius ervi* (Hymenoptera: Braconidae) females reared on different aphid host species

Gabriel I. Ballesteros<sup>1,2,*</sup>, Jürgen Gadau<sup>3,4,*</sup>, Fabrice Legeai<sup>5,6</sup>, Angelica Gonzalez-Gonzalez<sup>1,2</sup>, Blas Lavandero<sup>1</sup>, Jean-Christophe Simon<sup>6</sup> and Christian C. Figueroa<sup>1,2</sup>

Submitted 7 April 2017
Accepted 12 July 2017
Published 21 August 2017

Corresponding author
Christian C. Figueroa,
alfigueroa@utalca.cl

Academic editor
Mikhail Gelfand

Additional Information and
Declarations can be found on
page 20

DOI 10.7717/peerj.3640

cc Copyright
2017 Ballesteros et al.

Distributed under
Creative Commons CC-BY 4.0

OPEN ACCESS

## ABSTRACT

The molecular mechanisms that allow generalist parasitoids to exploit many, often very distinct hosts are practically unknown. The wasp *Aphidius ervi*, a generalist koinobiont parasitoid of aphids, was introduced from Europe into Chile in the late 1970s to control agriculturally important aphid species. A recent study showed significant differences in host preference and host acceptance (infectivity) depending on the host *A. ervi* were reared on. In contrast, no genetic differentiation between *A. ervi* populations parasitizing different aphid species and aphids of the same species reared on different host plants was found in Chile. Additionally, the same study did not find any fitness effects in *A. ervi* if offspring were reared on a different host as their mothers. Here, we determined the effect of aphid host species (*Sitobion avenae* versus *Acyrthosiphon pisum* reared on two different host plants alfalfa and pea) on the transcriptome of adult *A. ervi* females. We found a large number of differentially expressed genes (between host species: head: 2,765; body: 1,216; within the same aphid host species reared on different host plants: alfalfa versus pea: head 593; body 222). As expected, the transcriptomes from parasitoids reared on the same host species (pea aphid) but originating from different host plants (pea versus alfalfa) were more similar to each other than the transcriptomes of parasitoids reared on a different aphid host and host plant (head: 648 and 1,524 transcripts; body: 566 and 428 transcripts). We found several differentially expressed odorant binding proteins and olfactory receptor proteins in particular, when we compared parasitoids from different host species. Additionally, we found differentially expressed genes involved in neuronal growth and development as well as signaling pathways. These results point towards a significant rewiring of the transcriptome of *A. ervi* depending on aphid-plant complex where parasitoids develop, even if different biotypes of a certain aphid host species (*A. pisum*) are reared on the same host plant. This difference seems to persist even after the different wasp populations were reared on the same aphid host in the laboratory for more than 50 generations. This indicates that either the imprinting process is very persistent or there

<footnote n="1"><sup>1</sup> Instituto de Ciencias Biológicas, Universidad de Talca, Talca, Chile</footnote>

<footnote n="2"><sup>2</sup> Millennium Nucleus Centre in Molecular Ecology and Evolutionary Applications in the Agroecosystems, Universidad de Talca, Talca, Chile</footnote>

<footnote n="3"><sup>3</sup> School of Life Sciences, Arizona State University, Tempe, AZ, United States of America</footnote>

<footnote n="4"><sup>4</sup> Institute for Evolution and Biodiversity, Westfälische Wilhelms-Universität Münster, Münster, Germany</footnote>

<footnote n="5"><sup>5</sup> GenScale, INRIA Centre Rennes, Rennes, France</footnote>

<footnote n="6"><sup>6</sup> Institute of Genetics, Environment and Plant Protection, INRA, Le Rheu, France</footnote>

<footnote>* These authors contributed equally to this work.</footnote>
is enough genetic/allelic variation between *A. ervi* populations. The role of distinct molecular mechanisms is discussed in terms of the formation of host fidelity.

**Subjects** Animal Behavior, Entomology, Genomics

**Keywords** Phenotypic plasticity, Local adaptation, Parasitoid wasps, Transcriptome, Aphid pest control

## INTRODUCTION

Parasitoids are widely used in biological control programs, which are based on introduced, naturalized, natural or released parasitoid wasps (Starý et al., 1993). Most parasitoids used for controlling aphids are endoparasitoids, which lay eggs inside their host as part of their life cycle and eventually killing it. Endoparasitoids may attack many related host species (e.g., same family) (Loxdale & Harvey, 2016). However, they have at the same time the problem to overcome or avoid different defense mechanisms of their varied host species in order to survive through adulthood (Jones et al., 2015). One way to do this is the formation of host races, i.e., populations that are adapted to a specific host (Stireman et al., 2006). Alternatively, the parasite may change its “phenotype” depending on the host it encounters and parasitizes, i.e., it shows adaptive phenotypic plasticity (Crispo, 2008). A third alternative is that several host species can be parasitized by the parasitoid using the same strategy/phenotype. The former two mechanisms will lead to differences in host preference and host acceptance, affecting behavioral traits associated to host selection (Zepeda-Paulo et al., 2013) but should not necessarily affect or change parasitoid fitness (Rivero, 2000; Wang et al., 2016).

One of the most widely used species in biological control programs is *Aphidius ervi*, a worldwide distributed koinobiont endoparasitoid of several Macrosiphinae aphid species such as the pea aphid *Ayllorhosiphon pisum* (Henry et al., 2010; Stilmant et al., 2008) and the grain aphid *Sitobion avenae* (Cameron, Powell & Loxdale, 1984). The parasitoid *A. ervi* was introduced into Chile from Europe in the late 70's as part of an aphid biological control program in cereals. The introduced *A. ervi* wasps successfully parasitized both *A. pisum* on legumes (e.g., alfalfa) and *S. avenae* on cereals (e.g., wheat) (Starý, 1993; Zepeda-Paulo et al., 2013) although these two aphids differ in several important aspects (e.g., host range usage, body size and colour, semiochemicals present in the cuticle, cornicular secretions, defensive reactions, etc.) (Daza-Bustamante, Fuentes-Contreras & Niemeyer, 2003). The recent introduction of a small number of individuals and natural occurrence of *A. ervi* attacking different aphid species opens interesting questions regarding the molecular basis and evolution of host preference and host specific adaptations (host races). Adaptations to different hosts have been described and confirmed as a plausible speciation mechanism (ecological speciation; Abrahamson & Blair, 2008; Schlüter, 2000). Interestingly, whereas differences in host preference and acceptance have been described previously from parasitoids reared or collected from different aphid host species, no fitness effects were detected if those parasitoids were forced to lay their eggs in suitable hosts they were not reared on (Zepeda-Paulo et al., 2013). These results suggest that the preference
for the natal host (i.e., host fidelity) is not under direct selection and those parasitoids
may show adaptive phenotypic plasticity (Daza-Bustamante et al., 2002; Zepeda-Paulo et
al., 2013). Indeed, no host race specific differentiation has been detected in Chilean A. ervi
populations (Zepeda-Paulo et al., 2013) and a high gene flow predominantly mediated by
male dispersion was found between populations (Zepeda-Paulo et al., 2015). Therefore,
phenotypic plasticity should be playing a key role in the observed host fidelity (Henry,
Roitberg & Gillespie, 2008; Zepeda-Paulo et al., 2013). However, it is unclear how much
phenotypic plasticity the parasitoid needs to exploit different hosts or whether they can
use identical mechanisms and strategies when parasitizing different hosts. One way to test
between these options is to compare the transcriptomes of A. ervi females reared on different
hosts. If A. ervi uses the same strategy for parasitizing both hosts their transcriptome should
be very similar.

Despite its widespread use in applied and fundamental research, no genomic or
transcriptomic information are available for *A. ervi* (Colinet et al., 2014). This study uses
RNAseq to *de novo* assemble and annotate the first representative transcriptome for
*A. ervi*. Additionally, our experimental setup allowed us to identify differentially expressed
genes between *A. ervi* parasitizing two different host species (*Acyrthosiphon pisum* and
*Sitobion avenae*) for two different body parts (head and body) and the same host species
(*Acyrthosiphon pisum*) reared on two different host plants (alfalfa and pea). This will allow
us to determine how much phenotypic plasticity at the transcriptome level *A. ervi* shows
and whether it can use the same strategy when it parasitizes different hosts.

## MATERIALS AND METHODS

### Insect collection and rearing

Parasitized aphids (two host races of *Acyrthosiphon pisum* and *Sitobion avenae*) were
collected from fields of legumes and cereals in two different geographic zones in Chile:
Region de Los Rios (S39°51', W73°7') and Region del Maule (S35°24', W71°40'). The
alfalfa race of *A. pisum* (APA) was sampled on alfalfa (*Medicago sativa* L.) and the pea race
(APP) was sampled on pea (*Pisum sativum* L.) (*Peccoud et al., 2008*), while the grain aphid
*Sitobion avenae* (SA) was sampled on wheat (*Triticum aestivum* L.). Pea aphids (from both
alfalfa and pea races) were maintained in the laboratory on broad bean (*Vicia faba* L.) while
grain aphids were maintained on barley (*Hordeum vulgare* L.). Both host plants have been
used previously for aphid and parasitoid rearing in other studies (*Sepúlveda et al., 2017b;
Zepeda-Paulo et al., 2013*). Parasitized aphids were reared under laboratory conditions
that allowed continuous reproduction (20 °C, D16/N8 photoperiod) (*Zepeda-Paulo et
al., 2013*). *Aphidius ervi* parasitoids were collected as larvae from parasitized aphids,
recognizable as mummies and kept separated in vials until adult parasitoids emerged.
Species and sex of each emerging parasitoid was determined using a standard taxonomic
key (*Starý, 1995*). In order to establish inbred populations, a single, isolated naive *A. ervi*
virgin female was mated with a naive virgin male for 24 h in a petri dish with diluted
honey and water for sustenance. Mated females were then offered new aphid hosts *ad
libitum* in a separate cage. These inbred lineages were propagated for approximately 60
generations before samples were taken for the RNAseq experiments. All *A. ervi* parasitoids were reared on the same host aphid species from which they were originally collected (further on called natal host). Thus, three different and highly inbred *A. ervi* laboratory populations were established: (i) *A. ervi* population originally collected from *A. pisum* living on alfalfa (Ae-APA) (ii) *A. ervi* population originally collected from *A. pisum* living on pea (Ae-APP) and (iii) *A. ervi* population originally collected from *S. avenae* living on wheat (Ae-SA) (Fig. 1). Each week, new aphid infested plants were introduced into the *A. ervi* rearing cages for parasitoid population maintenance, together with vials containing diluted honey and water for adult parasitoid feeding. All aphid populations were free of known secondary endosymbionts; their presence was evaluated using the amplification of specific 16S rDNA from whole-body aphid DNA based on the set of known primers described by Peccoud et al. (2014). This method allows screening of different symbionts, including the protective bacteria *Hamiltonella defensa* and *Regiella insecticola*, and also *Serratia symbiotica*, *Rickettsia*, *Rickettsiella* and *Spiroplasma* (Sepúlveda et al., 2017a).

## RNA collection and sequencing

Adult female parasitoids were collected alive from each one of the three caged parasitoid populations (N = 20 per population) (Fig. 1) and stored in 1.5 ml centrifuge tubes containing RNALater (QIAGEN, Hilden, Germany) at -20 °C until dissection and RNA extraction. Heads and bodies were dissected on ice using a sterile scalpel and pooled in six different samples. For each sample, total RNA was extracted using the RNEasy Plant Mini Kit (QIAGEN, Hilden, Germany) following the manufacturer's instructions and eluted in 50 µl of RNAse free water. Total RNA was quantified by spectrophotometry (Epoch Microplate Spectrophotometer, Biotek) and fluorometry (Qubit 3.0, Qubit RNA Broad Range Assay Kit), and integrity checked in a Bioanalyzer 2100 RNA Nano Kit (Agilent, Santa Clara, CA, USA). Recovered total RNA was precipitated using 0.1 volumes of Sodium Acetate 3 M and 2 volumes of 100% Ethanol and shipped to Macrogen Korea for library preparation and sequencing. Ribosomal RNA was depleted from total RNA using the Ribo-Zero rRNA Removal Kit for enrichment of both insect mRNA and non poly-adenilated mRNA that might be present in *A. ervi* sequenced samples. Remaining RNA was used for library construction using TruSeq Stranded Total RNA Sample Preparation Kit (Illumina, San Diego, CA, USA), tagged, pooled and sequenced using an Illumina HiSeq 2000 (2 × 100 bp, Paired End libraries; Macrogen, Seoul, South Korea). Raw transcriptome data was deposited in NCBI's Sequence Read Archive database under BioProject ID: PRJNA377544.

## Illumina sequence processing and transcriptome assembly

Illumina RNA-seq libraries were quality checked with FastQC ver. 0.11.3 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) in order to assess the presence of adapters derived from sequencing, overrepresented kmers, read length and overall read quality scores. All libraries were processed with Trimmomatic ver. 0.35 (Bolger, Lohse & Usadel, 2014) to remove any remaining TruSeq adapter sequence and to eliminate low quality bases (Q < 3) from reads. After sequence processing, all remaining sequences
<figure><img src="image_3.png" /></figure>

**Figure 1** Sampling design for RNA sequencing and differential expression analysis in the aphid parasitoid wasp *Aphidius ervi*. (A) Parasitoid population rearing. (B) DE analysis between heads and bodies (three libraries/tissue). (C) DE analysis between parasitoid host races.
shorter than 36 bp long were also removed from all datasets. Clean Illumina datasets were pooled *in silico* by concatenating library files. Before assembly, ribosomal RNA reads were removed by mapping the libraries using Bowtie ver. 1.1.1 (Langmead et al., 2009) against a custom rRNA database created from insect ribosomal sequences downloaded from NCBI and keeping non-mapped reads. The remaining high quality reads were *de novo* assembled with Trinity ver. 2.0.6 using default parameters. Metrics for *de novo* assembly were obtained with QUAST ver. 2.3 (number of contigs, total length, N50, largest contig, %GC, etc.) (Gurevich et al., 2013) while transcriptome completeness was assessed by benchmarking the assembled transcriptome using BUSCO (benchmarking universal single-copy orthologs) v.1.1b1 (Simão et al., 2015). To determine whether this transcriptome encodes for one or more set of core genes conserved across a range of Arthropod species, a “completeness score” was calculated (Moreton, Izquierdo & Emes, 2016). A total of 2,675 near-universal single-copy orthologs from Arthropod species were used as reference core genes (available at busco.ezlab.org; Simão et al., 2015). Additionally, reference protein sequences of *A. ervi* were downloaded from both the Non Redundant (NR) and Transcriptome Shotgun Assembly (TSA) sequence databases available at NCBI (N = 422; Colinet et al., 2014). These sequences were used as a custom reference protein database for BLASTx alignments using the assembled transcriptome as query.

## Annotation and functional gene classification

Homology searches of contigs from the assembled *de novo* transcriptome were performed locally with BLASTX using the NR database (NCBI) as reference (April 2016 version), setting an *e*-value of 1e-5 as threshold. Any contig that showed homology with at least one gene or protein was designated as a hit contig, while contigs with no hits were disregarded for any of the follow up analysis. Contigs with top-hit to non-insect species (e.g., prokaryote, yeast, vertebrates, etc.) were removed from the assembly and stored separately for future analysis because we were focused only on insect genes. Insect species reference lists available in NCBI Taxonomy browser were used to identify contigs with top-hits for insects (search criteria “insecta”, 200 levels displayed). To further improve the accuracy of differentially expressed genes, all contigs aligning to the same protein were grouped using BLASTx homology results and were sorted by alignment BitScore. The sequence with the highest BitScore was considered as the best blast hit, selected and designed as an annotated contig (Ono et al., 2015). The insect-filtered, non-redundant contig fasta dataset was loaded into Blast2GO ver. 2.8 (Conesa & Götz, 2008) altogether with BLASTx results in XML format. We also performed InterPro annotation, Gene Ontology (GO) term assignment, enzyme code and pathway annotation using Kyoto Encyclopedia of Genes and Genomes (KEGG) term integrated into Blast2GO. Successfully annotated transcripts were categorized and assigned to GO terms from different GO categories (molecular function, cellular component and biological process). The final contig annotation table was obtained from the combination of “Top-Blast table” and “sequence table”, both exported from Blast2GO.
## Differential gene expression analysis

The insect filtered, non-redundant reference transcriptome was used as a basis for differential gene expression studies between tissues (separate for head and body) (Fig. 1). Note: each individual sample was based on 20 individual females. Additionally, differential expression studies were carried out in order to detect and describe unique expressed transcripts for parasitoid lines from different aphid hosts (Ae-APP = *A. pisum*—Pea; Ae-APA = *A. pisum*—Alfalfa; Ae-SA = *S. avenae*) (Fig. 1). Gene abundance estimation was performed by separately aligning the libraries to the reduced reference transcriptome using the align_and_estimate script included in Trinity (ver. 2.0.6). This script automated the reference transcript, performed library read alignment to the reference using Bowtie2 (ver. 2.2.4; Langmead et al., 2009), and estimated read abundances from mapping results per library with the RSEM package included in Trinity (ver. 2.0.6.; Li & Dewey, 2011, and it was used to interpret and analyze Bowtie2 mapping results. RSEM was used to combine each count matrix and to build up a raw transcript expression matrix and a TMM-normalized expression matrix (script abundance_estimates_to_matrix.pl); this raw counts matrix was further used for Differential expression (DE) analysis at tissue level (heads vs. bodies, three libraries per sample), aphid host species and host rearing plant species. DE analysis was performed with edgeR Bioconductor package implemented in R using the provided run_DE_analysis script in Trinity ver. 2.0.6. The package edgeR was selected as it has a relatively high sensitivity and specificity in DE analysis of pooled samples compared to other methods of analysis. Genes that had at least 4-fold change values with a FDR-corrected *p*-value of 0.01 or lower were considered as significantly differentially expressed between libraries/tissues. The annotation of DE contigs was performed by combining TMM-normalized expression profiles for each contig with the annotated transcriptome tables generated in Blast2GO together with BLASTx results. Uniquely, differentially expressed contigs were detected for all libraries. In the case of the SA lineage (both for heads and bodies), a contig was considered as overexpressed only if its expression profile was 4-fold times higher when compared to both APA and APP populations. For both *A. pisum* races, contigs were considered as overexpressed only if their fold-change was at least 4-fold times higher (FDR-corrected *p*-value < 0.01) compared to the SA race in both APA and APP populations. The GO term enrichment was performed using Fisher's exact test in Blast2GO.

## RESULTS

### *Aphidius ervi* reference transcriptome assembly

We generated a *de novo* transcriptome for *A. ervi* using transcriptomic datasets obtained from the sequencing of six Illumina libraries (NCBI SRA accession PRJNA377544). These datasets were obtained from pooled female adult parasitoids reared from different aphid hosts (APP = *A. pisum*—Pea; APA = *A. pisum*—Alfalfa; SA = *S. avenae*) (Table 1). All transcripts were concatenated to generate a reference transcriptome library that could act as further reference for other studies involving *A. ervi* (GenBank accession GFLW01000000). The RNAseq generated 431,237,778 raw reads (2 × 100 bp, paired end libraries). After
<table><caption>**Table 1** Summary of *Aphidius ervi* transcriptomic libraries and assembly statistics.</caption><thead><tr><th>Sequencing<br>Library sequenced</th><th>Raw reads</th><th>Filtered Reads</th></tr></thead><tbody><tr><td>Ae-APA</td><td>122,819,778</td><td>115,025,660</td></tr><tr><td>Ae-APP</td><td>124,329,988</td><td>113,456,324</td></tr><tr><td>Ae-SA</td><td>184,088,012</td><td>161,916,086</td></tr><tr><td>Total</td><td>431,237,778</td><td>390,398,070</td></tr><tr><td>Minus rRNA</td><td></td><td>237,214,294</td></tr></tbody></table>

<table><thead><tr><th>Assembly</th><th>Reference de novo</th><th>Blast Hits</th><th>Insect only</th></tr></thead><tbody><tr><td>Total transcript number</td><td>135,676</td><td>33,853</td><td>17,763</td></tr><tr><td>Total length</td><td>91,710,298</td><td>47,294,638</td><td>25,440,092</td></tr><tr><td>longest transcript</td><td>19,479</td><td>19,479</td><td>19,479</td></tr><tr><td>N50</td><td>1,516</td><td>2,027</td><td>1,983</td></tr><tr><td>%GC</td><td>30.15%</td><td>31.02%</td><td>30.67%</td></tr></tbody></table>

pre-processing and removal of ribosomal RNA reads, 237,214,294 reads remained. Using Trinity, filtered reads were assembled into 135,659 contigs (N50 of 1,516 bp, mean length of 675.95 bp). Among the assembled contigs, 38,567 were less than 300 bp (28.4%), 75,458 contigs were between 301 bp and 1,000 bp (55.6%), while 21,634 contigs had a size over 1,000 bp (16%). The assembly completeness was also evaluated using BUSCO (benchmarking universal single-copy orthologs) showing that 76.6% complete conserved genes were found in our assembly, 13.7% corresponded to fragmented conserved genes while only 9.8% of single-copy ortholog genes were missing. Our results are similar to BUSCO metrics reported for other *de novo* insect transcriptomes assemblies such as the Western tarnished plant bug *Lygus hesperus* assembly performed with Illumina datasets (Tassone et al., 2016). Additionally, all previously known *A. ervi* sequences (N = 422) retrieved from NCBI were present in our transcriptomic assembly.

## Sequence annotation

The BLASTx alignments revealed that 33,853 contigs were annotated to a known protein within the NR database (24.9% of total contigs) (Table 1). Most transcript sequences with protein hits matched to other braconid endoparasitoids such as *Diachasma alloeum* (parasitoid of the apple maggot *Rhagoletis pomonella*), *Fopius arisanus* (parasitoid of Tephritid fruit flies) and *Microplitis demolitor* (parasitoid of noctuid larvae) (Fig. 2), all species for which *de novo* transcriptomes have been published (Burke & Strand, 2012; Calla et al., 2015). This dataset was further filtered to remove redundant contigs by using best-blast hit criteria (see ‘Materials and Methods’). Additionally, annotated sequences from non-insect organisms were removed; both steps filtered out 16,090 contigs. The remaining 17,763 contigs were the basis for all follow up analyses (annotation and differential expression). A total of 10,492 contigs (59.1%) could be annotated based on their sequence homology with GO terms. As contigs can be assigned to more than one GO category, 14,614 contigs were assigned to biological process, 7,945 contigs were classified under
<figure><img src="image_4.png" /><figcaption>**Figure 2** Species distribution of unigene sequences of the aphid parasitoid wasp *Aphidius ervi* transcripts to other insect species using homologous BLASTx hits and NR-NCBI database.</figcaption></figure>

molecular function, and 5,976 were classified in cellular component (Fig. 3). This GO term distribution is congruent with other insect transcriptomes already sequenced (Hwang et al., 2016).

**Transcriptomic differences between tissues and function of transcripts with different expression levels**

We used our filtered transcriptome as a reference dataset to perform Differential Expression (DE) analysis for *A. ervi* reared on the three different hosts. Our results indicate that there is a good correlation within head samples and within body samples, while clear differential expression patterns are observed for different body parts (Fig. 4; Fig. S1). Differential expression analysis showed that 6,389 transcripts are being differentially expressed between head and body samples (N = 3,445 up-regulated in heads, N = 2,944 up-regulated in bodies). The GO Term enrichment analysis was performed between differentially expressed genes and indicated the most enriched GO terms for bodies and heads (Fig. 4). “Signal transduction” and “signal transducer activity” were among the most enriched GO terms for heads, hence indicating that signaling and response stimulus associated transcripts are prevalent in the head transcriptomes, which is expected as RNA was extracted from whole heads including the brain and chemical sense organs such as antennae (Glaser et al., 2015). In the case of bodies, “ribosome biogenesis” and “peptidase activity”, were, amongst
<figure><img src="image_5.png" /><figcaption>**Figure 3** Gene Ontology (GO) annotations for the reference transcriptome of *Aphidius ervi* separated by GO categories (Biological process: blue. Molecular function: red. Cellular component: green).</figcaption></figure>

others, the most enriched GO terms found. As we extracted RNA from headless bodies without further dissection, the enriched “peptidase activity” could probably reference to venom proteins such as serine proteases or gamma-glutamyl transpeptidases, which have also been described previously for *A. ervi* (*Colinet et al., 2014*). These venom proteins are injected into the host at oviposition and would have a role in the modulation of the aphid physiology by inducing apoptosis of host ovaries and arresting host reproduction (*Colinet et al., 2014; Falabella et al., 2007*).

**Transcriptomic differences between parasitoid lines and function of transcripts with different expression levels**

At body level, we found 239 transcripts with differential expression patterns between *A. ervi* populations that were originally collected from *A. pisum* on alfalfa and pea, respectively (149 up-regulated genes for Ae-AP and 90 up-regulated genes in Ae-SA; *q*-value < 0.01, fold-change > 4), while at head level 390 transcripts showed differential expression (219 upregulated in Ae-AP and 171 upregulated in Ae-SA). All differentially expressed transcripts were annotated using the results from BlastX alignments while GO terms were assigned using Blast2GO mapping results, which in turn was used to identify the functions of genes displaying different expression patterns between populations and tissues (Figs. 5 and 6); this approach was used because GO term enrichment analysis failed to find any enriched
<figure><img src="image_6.png" /></figure>

**Figure 4** Differential GO term distribution between *Aphidius ervi* bodies and heads (Blast2GO Fisher's exact test with FDR correction). Reference set: full *A. ervi* transcriptome. Test set: % of sequences associated to GO-enriched terms. (A) enriched GO-terms in bodies. (B) enriched GO-terms in heads. CC, Cellular Component; MF, Molecular Function; BP, Biological Process.

term. Top 20 differentially expressed gene lists (ranked first according to fold-change) are reported in Table 2 (Complete lists in File S1).

## Identification of putative chemosensory and olfaction-related genes

Olfaction plays a crucial role in insect behavior such as mate recognition, foraging, host location and host discrimination or finding shelter in complex environments (Suh, Bohbot & Zwiebel, 2014). Behavioral differences in host preference and host acceptance based on aphid host species and plants have been reported previously for *A. ervi* (Zepeda-Paulo et al., 2013). This ability to differentiate between hosts almost certainly involves chemical signal perception (Takemoto, Kainoh & Takabayashi, 2011). As olfactory behavior differences depend on both odorant recognition and signal propagation/processing, genes coding for odorant perception such as Odorant Binding Proteins (OBPs), Chemosensory Proteins (CSPs) and Odorant Receptors (ORs) and genes coding for signal propagation such as voltage-gated sodium channels are prime candidate genes underlying the observed differences in host-preference in *A. ervi* (Zepeda-Paulo et al., 2013). This variation in the ability to perceive and respond to chemosensory cues would provide a target for
<table><caption>**Table 2** Top 20 genes with differential expression patterns between *Aphidius ervi*-AP and *A. ervi*-SA on two different tissues (heads and bodies).</caption><thead><tr><th>ID</th><th>Library</th><th>Sequence description</th><th>Log2-fold change</th><th>P-Value</th><th>FDR-adjusted p-value</th></tr></thead><tbody><tr><td>TR52610-c0_g1_i1</td><td>Ae-AP Body</td><td>sodium hydrogen exchanger 7 isoform x4</td><td>13.57</td><td>6.48E-20</td><td>3.73E-16</td></tr><tr><td>TR27559-c2_g1_i1</td><td>Ae-AP Body</td><td>rho gtpase-activating protein 190 isoform x1</td><td>13.27</td><td>1.18E-19</td><td>4.82E-16</td></tr><tr><td>TR36885-c0_g2_i1</td><td>Ae-AP Body</td><td>oxidoreductase glyr1 homolog</td><td>12.17</td><td>2.86E-16</td><td>3.71E-13</td></tr><tr><td>TR42270-c3_g5_i10</td><td>Ae-AP Body</td><td>calcium-activated potassium channel slowpoke</td><td>12.12</td><td>3.43E-16</td><td>6.80E-13</td></tr><tr><td>TR27536-c4_g1_i1</td><td>Ae-AP Body</td><td>disco-interacting protein 2 isoform x1</td><td>11.84</td><td>3.76E-10</td><td>7.65E-08</td></tr><tr><td>TR42293-c6_g1_i8</td><td>Ae-AP Body</td><td>lon protease mitochondrial isoform x1</td><td>11.44</td><td>3.38E-14</td><td>2.23E-11</td></tr><tr><td>TR52097-c0_g1_i3</td><td>Ae-AP Body</td><td>synaptojanin-1 isoform x1</td><td>11.14</td><td>2.84E-13</td><td>1.85E-10</td></tr><tr><td>TR30823-c0_g1_i18</td><td>Ae-AP Body</td><td>zinc finger protein rotund isoform x3</td><td>10.98</td><td>3.93E-12</td><td>9.37E-10</td></tr><tr><td>TR55124-c1_g1_i1</td><td>Ae-AP Body</td><td>ryanodine receptor isoform x6</td><td>10.95</td><td>2.18E-11</td><td>6.15E-09</td></tr><tr><td>TR55070-c8_g1_i6</td><td>Ae-AP Body</td><td>dynamin isoform x2</td><td>10.91</td><td>3.11E-12</td><td>7.91E-10</td></tr><tr><td>TR4006-c7_g2_i9</td><td>Ae-AP Body</td><td>down syndrome cell adhesion molecule-like protein dscam2 isoform x30</td><td>10.88</td><td>1.15E-12</td><td>4.56E-10</td></tr><tr><td>TR13038-c0_g3_i4</td><td>Ae-AP Body</td><td>excitatory amino acid transporter isoform x1</td><td>10.84</td><td>3.07E-12</td><td>8.35E-10</td></tr><tr><td>TR42270-c3_g5_i5</td><td>Ae-AP Body</td><td>calcium-activated potassium channel slowpoke isoform x7</td><td>10.75</td><td>8.45E-11</td><td>2.06E-08</td></tr><tr><td>TR37837-c6_g1_i2</td><td>Ae-AP Body</td><td>inorganic phosphate cotransporter isoform x1</td><td>10.74</td><td>6.56E-12</td><td>2.22E-09</td></tr><tr><td>TR28738-c0_g3_i2</td><td>Ae-AP Body</td><td>embryonic polarity protein dorsal-like isoform x3</td><td>10.66</td><td>4.77E-12</td><td>1.46E-09</td></tr><tr><td>TR20185-c0_g1_i1</td><td>Ae-AP Body</td><td>transposable element p transposase</td><td>10.64</td><td>1.09E-11</td><td>3.53E-09</td></tr><tr><td>TR29865-c1_g2_i1</td><td>Ae-AP Body</td><td>cordon-bleu 1</td><td>10.57</td><td>2.74E-10</td><td>5.85E-08</td></tr><tr><td>TR24635-c2_g2_i3</td><td>Ae-AP Body</td><td>zinc finger cchc domain-containing protein 4</td><td>10.48</td><td>6.09E-11</td><td>1.16E-08</td></tr><tr><td>TR16911-c9_g1_i6</td><td>Ae-AP Body</td><td>dynein heavy cytoplasmic isoform x2</td><td>10.26</td><td>2.05E-08</td><td>2.66E-06</td></tr><tr><td>TR19336-c1_g1_i6</td><td>Ae-AP Body</td><td>uncharacterized protein LOC100740589 isoform X3</td><td>10.18</td><td>3.73E-07</td><td>3.52E-05</td></tr><tr><td>TR41810-c12_g5_i10</td><td>Ae-SA Body</td><td>vinculin isoform x9</td><td>11.09</td><td>3.07E-12</td><td>1.12E-09</td></tr><tr><td>TR42270-c3_g5_i12</td><td>Ae-SA Body</td><td>calcium-activated potassium channel slowpoke isoform x6</td><td>10.99</td><td>5.80E-12</td><td>1.95E-09</td></tr><tr><td>TR21226-c0_g1_i1</td><td>Ae-SA Body</td><td>nfu1 iron-sulfur cluster scaffold mitochondrial-like</td><td>10.94</td><td>8.07E-12</td><td>2.58E-09</td></tr><tr><td>TR19336-c1_g1_i10</td><td>Ae-SA Body</td><td>fh1 fh2 domain-containing protein 3 isoform x4</td><td>10.83</td><td>1.69E-11</td><td>4.93E-09</td></tr><tr><td>TR31615-c2_g1_i1</td><td>Ae-SA Body</td><td>PREDICTED: uncharacterized protein LOC105456969</td><td>10.75</td><td>2.70E-11</td><td>7.44E-09</td></tr><tr><td>TR31615-c7_g1_i1</td><td>Ae-SA Body</td><td>ubiquitin carboxyl-terminal hydrolase 17-like partial</td><td>10.46</td><td>1.66E-10</td><td>3.67E-08</td></tr><tr><td>TR32717-c6_g1_i1</td><td>Ae-SA Body</td><td>paired amphipathic helix protein sin3a</td><td>9.93</td><td>3.40E-09</td><td>5.60E-07</td></tr><tr><td>TR4388-c16_g1_i1</td><td>Ae-SA Body</td><td>whirlin isoform x1</td><td>9.66</td><td>1.71E-08</td><td>2.32E-06</td></tr><tr><td>TR32865-c1_g1_i2</td><td>Ae-SA Body</td><td>a disintegrin and metalloproteinase with thrombospondin motifs 16 isoform x1</td><td>9.64</td><td>2.06E-08</td><td>2.71E-06</td></tr><tr><td>TR13691-c0_g1_i3</td><td>Ae-SA Body</td><td>membrane metallo-endopeptidase-like partial</td><td>9.61</td><td>7.62E-13</td><td>3.88E-10</td></tr><tr><td>TR45167-c1_g2_i2</td><td>Ae-SA Body</td><td>kielin chordin-like protein isoform x2</td><td>9.36</td><td>9.44E-08</td><td>1.07E-05</td></tr><tr><td>TR9096-c9_g1_i6</td><td>Ae-SA Body</td><td>piezo-type mechanosensitive ion channel component 1 isoform x2</td><td>9.28</td><td>6.66E-12</td><td>2.21E-09</td></tr></tbody></table>

(continued on next page)
<table><caption>Table 2 (continued)</caption><thead><tr><th>ID</th><th>Library</th><th>Sequence description</th><th>Log2-fold change</th><th>P-Value</th><th>FDR-adjusted p-value</th></tr></thead><tbody><tr><td>TR49009-c9_g1_i7</td><td>Ae-SA Body</td><td>voltage-dependent calcium channel subunit alpha-2 delta-3 isoform x2</td><td>9.26</td><td>1.66E-07</td><td>1.76E-05</td></tr><tr><td>TR47041-c0_g1_i4</td><td>Ae-SA Body</td><td>PREDICTED: uncharacterized protein LOC106789540 isoform X2</td><td>9.21</td><td>2.10E-07</td><td>2.14E-05</td></tr><tr><td>TR27559-c2_g1_i4</td><td>Ae-SA Body</td><td>rho gtpase-activating protein 190 isoform x1</td><td>9.02</td><td>6.48E-07</td><td>5.69E-05</td></tr><tr><td>TR23225-c0_g2_i6</td><td>Ae-SA Body</td><td>serine threonine-protein kinase ick-like isoform x2</td><td>8.85</td><td>1.49E-06</td><td>1.17E-04</td></tr><tr><td>TR37577-c0_g1_i2</td><td>Ae-SA Body</td><td>probable 28s rrna (cytosine-c )-methyltransferase</td><td>8.79</td><td>2.00E-06</td><td>1.52E-04</td></tr><tr><td>TR48980-c1_g2_i1</td><td>Ae-SA Body</td><td>creb-binding protein isoform x5</td><td>8.77</td><td>2.55E-10</td><td>4.86E-08</td></tr><tr><td>TR45475-c0_g2_i2</td><td>Ae-SA Body</td><td>ras-related protein m-ras-like</td><td>8.61</td><td>5.06E-06</td><td>3.38E-04</td></tr><tr><td>TR34000-c0_g1_i1</td><td>Ae-SA Body</td><td>neprilysin-2 isoform x1</td><td>8.56</td><td>9.29E-10</td><td>1.50E-07</td></tr><tr><td>TR3953-c4_g1_i15</td><td>Ae-AP Head</td><td>sorbin and sh3 domain-containing protein 1 isoform x3</td><td>14.77</td><td>5.90E-23</td><td>3.41E-19</td></tr><tr><td>TR27559-c2_g1_i1</td><td>Ae-AP Head</td><td>rho gtpase-activating protein 190 isoform x1</td><td>14.07</td><td>6.39E-22</td><td>1.19E-18</td></tr><tr><td>TR42270-c3_g5_i10</td><td>Ae-AP Head</td><td>calcium-activated potassium channel slowpoke</td><td>14.04</td><td>8.41E-21</td><td>1.52E-17</td></tr><tr><td>TR16911-c9_g1_i6</td><td>Ae-AP Head</td><td>dynein heavy cytoplasmic isoform x2</td><td>13.97</td><td>1.95E-20</td><td>2.81E-17</td></tr><tr><td>TR27544-c9_g3_i1</td><td>Ae-AP Head</td><td>bromodomain adjacent to zinc finger domain protein 2b-like isoform x7</td><td>13.31</td><td>3.90E-18</td><td>2.42E-15</td></tr><tr><td>TR41810-c12_g5_i6</td><td>Ae-AP Head</td><td>vinculin isoform x5</td><td>12.85</td><td>1.59E-18</td><td>1.10E-15</td></tr><tr><td>TR4006-c7_g2_i9</td><td>Ae-AP Head</td><td>down syndrome cell adhesion molecule-like protein dscam2 isoform x30</td><td>12.76</td><td>2.79E-18</td><td>1.77E-15</td></tr><tr><td>TR36885-c0_g2_i1</td><td>Ae-AP Head</td><td>oxidoreductase glyr1 homolog</td><td>12.74</td><td>1.20E-17</td><td>6.36E-15</td></tr><tr><td>TR52610-c0_g1_i1</td><td>Ae-AP Head</td><td>sodium hydrogen exchanger 7 isoform x4</td><td>12.7</td><td>2.18E-17</td><td>8.18E-15</td></tr><tr><td>TR9088-c5_g1_i4</td><td>Ae-AP Head</td><td>heterogeneous nuclear ribonucleoprotein k</td><td>12.16</td><td>1.66E-16</td><td>6.35E-14</td></tr><tr><td>TR52097-c0_g1_i3</td><td>Ae-AP Head</td><td>synaptojanin-1 isoform x1</td><td>12.07</td><td>1.21E-15</td><td>3.58E-13</td></tr><tr><td>TR53811-c0_g2_i6</td><td>Ae-AP Head</td><td>e3 ubiquitin-protein ligase hectd1 isoform x4</td><td>11.76</td><td>3.20E-15</td><td>8.91E-13</td></tr><tr><td>TR23225-c0_g2_i1</td><td>Ae-AP Head</td><td>serine threonine-protein kinase ick-like isoform x1</td><td>11.75</td><td>2.83E-15</td><td>7.29E-13</td></tr><tr><td>TR8499-c8_g1_i3</td><td>Ae-AP Head</td><td>netrin receptor unc5c</td><td>11.64</td><td>7.35E-12</td><td>7.20E-10</td></tr></tbody></table>

(continued on next page)
<table><caption>Table 2 (continued)</caption><thead><tr><th>ID</th><th>Library</th><th>Sequence description</th><th>Log2-fold change</th><th>P-Value</th><th>FDR-adjusted p-value</th></tr></thead><tbody><tr><td>TR13228-c0_g1_i1</td><td>Ae-AP Head</td><td>proteasome subunit alpha type-3</td><td>11.63</td><td>4.91E-14</td><td>1.03E-11</td></tr><tr><td>TR45093-c0_g1_i2</td><td>Ae-AP Head</td><td>adamts-like protein 4 isoform x2</td><td>11.47</td><td>3.37E-14</td><td>5.61E-12</td></tr><tr><td>TR11878-c8_g1_i2</td><td>Ae-AP Head</td><td>e3 ubiquitin-protein ligase nedd-4 isoform x1</td><td>11.03</td><td>1.05E-10</td><td>7.80E-09</td></tr><tr><td>TR45475-c0_g2_i1</td><td>Ae-AP Head</td><td>ras-related protein m-ras-like</td><td>11</td><td>3.97E-13</td><td>6.57E-11</td></tr><tr><td>TR16826-c0_g1_i1</td><td>Ae-AP Head</td><td>prenylcysteine oxidase-like</td><td>10.48</td><td>9.78E-11</td><td>1.26E-08</td></tr><tr><td>TR53868-c5_g2_i3</td><td>Ae-AP Head</td><td>caax prenyl protease 1 homolog</td><td>10.37</td><td>4.92E-08</td><td>3.92E-06</td></tr><tr><td>TR53728-c7_g1_i11</td><td>Ae-SA Head</td><td>sodium channel protein para isoform x10</td><td>13.95</td><td>4.98E-21</td><td>6.31E-18</td></tr><tr><td>TR49009-c9_g1_i5</td><td>Ae-SA Head</td><td>voltage-dependent calcium channel subunit alpha-2 delta-3 isoform x2</td><td>13.57</td><td>6.80E-20</td><td>6.57E-17</td></tr><tr><td>TR37911-c0_g2_i7</td><td>Ae-SA Head</td><td>focal adhesion kinase 1 isoform x1</td><td>12.6</td><td>4.95E-17</td><td>1.83E-14</td></tr><tr><td>TR43558-c0_g1_i1</td><td>Ae-SA Head</td><td>PREDICTED: uncharacterized protein LOC107045241</td><td>12.54</td><td>7.85E-17</td><td>2.77E-14</td></tr><tr><td>TR9014-c15_g1_i2</td><td>Ae-SA Head</td><td>a disintegrin and metalloproteinase with thrombospondin motifs 8 isoform x8</td><td>12.22</td><td>6.79E-16</td><td>1.88E-13</td></tr><tr><td>TR49009-c9_g1_i4</td><td>Ae-SA Head</td><td>voltage-dependent calcium channel subunit alpha-2 delta-3 isoform x1</td><td>12.01</td><td>2.64E-15</td><td>6.35E-13</td></tr><tr><td>TR31615-c2_g1_i1</td><td>Ae-SA Head</td><td>PREDICTED: uncharacterized protein LOC105456969</td><td>11.92</td><td>5.07E-15</td><td>1.12E-12</td></tr><tr><td>TR27559-c2_g1_i4</td><td>Ae-SA Head</td><td>rho gtpase-activating protein 190 isoform x1</td><td>11.36</td><td>3.46E-14</td><td>6.76E-12</td></tr><tr><td>TR23225-c0_g2_i6</td><td>Ae-SA Head</td><td>serine threonine-protein kinase ick-like isoform x2</td><td>11.24</td><td>7.57E-14</td><td>1.37E-11</td></tr><tr><td>TR42270-c3_g5_i15</td><td>Ae-SA Head</td><td>calcium-activated potassium channel slowpoke isoform x16</td><td>11.15</td><td>1.37E-13</td><td>2.39E-11</td></tr><tr><td>TR28916-c0_g3_i4</td><td>Ae-SA Head</td><td>liprin-beta-1 isoform x4</td><td>11.15</td><td>1.39E-13</td><td>2.42E-11</td></tr><tr><td>TR5686-c0_g1_i1</td><td>Ae-SA Head</td><td>RNA-directed DNA polymerase from mobile element jockey-like</td><td>11.09</td><td>2.08E-13</td><td>3.48E-11</td></tr><tr><td>TR4388-c16_g1_i1</td><td>Ae-SA Head</td><td>whirlin isoform x1</td><td>11.05</td><td>2.58E-13</td><td>4.25E-11</td></tr><tr><td>TR52871-c3_g1_i2</td><td>Ae-SA Head</td><td>carotenoid isomerooxygenase-like</td><td>10.91</td><td>6.65E-13</td><td>1.02E-10</td></tr><tr><td>TR5656-c4_g1_i1</td><td>Ae-SA Head</td><td>nuclear protein localization protein 4 homolog isoform x2</td><td>10.66</td><td>3.34E-12</td><td>4.52E-10</td></tr><tr><td>TR19336-c1_g1_i10</td><td>Ae-SA Head</td><td>fh1 fh2 domain-containing protein 3 isoform x4</td><td>10.26</td><td>4.36E-11</td><td>4.98E-09</td></tr><tr><td>TR41810-c12_g5_i10</td><td>Ae-SA Head</td><td>vinculin isoform x9</td><td>9.99</td><td>1.02E-14</td><td>2.40E-12</td></tr><tr><td>TR30823-c0_g1_i1</td><td>Ae-SA Head</td><td>zinc finger protein 853-like isoform x1</td><td>9.93</td><td>3.34E-10</td><td>3.24E-08</td></tr><tr><td>TR30823-c0_g1_i9</td><td>Ae-SA Head</td><td>zinc finger protein rotund isoform x7</td><td>9.86</td><td>5.27E-10</td><td>4.90E-08</td></tr><tr><td>TR53880-c6_g1_i1</td><td>Ae-SA Head</td><td>elks rab6-interacting cast family member 1 isoform x2</td><td>9.76</td><td>1.87E-07</td><td>1.27E-05</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 5** GO term distribution—biological process for genes with different expression patterns between *Aphidius ervi* libraries.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6** GO term distribution—Molecular function for genes with different expression patterns between *Aphidius ervi* libraries.</figcaption></figure>
<table><caption>**Table 3** Expression patterns of transcripts involved in chemoperception between *Aphidius ervi* populations reared on their natal hosts.</caption><thead><tr><td>Putative annotation</td><td>Log2-fold change</td><td>FDR-adjusted p-value</td><td>Higher in</td></tr></thead><tbody><tr><td>Odorant Receptor 13a like</td><td>4.13</td><td>2.26E-003</td><td>Ae-AP</td></tr><tr><td>Odorant Receptor 20</td><td>8.03</td><td>1.45E-005</td><td>Ae-AP</td></tr><tr><td>Odorant Receptor 98</td><td>3.24</td><td>3.37E-003</td><td>Ae-AP</td></tr><tr><td>Odorant Receptor OR2-like</td><td>8.37</td><td>3.22E-004</td><td>Ae-AP</td></tr><tr><td>Odorant Receptor 13 -isoform x2</td><td>3.76</td><td>5.98E-005</td><td>Ae-AP</td></tr><tr><td>Odorant Binding protein 83</td><td>4.40</td><td>3.78E-004</td><td>Ae-AP</td></tr><tr><td>Odorant Binding protein 69</td><td>3.36</td><td>2.34E-004</td><td>Ae-AP</td></tr><tr><td>Glutamate receptor kainate 2</td><td>3.75</td><td>2.10E-003</td><td>Ae-AP</td></tr><tr><td>Odorant Binding protein 56d like</td><td>7.63</td><td>2.10E-003</td><td>Ae-SA</td></tr></tbody></table>

adaptive evolution (Arya et al., 2015) or phenotypic plasticity. The identification of genes potentially involved in olfactory behavior in *A. ervi* and the highlighted gene expression differences between wasp populations exploiting different aphid lineages provide the first clues to understand the molecular basis of host-fidelity (Li et al., 2013; Zepeda-Paulo et al., 2013; Sepúlveda et al., 2017b). Homology analyses using NR database identified 91 contigs belonging to gene families involved in insect chemoperception such as OBPs (10 transcripts), CSPs (2 transcripts), SNMPs (1 transcript), ORs (76 transcripts, including the conserved odorant co-receptor, Orco) and ionotropic receptors (IRs; 2 transcripts). Compared to other parasitoids in which these gene families had been annotated using genomic approaches or transcriptomic sequencing from antennal tissues, we found less ORs (76 ORs, 10 OBPs, 2 CSPs, 1 SNMP) compared to *Cotesia chilonis* (117 ORs, 8 OBPs, 2 CSPs, 3 SNMPs) (Qi et al., 2015), but more than *Sclerodermus sp.* (8 ORs, 10 OBPs, 10 CSPs, 2 SNMPs) (Zhou et al., 2015). However, these numbers have to be taken with caution because these gene families are extremely difficult to annotate automatically. We found several significantly differentially expressed genes involved in chemical perception between our three *A. ervi* populations, two OBPs and five ORs were expressed at a higher level in the Ae-AP host-race while only one OBP had a higher expression value on the Ae-SA host-race (Table 3). Additionally, we found higher expression of an IR (the glutamate receptor kainate 2) in the Ae-AP host-race relative to Ae-SA host-race (Table 3). Interestingly, increased expression levels of glutamate receptor kainate 2 have been linked to olfactory responses in the salmon louse *Caligus rogercresseyi* (Núñez Acuña et al., 2014).

## Cellular signaling and neural development

Within the *A. ervi* reference transcriptome, genes coding for proteins participating in neuronal development and synaptic function of the nervous system were also found, which include the Rho family of GTPases. As Rho signaling activity plays a key role in neural morphological plasticity through dendritic reorganization and structural remodeling of synapses (Tolias, Duman & Um, 2011), it has been linked to long-term memory formation, including olfactory learning (Dobrin & Fahrbach, 2012). Additionally, variation in olfactory
behavior showed by different lines of *Drosophila melanogaster* has been linked to variants in genes involved in nervous systems development and function, such as Rho signaling (Arya et al., 2015).

Remarkably, we found a higher expression for Rho GTPase activating protein 190 (Rho-GAP 190) (Table 2) and Rho guanine nucleotide exchange factor 7 (Rho-GEF 7) (File S2), in the heads of Ae-AP compared to Ae-SA heads. We also found a higher expression for a gene coding for a sodium channel (sodium channel protein 60e) in Ae-AP heads compared to Ae-SA heads (File S1). In *D. melanogaster*, this sodium channel participates in processing olfactory information and regulates olfactory acuity, so its reduced expression impairs olfaction (Kulkarni et al., 2002; Zhang et al., 2013). Furthermore, this sodium channel in *D. melanogaster* is particularly expressed in olfactory organs (third antennal segment and maxillary palps) and brains (Gosselin-Badaroudine et al., 2016; Kulkarni et al., 2002). Other protein involved in cellular signaling is Calmodulin (CaM), a highly-conserved protein that contains four EF-hand domains that allow binding Ca<sup>+2</sup> ions (Park et al., 2008). Conformational changes in CaM allow its interaction with several target proteins including Orco (through Orco's CaM binding motif) (Bahk & Jones, 2016) and modulate insect OR function (Mukunda et al., 2014); repetitive subthreshold odor stimulation of olfactory neurons sensitizes ORs, and inhibition of CaM expression abolishes sensitization (Mukunda et al., 2016). Higher expression for a transcript coding for a CaM protein was found in Ae-AP heads compared to Ae-SA heads (neo-calmodulin-like isoform x4; 4.76 logFC; File S1).

DISCUSSION

Our study provides the first comprehensive reference transcriptome for the aphid parasitoid wasp *Aphidius ervi* obtained from females reared on three different aphid hosts: two host-races of the pea aphid *A. pisum* and the grain aphid *S. avenae*. From the RNA sequencing of head and body of these parasitoid lines, we were able to compare expression profiles and identify putative genes involved in host fidelity (*Sepúlveda et al., 2017b*). Enriched GO terms in a set of significantly up-regulated genes, suggest possible gene regulatory networks responsible for the observed differences at the phenotypic level. Not surprisingly, a comparison between head and body (thorax/gaster) tissues revealed that the activity related with stimulus perception and processing (signal transduction, transmembrane transport activity among others; Fig. 4.) are highly enriched in head tissue. On the other hand, GO terms enriched in bodies such as peptidase activity (i.e Neprilysin) are probably associated with the presence of transcripts coding for venom proteins, as *A. ervi* females have glands for venom production that is injected at oviposition and enriched in peptidases (*Colinet et al., 2014; Falabella et al., 2007; Nguyen et al., 2013*).

Behavioral experiments using the same laboratory populations from which the individuals of our transcriptome experiment were taken, showed that *A. ervi* collected from different aphid species (*A. pisum* and *S. avenae*) and host-races of the pea aphids (pea versus alfalfa) differed in several infectivity traits (host preference) but not in virulence (fitness) (Zepeda-Paulo et al., 2013). This is despite the lack of any detectable genetic structuring in Chilean *A. ervi* natural populations collected from different aphid hosts
species in the field (Daza-Bustamante et al., 2002; Zepeda-Paulo et al., 2015). The absence of genetic differentiation is explained by the very recent (1970ies) and single introduction event of *A. ervi* as a biological control agent for aphid pests in Chile (Zepeda-Paulo et al., 2016; Zúñiga et al., 1986). The proximal mechanisms of *A. ervi* infectivity (locating, searching and accepting an aphid host) most likely involve chemical cues from either the host-plant complex such as blends of host-induced volatiles produced by parasitized plants (Sasso et al., 2007; Takemoto & Takabayashi, 2015) or directly from the host (e.g., cuticle, cornicle secretions, faeces, exuviae, sex pheromones) (Powell & Wright, 1988). These cues are most likely learned during parasitoid development inside the host (Henry, Roitberg & Gillespie, 2008), during parasitoid emergency (Gutiérrez-Ibáñez, Villagra & Niemeyer, 2007) or imprinted rather than based on genotyping differences in preference, i.e., innate (Antolin, Bjorksten & Vaughn, 2006). Although the influence of possible genotypic variation between *A. ervi* populations can currently not be ruled out, as there could be a genetic basis accounting for variation in olfactory behavior as is the case of *D. melanogaster* (Arya et al., 2015), further studies need to be carried out to determine whether variations in olfactory behavior are inherited in *A. ervi* or they are just learned based on phenotypic plasticity and differential gene expression.

We know from many examples that the exposure of insects to different environments during juvenile development can lead to substantial differences in the transcriptome of adult individuals (Berens, Hunt & Toth, 2015; Schrader et al., 2015). Hence, in our comparison between head transcriptomes of *A. ervi* reared on different aphid species we focused on differential expression of key components of olfactory and learning pathways (e.g., peripheral system: olfactory receptors, odorant binding proteins, ionotropic receptors, gustatory receptors; nervous system function and development). We identified in silico a total of 91 unigenes possessing high-sequence identities with chemosensation-related genes, including IRs, ORs, OBPs, CSPs, SNMPs and Orco, which is similar to what has been reported for other parasitoids species (Qi et al., 2015; Zhou et al., 2015), so our approach proved to be effective for detecting and annotating genes coding for proteins associated with olfaction. As changes in olfactory sensitivity could be driven by mutations in key genes, gene gains and losses, and/or variation in gene expression (Glaser et al., 2015), and due to the lack of genetic differences between parasitoids coming from different hosts, we focused our analysis on gene expression. We found 3 OBPs and 5 ORs that were differentially expressed between head samples of Ae-SA and Ae-AP. Interestingly, both Ae-AP populations (APA and APP) had a higher number of up-regulated chemosensory genes and neuronal-related genes compared to Ae-SA (Table 3–additional file 1), which is similar to that observed in other insect species capable of using different hosts, such as in the Mediterranean corn borer *Sesamia nonagrioides* (Glaser et al., 2015).

We also studied the expression profiles for genes involved in both cellular signaling and neural development. Surprisingly, genes coding for proteins involved in neuronal morphology re-modelling (Tolias, Duman & Um, 2011) had higher transcriptional levels in Ae-AP heads, including Rho-GTPase activating proteins and Rho guanine nucleotide exchange factor, both of which have been described as regulators of Rho signaling activity and have been linked to olfactory learning and long-term memory formation in *A. mellifera*
(Dobrin & Fahrbach, 2012), while variants in genes involved in neural development and signaling (e.g., Rho proteins) are related to variation in olfactory behavior in *D. melanogaster* (Arya et al., 2015). Regarding signal transduction, we found higher expression levels in Ae-AP heads for both the sodium channel protein 60e, which participates in processing olfactory information and olfactory acuity (Kulkarni et al., 2002; Zhang et al., 2013), and Calmodulin (CaM), which interacts with Orco through binding of Ca<sup>+2</sup> ions and modulates insect OR function (Mukunda et al., 2014).

Taken together, our transcriptional evidence coincides with the observation that *A. ervi*—AP shows a higher discrimination in terms of host preference than *A. ervi*—SA (Zepeda-Paulo et al., 2013). This differential expression of candidate chemosensory genes, signaling genes and neuronal development genes may explain host-fidelity in *A. ervi*, and could be a signature of adaptive phenotypic plasticity to different host and host-plant induced environments (Glaser et al., 2015). The genes emphasized in this study deserve special attention for future research in order to prove their role in host preference and host selection by *A. ervi* (i.e., host fidelity). Further studies should also consider detailed analysis both at gene sequence level between *A. ervi* populations, and alternative splicing of these coding genes, as in the case of the tarnished plant bug (*Lygus lineolaris*), where alternative splicing may contribute to the divergence of OBPs (Oppenheim et al., 2015).

More generally, our results provide some input for discussion on the impact of phenotypic plasticity on evolutionary changes. It is indeed puzzling that, although the population does not seem to be subdivided but rather genetically homogenous (most likely due to the relative recent introduction of a limited number of individuals into Chile), we can observe distinct and persistent phenotypic and transcriptional differences between parasitoids coming from different aphid host species and host races. It has been put forward that a developmental reorganization from ancestral phenotypes due to new environmental input or conditions (host races or host biotypes in this case) does not necessarily require new mutations to produce novel or distinct phenotypes (West-Eberhard, 2003). The *A. ervi* populations in Chile may provide an example for developmental reorganization because it seems unlikely that they evolved novel adaptive mutations, as this biocontrol agent was intentionally introduced in few numbers in Chile just about 40 years ago. The next evolutionary step towards speciation based on host races would/could be genetic accommodation, i.e., developmental variants (in our case host races). When novel genetic variants are fixed within populations, those variants could limit phenotypic plasticity, so individuals can no longer easily switch between different phenotypes due to trade-offs and constraints (i.e., reduced fitness if they are not reared on their preferred hosts). It would be interesting to see whether and how many of the differentially expressed genes underlying the observed host related differences in *A. ervi* also show fixed genetic differences.

## ACKNOWLEDGEMENTS

We thank Francisca Zepeda-Paulo for providing *A. ervi* populations, Sebastián Ortiz-Martinez and Daniela Sepúlveda for help in aphids and parasitoid rearing, Cinthya Villegas for help in parasitoid rearing and sexual identification, Lucia Briones for help in laboratory RNA extractions and Ruben Palma for help in data analysis.
# ADDITIONAL INFORMATION AND DECLARATIONS

## Funding

This work was funded by FONDECYT grant 1130483 to CC Figueroa and Iniciativa Científica Milenio (ICM) NC120027. GI Ballesteros was funded by CONICYT PhD fellowship 21120371. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## Grant Disclosures

The following grant information was disclosed by the authors:

FONDECYT grant: 1130483.

Iniciativa Científica Milenio (ICM): NC120027.

CONICYT PhD fellowship: 21120371.

## Competing Interests

The authors declare there are no competing interests.

## Author Contributions

* Gabriel I. Ballesteros conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, wrote the paper, prepared figures and/or tables, reviewed drafts of the paper.

* Jürgen Gadau conceived and designed the experiments, analyzed the data, contributed reagents/materials/analysis tools, wrote the paper, prepared figures and/or tables, reviewed drafts of the paper.

* Fabrice Legeai contributed reagents/materials/analysis tools, reviewed drafts of the paper.

* Angelica Gonzalez-Gonzalez reviewed drafts of the paper, bioinformatics support.

* Blas Lavandero A. *ervi* collections/contributed reagents/materials/analysis tools, reviewed drafts of the paper.

* Jean-Christophe Simon reviewed drafts of the paper.

* Christian C. Figueroa conceived and designed the experiments, analyzed the data, contributed reagents/materials/analysis tools, wrote the paper, reviewed drafts of the paper.

## DNA Deposition

The following information was supplied regarding the deposition of DNA sequences:
Transcriptomic assembly is available with GenBank accession GFLW01000000. Raw data is available under BioProject ID: PRJNA377544. Filtered transcriptomic assembly is available at: Ballesteros, Gabriel (2017): *Aphidius ervi, de novo* transcriptome assembly, reduced set. figshare.

https://doi.org/10.6084/m9.fishshare.4816939.v2.

## Data Availability

The following information was supplied regarding data availability:
Raw transcriptomic data is available under BioProject ID: PRJNA377544.
Table S1 is available at Ballesteros, Gabriel (2017): differentially expressed gene lists for
*A. ervi*. figshare.

https://doi.org/10.6084/m9 figshare.4816942.v4.

Table S2 and Fig. S1 are available at: Ballesteros, Gabriel (2017): Read count table mapped and sample correlation matrix—*A. ervi* transcriptome study. figshare.

https://doi.org/10.6084/m9 figshare.5116270.v1.

Fasta file of filtered transcriptome assembly is available at: Ballesteros, Gabriel (2017):
*Aphidius ervi*, de novo transcriptome assembly, reduced set. figshare.

https://doi.org/10.6084/m9 figshare.4816939.v2.

Annotation of transcriptome is available at: Ballesteros, Gabriel (2017): Annotation
table for *A. ervi* transcriptome, filtered dataset. figshare.

https://doi.org/10.6084/m9 figshare.4822069.v1.

**Supplemental Information**

Supplemental information for this article can be found online at http://dx.doi.org/10.7717/
peerj.3640#supplemental-information.

REFERENCES

**Abrahamson WG, Blair CP. 2008.** Sequential radiation through host-race formation: herbivore diversity leads to diversity in natural enemies. In: *Specialization, speciation, and radiation: the evolutionary biology of herbivorous insects*. Berkeley: University of California Press, 188–202.

**Antolin MF, Bjorksten TA, Vaughn TT. 2006.** Host-related fitness trade-offs in a presumed generalist parasitoid, *Diaeretiella rapae* (Hymenoptera: Aphidiidae).
*Ecological Entomology* **31**:242–254 DOI 10.1111/j.1365-2311.2006.00769.x.

**Arya GH, Magwire MM, Huang W, Serrano-Negron YL, Mackay TFC, Anholt RRH. 2015.** The genetic basis for variation in olfactory behavior in *Drosophila melanogaster*. *Chemical Senses* **40**:233–243 DOI 10.1093/chemse/bjv001.

**Bahk S, Jones WD. 2016.** Insect odorant receptor trafficking requires calmodulin. *BMC Biology* **14**:83 DOI 10.1186/s12915-016-0306-x.

**Berens AJ, Hunt JH, Toth AL. 2015.** Nourishment level affects caste-related gene expression in *Polistes* wasps. *BMC Genomics* **16**:235 DOI 10.1186/s12864-015-1410-y.

**Bolger AM, Lohse M, Usadel B. 2014.** Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics* **30**:2114–2120 DOI 10.1093/bioinformatics/btu170.

**Burke GR, Strand MR. 2012.** Deep sequencing identifies viral and wasp genes with potential roles in replication of *Microplitis demolitor* Bracovirus. *Journal of Virology* **86**:3293–3306 DOI 10.1128/JVI.06434-11.

**Calla B, Sim SB, Hall B, Derego T, Liang GH, Geib SM. 2015.** Transcriptome of the egg parasitoid *Fopius arisanus*: an important biocontrol tool for Tephritid fruit fly suppression. *Gigascience* **4**:1–5 DOI 10.1186/s13742-015-0075-4.

**Cameron PJ, Powell W, Loxdale HD. 1984.** Reservoirs for *Aphidius ervi* Hali-day (Hymenoptera: Aphidiidae), a polyphagous parasitoid of cereal aphids
(Hemiptera: Aphididae). *Bulletin of Entomological Research* **74**:647–656
DOI 10.1017/S0007485300014024.

Colinet D, Anselme C, Deleury E, Mancini D, Poulain J, Azéma-Dossat C, Belghazi M, Tares S, Pennacchio F, Poirié M, Gatti J. 2014. Identification of the main venom protein components of *Aphidius ervi*, a parasitoid wasp of the aphid model *Acyrthosiphon pisum*. *BMC Genomics* **15**:342 DOI 10.1186/1471-2164-15-342.

Conesa A, Götz S. 2008. Blast2GO: a comprehensive suite for functional analysis in plant genomics. *International Journal of Plant Genomics* **2008**:619832 DOI 10.1155/2008/619832.

Crispo E. 2008. Modifying effects of phenotypic plasticity on interactions among natural selection, adaptation and gene flow. *Journal of Evolutionary Biology* **21**:1460–1469 DOI 10.1111/j.1420-9101.2008.01592.x.

Daza-Bustamante P, Fuentes-Contreras E, Niemeyer HM. 2003. Acceptance and suitability of *Acyrthosiphon pisum* and *Sitobion avenae* as hosts of the aphid parasitoid *Aphidius ervi* (Hymenoptera: Braconidae). *European Journal of Entomology* **100**:49–53 DOI 10.14411/eje.2003.010.

Daza-Bustamante P, Fuentes-Contreras E, Rodriguez LC, Figueroa CCC, Niemeyer HMM, Rodríguez LC, Figueroa CC, Niemeyer HM. 2002. Behavioural differences between *Aphidius ervi* populations from two tritrophic systems are due to phenotypic plasticity. *Entomologia Experimentalis et Applicata* **104**:321–328 DOI 10.1046/j.1570-7458.2002.01018.x.

Dobrin SE, Fahrbach SE. 2012. Rho GTPase activity in the honey bee mushroom bodies is correlated with age and foraging experience. *Journal of Insect Physiology* **58**:228–234 DOI 10.1016/j.jinsphys.2011.11.009.

Falabella P, Riviello L, Caccialupi P, Rossodivita T, Teresa Valente M, Luisa De Stradis M, Tranfaglia A, Varriccio P, Gigliotti S, Graziani F, Malva C. 2007. A γ-glutamyl transpeptidase of *Aphidius ervi* venom induces apoptosis in the ovaries of host aphids. *Insect Biochemistry and Molecular Biology* **37**:453–465 DOI 10.1016/j.ibmb.2007.02.005.

Glaser N, Gallot A, Legeai F, Harry M, Kaiser L, Le Ru B, Calatayud P-A, Jacquin-Joly E. 2015. Differential expression of the chemosensory transcriptome in two populations of the stemborer *Sesamia nonagrioides*. *Insect Biochemistry and Molecular Biology* **65**:28–34 DOI 10.1016/j.ibmb.2015.07.008.

Gosselin-Badaroudine P, Moreau A, Simard L, Cens T, Rousset M, Collet C, Charnet P, Chahine M. 2016. Biophysical characterization of the honeybee DSC1 orthologue reveals a novel voltage-dependent Ca²⁺ channel subfamily: CaV 4. *Journal of General Physiology* **148**:Jgp.201611614 DOI 10.1085/jgp.201611614.

Gurevich A, Saveliev V, Vyahhi N, Tesler G. 2013. QUAST: quality assessment tool for genome assemblies. *Bioinformatics* **29**:1072–1075 DOI 10.1093/bioinformatics/btt086.

Gutiérrez-Ibáñez C, Villagra CA, Niemeyer HM. 2007. Pre-pupation behaviour of the aphid parasitoid *Aphidius ervi* (Haliday) and its consequences for pre-imaginal learning. *Naturwissenschaften* **94**:595–600 DOI 10.1007/s00114-007-0233-3.
Henry LM, May N, Acheampong S, Gillespie DR, Roitberg BD. 2010. Host-adapted par- asitoids in biological control: does source matter? *Ecological Applications* 20:242–250
DOI 10.1890/08-1869.1.

Henry LM, Roitberg BD, Gillespie DR. 2008. Host-range evolution in *Aphidius* par- asitoids: fidelity, virulence and fitness trade-offs on an ancestral host. *Evolution* 62:689–699 DOI 10.1111/j.1558-5646.2007.00316.x.

Hwang H-J, Patnaik BB, Kang SW, Park SY, Wang TH, Park EB, Chung JM, Song DK, Patnaik HH, Kim C, Kim S, Lee JB, Jeong HC, HS Park, Han YS, Lee YS. 2016. RNA sequencing, *de novo* assembly, and functional annotation of an endangered Nymphalid butterfly, *Fabriciana nerippe* Felder, 1862. *Entomological Research* 46(2):148–161 DOI 10.1111/1748-5967.12164.

Jones TS, Bilton AR, Mak L, Sait SM. 2015. Host switching in a generalist parasitoid: contrasting transient and transgenerational costs associated with novel and original host species. *Ecology and Evolution* 5:459–465 DOI 10.1002/ece3.1333.

Kulkarni NH, Yamamoto AH, Robinson KO, Mackay TFC, Anholt RRH. 2002. The DSC1 channel, encoded by the smi60E locus, contributes to odor-guided behavior in *Drosophila melanogaster*. *Genetics* 161:1507–1516.

Langmead B, Trapnell C, Pop M, Salzberg SL. 2009. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. *Genome Biology* 10:1–10 DOI 10.1186/gb-2009-10-3-r25.

Li B, Dewey CN. 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. *BMC Bioinformatics* 12:323 DOI 10.1186/1471-2105-12-323.

Li H, Zhang H, Guan R, Miao X. 2013. Identification of differential expression genes associated with host selection and adaptation between two sibling insect species by transcriptional profile analysis. *BMC Genomics* 14:582 DOI 10.1186/1471-2164-14-582.

Loxdale HD, Harvey JA. 2016. The “generalism” debate: misinterpreting the term in the empirical literature focusing on dietary breadth in insects. *Biological Journal of the Linnean Society* 119:265–282 DOI 10.1111/bij.12816.

Moreton J, Izquierdo A, Emes RD. 2016. Assembly, assessment, and availability of *De novo* generated eukaryotic transcriptomes. *Frontiers in Genetics* 6:1–9 DOI 10.3389/fgene.2015.00361.

Mukunda L, Miazzi F, Kaltofen S, Hansson BS, Wicher D. 2014. Calmodulin modulates insect odorant receptor function. *Cell Calcium* 55:191–199 DOI 10.1016/j.ceca.2014.02.013.

Mukunda L, Miazzi F, Sargsyan V, Hansson BS, Wicher D. 2016. Calmodulin affects sensitization of *Drosophila melanogaster* odorant receptors. *Frontiers in Cellular Neuroscience* 10:1–11 DOI 10.3389/fncel.2016.00028.

Nguyen TTA, Magnoli I, Cloutier C, Michaud D, Muratori F, Hance T. 2013. Early presence of an enolase in the oviposition injecta of the aphid parasitoid *Aphidius ervi* analyzed with chitosan beads as artificial hosts. *Journal of Insect Physiology* 59:11–18 DOI 10.1016/j.jinsphys.2012.11.002.
Núñez Acuña G, Valenzuela-Muñoz V, Marambio JP, Wadsworth S, Gallardo-Escárate C. 2014. Insights into the olfactory system of the ectoparasite *Caligus rogercresseyi*: molecular characterization and gene transcription analysis of novel ionotropic receptors. *Experimental Parasitology* 145:99–109 DOI 10.1016/j.exppara.2014.08.003.

Ono H, Ishii K, Kozaki T, Ogiwara I, Kanekatsu M, Yamada T. 2015. Removal of redundant contigs from *de novo* RNA-Seq assemblies via homology search improves accurate detection of differentially expressed genes. *BMC Genomics* 16:1031 DOI 10.1186/s12864-015-2247-0.

Oppenheim SJ, Baker RH, Simon S, DeSalle R. 2015. We can't all be supermodels: the value of comparative transcriptomics to the study of non-model insects. *Insect Molecular Biology* 24:139–154 DOI 10.1111/imb.12154.

Park HY, Kim SA, Korlach J, Rhoades E, Kwok LW, Zipfel WR, Waxham MN, Webb W.W, Pollack L. 2008. Conformational changes of calmodulin upon Ca2+ binding studied with a microfluidic mixer. *Proceedings of the National Academy of Sciences of the United States of America* 105:542–547 DOI 10.1073/pnas.0710810105.

Peccoud J, Bonhomme J, Mahéo F, De la Huerta M, Cosson O, Simon JC. 2014. Inheritance patterns of secondary symbionts during sexual reproduction of pea aphid biotypes. *Insect Science* 21(3):291–300 DOI 10.1111/1744-7917.12083.

Peccoud J, Figueroa CC, Silva AX, Ramirez CC, Mieuzet L, Bonhomme J, Stoeckel S, Plantegenest M, Simon JC. 2008. Host range expansion of an introduced insect pest through multiple colonizations of specialized clones. *Molecular Ecology* 17:4608–4618 DOI 10.1111/j.1365-294X.2008.03949.x.

Powell W, Wright AF. 1988. The abilities of the aphid parasitoids *Aphidius ervi* Haliday and *A. rhopalosiphi* De Stefani Perez (Hymenoptera: Braconidae) to transfer between different known host species and the implications for the use of alternative hosts in pest control strategies. *Bulletin of Entomological Research* 78:683–693 DOI 10.1017/S0007485300015546.

Qi Y, Teng Z, Gao L, Wu S, Huang J, Ye G, Fang Q. 2015. Transcriptome analysis of an endoparasitoid wasp *Cotesia chilonis* (Hymenoptera: Braconidae) reveals genes involved in successful parasitism. *Archives of Insect Biochemistry and Physiology* 88:203–221 DOI 10.1002/arch.21214.

Rivero A. 2000. The relationship between host selection behaviour and offspring fitness in a koinobiont parasitoid. *Ecological Entomology* 25:467–472 DOI 10.1046/j.1365-2311.2000.00276.x.

Sasso R, Iodice L, Cristina Digilio M, Carretta A, Ariati L, Guerrieri E. 2007. Host-locating response by the aphid parasitoid *Aphidius ervi* to tomato plant volatiles. *Journal of Plant Interactions* 2:175–183 DOI 10.1080/17429140701591951.

Schluter D. 2000. *The ecology of adaptive radiation*. 1st edition. Oxford: Oxford University Press.

Schrader L, Simola DF, Heinze J, Oettler J. 2015. Sphingolipids, transcription factors, and conserved toolkit genes: developmental plasticity in the ant cardiocondyla obscurior. *Molecular Biology and Evolution* 32:1474–1486 DOI 10.1093/molbev/msv039.
**Sepúlveda DA, Zepeda-Paulo F, Ramírez CC, Lavandero B, Figueroa CC. 2017a.** Diversity, frequency and geographic distribution of facultative bacterial endosymbionts in introduced aphid pests. *Insect Science* **24**:511–521 DOI 10.1111/1744-7917.12313.

**Sepúlveda DA, Zepeda-Paulo F, Ramírez CC, Lavandero B, Figueroa CC. 2017b.** Loss of host fidelity in highly inbred populations of the parasitoid wasp *Aphidius ervi* (Hymenoptera: Braconidae). *Journal of Pest Science* **90**:649–658 DOI 10.1007/s10340-016-0798-8.

**Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. 2015.** BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. *Bioinformatics* **31**:3210–3212 DOI 10.1093/bioinformatics/btv351.

**Starý P. 1993.** The fate of released parasitoids (Hymenoptera: Braconidae, Aphidiinae) for biological control of aphids in Chile. *Bulletin of Entomological Research* **83**:633–639 DOI 10.1017/S0007485300040062.

**Starý P. 1995.** The Aphidiidae of Chile (Hymenoptera, Ichneumonoidea, Aphidiidae). *Deutsche Entomologische Zeitschrift* **42**:113–138.

**Starý P, Gerding M, Norambuena H, Remaudiere G. 1993.** Environmental-Research on aphid parasitoid biocontrol agents in Chile (Hym, Aphidiidae, Hom, Aphidoidea). *Journal of Applied Entomology* **115**:292–306 DOI 10.1111/j.1439-0418.1993.tb00394.x.

**Stireman JO, Nason JD, Heard SB, Seehawer JM. 2006.** Cascading host-associated genetic differentiation in parasitoids of phytophagous insects. *Proceedings of the Royal Society B: Biological Sciences* **273**:523–530 DOI 10.1098/rspb.2005.3363.

**Stilmant D, Bellinghen C, Hance T, Boivin G. 2008.** Host specialization in habitat specialists and generalists. *Oecologia* **156**:905–912 DOI 10.1007/s00442-008-1036-8.

**Suh E, Bohbot JD, Zwiebel LJ. 2014.** Peripheral olfactory signaling in insects. *Current Opinion in Insect Science* **6**:86–92 DOI 10.1016/j.cois.2014.10.006.

**Takemoto H, Kainoh Y, Takabayashi J. 2011.** Learning of plant volatiles by aphid parasitoids: timing to learn. *Journal of Plant Interactions* **6**:137–140 DOI 10.1080/17429145.2010.548607.

**Takemoto H, Takabayashi J. 2015.** Parasitic wasps *Aphidius ervi* are more attracted to a blend of host-induced plant volatiles than to the independent compounds. *Journal of Chemical Ecology* DOI 10.1007/s10886-015-0615-5.

**Tassone EE, Geib SM, Hall B, Fabrick JA, Brent CS, Hull JJ. 2016.** *De novo* construction of an expanded transcriptome assembly for the western tarnished plant bug, *Lygus hesperus*. *Gigascience* **5**:1–6 DOI 10.1186/s13742-016-0109-6.

**Tolias KF, Duman JG, Um K. 2011.** Control of synapse development and plasticity by Rho GTPase regulatory proteins. *Progress in Neurobiology* **94**:133–148 DOI 10.1016/j.pneurobio.2011.04.011.

**Wang XG, Kaçar G, Biondi A, Daane KM. 2016.** Life-history and host preference of *Trichopria drosophilae*, a pupal parasitoid of spotted wing drosophila. *BioControl* **61**:387–397 DOI 10.1007/s10526-016-9720-9.

**West-Eberhard MJ. 2003.** *Developmental plasticity and evolution*. 1st edition. New York: Oxford University Press.
Zepeda-Paulo F, Dion E, Lavandero B, Mahéo F, Outreman Y, Simon JC, Figueroa CC. 2016. Signatures of genetic bottleneck and differentiation after the introduction of an exotic parasitoid for classical biological control. *Biological Invasions* 18:565–581 DOI 10.1007/s10530-015-1029-6.

Zepeda-Paulo F, Lavandero B, Mahéo F, Dion E, Outreman Y, Simon JC, Figueroa CC. 2015. Does sex-biased dispersal account for the lack of geographic and host-associated differentiation in introduced populations of an aphid parasitoid? *Ecology and Evolution* 5:2149–2161 DOI 10.1002/ece3.1504.

Zepeda-Paulo FA, Ortiz-Martínez SA, Figueroa CC, Lavandero B. 2013. Adaptive evolution of a generalist parasitoid: Implications for the effectiveness of biological control agents. *Evolutionary Applications* 6:983–999 DOI 10.1111/eva.12081.

Zhang T, Wang Z, Wang L, Luo N, Jiang L, Liu Z, Wu CF, Dong K. 2013. Role of the DSC1 channel in regulating neuronal excitability in *Drosophila melanogaster*: extending nervous system stability under stress. *PLOS Genetics* 9(3):e1003327 DOI 10.1371/journal.pgen.1003327.

Zhou C-X, Min S-F, Yan-Long T, Wang M-Q. 2015. Analysis of antennal transcriptome and odorant binding protein expression profiles of the recently identified parasitoid wasp, *Sclerodermus* sp. *Comparative Biochemistry and Physiology - Part D: Genomics and Proteomics* 16:10–19 DOI 10.1016/j.cbd.2015.06.003.

Zúñiga E, Van den Bosch R, Drea J, Gruber F. 1986. The biological control project against the cereal aphids (Hom., Aphididae) in Chile. 2: exploration, importation and quarantine of predator and parasitoid species. Agric. Tec. *Agricultura Tecnica* 46:479–487.