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

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 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 1^e−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.

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 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 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 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 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 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 olfactorybehavior 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⁺² 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).