Transcriptomics of an extended phenotype: parasite manipulation of wasp social behaviour shifts expression of caste-related genes

1. Background

Parasitism is a widespread and highly successful life-history strategy. Parasites can reduce the fitness of their hosts, often causing profound changes in host physiology and behaviour [1–3]. In some cases, referred to as ‘adaptive host manipulation’, phenotypic changes in the host may benefit the parasite [4]. There are several striking cases in which parasites induce unusual and novel behaviours that are not typically part of the host's normal behavioural repertoire [5]. Alternatively, parasites may cause more subtle alterations to host behaviour, such as changes in frequency or timing of otherwise normal behaviour [6,7]. There are many open questions about the mechanisms by which parasites accomplish host behavioural manipulation, and how novel parasite-induced host behaviours evolve. Advances in genomics offer an exciting opportunity to examine the mechanistic underpinnings and evolution of adaptive host manipulation [8].

Phenotypic plasticity is an important driver of evolutionary change and phenotypic novelty [9]. One strategy a parasite might take to shift its host's phenotype would be to exploit existing host phenotypic plasticity (and associated genetic pathways) in ways that benefit the parasite. This has the potential to be an efficient strategy for parasite manipulation of host phenotypes. Previous studies have documented a wide array of host behavioural, physiological and gene expression changes associated with parasitization [2,10–15]. However, few studies have examined changes in host gene expression through the lens of adaptive parasite manipulation, also viewed as an ‘extended phenotype’ of the parasite's genome [16]. Because phenotypic plasticity is the result of shifts in gene expression, one powerful way to examine how parasites affect their hosts is using transcriptomics, which is now easily applicable to ecological model systems [17]. Although previous studies have examined the effects of parasites on their hosts' transcriptomes [4,12–15], none have specifically examined how and if parasites influence host behaviour via hijacking inherent host phenotypic plasticity.

To begin to explore the hypothesis that parasites can manipulate hosts by shifting plastic phenotypes, we focus on transcriptome-wide gene expression in the well-studied case of paper wasp hosts, Polistes dominula, and their parasites, Xenos vesparum (Insecta: Strepsiptera). This system affords an exciting opportunity to study mechanisms and evolution of host manipulation for several reasons [18,19]. First, P. dominula colonies are often heavily parasitized by X. vesparum, which are obligate endoparasitoids, meaning that they must develop within the host and also lead to reproductive death via castration [20,21]. Second, P. dominula are a model species with substantial behavioural and developmental plasticity in the form of flexible social castes [22,23]. Third, X. vesparum have well-documented effects on their hosts, including behavioural and physiological changes that suggest they create ‘wrong caste’ phenotypes in hosts (described below) [18,19]. Fourth, P. dominula have new genomic data available, including a genome and caste-related transcriptome [18], and identified genes and pathways related to caste differences [24–29].

The typical colony cycle of temperate Polistes wasps involves several phases, and in each, there are distinct groups of females showing extensive phenotypic plasticity in physiology and behaviour (figure 1). Colonies are founded in the spring by one or a few cooperating founding females, who build the nest and rear a first set of offspring, and all females, who become workers upon adult emergence [30]. Upon worker emergence, the founding female becomes the primary reproductive (queen) of the colony, while the workers perform tasks related to brood care and colony growth [31]. Later in the colony cycle, larvae are reared by workers and emerge as males or female ‘gynes’—non-working pre-overwintering queens with large fat stores [32]. Gynes will leave the colony in fall to form extranidal aggregations with other gynes, where they overwinter until they disperse to found new colonies the following spring [30].

• Download figure
• Open in new tab
• Download PowerPoint

Xenos vesparum infest P. dominula larvae in various host larval instars in multiple colony phases [33], including early in the season when the founding queen is laying worker-destined eggs [33–37] (figure 1). Xenos vesparum remains endoparasitic until after the host emerges as an adult [36]. Then, X. vesparum extrudes its anterior region between the host's abdominal segments. The male parasitoid emerges as a fully motile, winged adult that will disperse to mate [21,33,38], whereas the female remains neotenic, living permanently endoparasitic, with an extruded cephalothorax with no eyes or appendages [20,21,39].

A parasitized female wasp (hereafter ‘stylopized female’, after suborder Stylopidia) undergoes a dramatic physio-behavioural deviation [40] (figure 1; summarized in electronic supplementary material, table S1), including loss of ovaries, asocial behaviour and prolonged lifespan, when infested. Early emerging stylopized females desert their natal colonies in mid-summer, during peak colony activity, and form extranidal aggregations with stylopized females from other colonies [19,41]. This aberrant host behaviour occurs close to the time of male parasitoid emergence and appears to help facilitate mate location by the parasitoids [38]. Later, stylopized females are joined by unstylopized aggregating gynes, forming overwintering clusters. Stylopized females remain with the unstylopized gynes until the next spring when the latter become nest-founding queens [19,40,41]. The stylopized females do not found nests, and instead, transmit parasitoid larvae by visiting other newly founded P. dominula nests or via shared floral forage sites [33,36,37].

This study focuses on the unusual nest desertion and aggregation behaviour of stylopized females. Because the stylopized females described here typically emerge in early summer, they should be workers, engaging in colony work such as nest building, foraging and brood care. However, they exhibit behaviour that is typical of unstylopized gynes—lack of nest work, nest desertion and formation of pre-overwintering aggregations. Despite our detailed understanding of the behavioural effects of stylopization, the molecular mechanisms underlying this behavioural deviation are unknown.

Here, we provide transcriptomic data underlying the extended phenotype of X. vesparum parasitoids in the context of manipulation of their P. dominula hosts' behaviour. We used high-throughput RNA-sequencing to measure global brain gene expression patterns in unstylopized workers, unstylopized gynes in pre-winter aggregations and early emerging stylopized females. Importantly, we measured gene expression in stylopized females before they deserted their nests; thus, we were able to capture gene expression patterns that preceded the full expression of the parasitoid-induced behaviour and are more likely to be causal to, rather than a consequence of, the shifted behaviour. First, we hypothesized that the overall pattern of brain gene expression in stylopized females would be shifted from towards gyne-like patterns. Second, we also hypothesized that the parasitoid-induced shift in gene expression would involve caste-related genes, i.e. those that show differential expression between unstylopized workers and gynes.

2. Material and methods

(a) Sample collection

Wasps were collected at three sites in Tuscany, Italy, including Siena, Impruneta and Sesto, in the spring and summer of 2009. Upon collection, wasp heads were immediately removed and placed in RNAlater^© (Life Technologies, AM7023), stored at 4°C for 24 h (as per the manufacturers' instructions), then transferred to −20°C and kept frozen until brain dissections. Bodies were frozen for later dissection. We collected three groups of females: 21 unstylopized workers from nine nests in June (herein referred to as ‘Workers, W’); 20 unstylopized gynes from pre-overwintering aggregations in August (Gynes, G) and 15 stylopized females from six nests in June (Stylopized, S). Stylopized females were collected from a subset of the same nests as workers.

(b) Dissections and RNA extractions

Abdomens were removed from the wasps' bodies and dissected. To identify caste status and parasitoid presence, the following information was recorded for each wasp: ovary development score [42], fat body size score [42] and presence of X. vesparum, including number, developmental stage and sex of the parasitoids [43]. We selected 15 wasps that were stylopized by at least one neotenic female X. vesparum (although some were stylopized by more than one female or a female and a male).

Brains from each individual were dissected in RNAlater^© at room temperature and then stored at −80°C until RNA extraction. Total RNA from each brain was extracted using the RNeasy^© Mini Kit (Qiagen, 74 104), including a 15 min DNase I (Qiagen, 79 254) step before ethanol washes. RNA was quantified using a NanoDrop 1000^© (Thermo Fisher Scientific) spectrophotometer and Qubit Fluorometer (Life Technologies). Samples with 260 : 230/260 : 280 < 1.9 were not included in the transcriptomic analysis. A subset of samples were quantified by Bioanalyzer (Agilent), which also served to verify quality based on RNA integrity numbers. Following quantification, RNA was stored at −80°C until later use. One set was allocated exclusively for transcriptome sequencing and another for qRT-PCR validation. RNA from 14 pooled samples were used for sequencing (five gynes, five workers and four stylopized). The composition of each pooled sample varied slightly, from two to four individual brains, to provide enough RNA for sequencing.

(c) RNA-sequencing

At the Pennsylvania State University Genomics Core Facility, each of the pooled RNA samples was used for a standard SOLiD mRNA library preparation with barcoding. We performed two rounds of oligo(dT) selection of the poly(A) RNA with the MicroPoly (A) Purist™ Kit (Ambion LifeTech). cDNA libraries were prepared starting from 100 to 500 ng poly(A) RNA and using the SOLiD Total RNA-Seq Kit according to the manufacturer's instructions (Applied Biosystems/Life Technologies). Sequencing was performed on the Applied Biosystems (Life Technologies) SOLiDv4 according to the manufacturer's instructions. All samples were loaded into a single flow cell of a SOLiD 5500xl Sequencer, producing 76 base single-end reads (total number of reads per sample listed in electronic supplementary material, table S2).

To minimize the noise from sequencing mistakes, raw reads were processed using the Perl script (SOLiD_preprocess_filter_v2.pl). Reads with quality scores below 10 were discarded. A comparison between two popular colour space alignment tools, BFAST and SHRiMP, led us to use BFAST for improved mapping rates (electronic supplementary material, table S2). Using BFAST, we first built two reference genomes both in colour space and in nucleotide space. Then, we created indexes for the reference genome in colour space and found candidate alignment locations (CALs). Next, local alignments were performed in CALs using the P. dominula reference genome [24] in nucleotide space. Finally, low-quality alignments were filtered by a post-processing step.

(d) Analysis of gene expression

We used HTSeq (http://www-huber.embl.de/users/anders/HTSeq/) for counting the number of reads that mapped to 11 506 transcripts based on the P. dominula genome annotation r1.0 [24]. These counts were used by edgeR for differential expression analyses. We filtered low expressed transcripts by the following criterion. We calculated transcript CPM (counts per million mapped reads) for each sample and preserved transcripts which satisfied CPM > 3 in at least three samples. By using 3 as a cut-off, 8484 (approx. 74%) transcripts remained for the differential expression analyses.

Raw read counts were analysed with R using the edgeR package from Bioconductor [44], using the standard protocol outlined in the edgeR user guide, except we used a more stringent CPM cut-off. The edgeR user guide recommends keeping a gene when at least two samples with CPM > 1. However, considering our relatively high sample size and read counts, we set our threshold for keeping a gene as at least three samples with CPM > 3. Normalization was performed using the TMM (Trimmed Mean of M-values) method. First, we wished to examine which of the three wasp phenotypes had distinct expression profiles versus all of the others, and also which phenotype explained most of the difference in global gene expression. For this purpose, we adopted a generalized linear model (GLM) fit to the count data and identified differentially expressed genes using planned linear contrasts, as described in Mikheyev & Linksvayer [45]. Second, we analysed focal pairwise comparisons to identify groups of genes that were significantly differentially expressed between wasp phenotypes. We compared the following wasps: workers versus gynes, stylopized versus workers and stylopized versus gynes. For this approach, read counts were log₂ transformed and corrected for skew. We detected differential levels of gene expression using a modified Fisher's exact test, which took into account the dispersion and the multiple samples. Finally, raw p-values for each gene were corrected for multiple comparisons with an false discover rate (FDR) cut-off of 0.05. The HTSeq raw read counts for the 11 506 transcripts, the R scripts for processing the data and creating the Venn diagram, data files for the condition-averaged raw reads counts for the 8484 genes and FDR values for each pairwise comparison are available on the Github repository (https://github.com/ruolin/WaspsProject).

For global analyses of gene expression, we used hierarchical clustering (Ward method) in JMP Pro v. 10.0 (SAS, Cary, NC) and principal component analysis in R. To perform Gene Ontology (GO) analyses, we obtained Drosophila melanogaster orthologues with BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) for all P. dominula transcripts, including those that were significantly differentially expressed between treatments (see electronic supplementary material, table S4 for a complete list). We used GO functional annotation clustering in DAVID v. 6 [46,47] with medium stringency; significance was assessed at p < 0.05 and with a Benjamini–Hochberg correction for multiple testing. To identify overrepresented biological functions (enrichment analysis), we compared the annotation composition in our list of differentially expressed genes with that of a population background composed by all the P. dominula transcripts (represented in the RNA-sequencing data) with Drosophila orthologues.

(e) Real-time qRT-PCR validation

We selected two genes to validate RNA-sequencing data from transcripts found differentially expressed in each phenotype: P. dominula defensin precursor (Defensin; GU327374.1) and P. dominula Immune-response protein 30 (IRP30; JN181874.1). These genes were chosen for validation because they had large and robust expression differences between all three experimental groups. RNA from additional P. dominula samples not already used for RNA-sequencing was used for validation. We only had a small number of samples available for this validation exercise (n_W = 4, n_G = 3 and n_S = 5). cDNA was transcribed from individual RNA samples, spiked with an external control sequence (mCherry [48], with SuperScript III First-Strand Synthesis Kit (Invitrogen)). qRT-PCR was performed, using Power SYBR^® Green PCR Master mix (Life Technologies), on the CFX384 Touch™ Real-Time PCR Detection System.

Control genes were identified from the P. dominula genome [24], using best BLAST hits to control genes used in other studies and related taxa: these genes were homologous to other hymenopteran Actin and Elongation-factor 1α (EF1α) sequences. We also used an external control mCherry (Discosoma sp., [48]). Primer sequences for each of the control and focal genes were designed using Primer Quest (Integrated DNA Technologies). Full primer sequences and additional protocol details are found in electronic supplementary material, table S4.

3. Results

(a) Sequencing data

Fourteen libraries of SOLiD RNA-seq generated 904 464 442 single-end reads of 76 bases in length. The majority (72%) of reads (657 626 532) passed quality control steps. Most (11) of the samples preserved nearly 90% reads after filtering. There were four samples retaining approximately 30% reads after filtering (electronic supplementary material, table S3), and these were kept in the analysis because they were evenly distributed across the sample groups. Using BFAST [49] to align SOLiD RNA-seq reads, 502 657 739 (76.4%) reads were uniquely mapped to the reference genome of P. dominula [24].

(b) Gene expression patterns

We based estimates of gene expression only on uniquely mapped reads aligned to the P. dominula gene annotation r1.0 [24]. Overall, we identified 367 differentially expressed transcripts (DETs) (see electronic supplementary material, table S5). Hierarchical clustering analysis of all DETs (figure 2a) revealed that the overall gene expression patterns of stylopized females were more similar to those of gynes than workers, despite the fact that stylopized females were collected on worker phase nests alongside other normal workers. Gyne-shifted gene expression patterns of stylopized females were also apparent from a principal component analysis on the 367 DETs (see electronic supplementary material, figure S1).

• Download figure
• Open in new tab
• Download PowerPoint

Overall, gynes were the most divergent group in terms of brain gene expression. Gynes had a total of 282 DETs (142 were upregulated and 140 were downregulated, p < 0.05, FDR-adjusted) relative to the other two groups. Stylopized wasps were the least divergent in brain gene expression relative to the other two groups (59 DETs: 53 upregulated and six downregulated), which agree with the idea that this phenotype has intermediate features between worker and gyne phenotypes. Finally, workers had 94 DETs relative to the other two wasp phenotypes (39 upregulated and 55 downregulated).

Analysis via focal pairwise comparisons (figure 2b) identified 305 caste-related DETs between workers and gynes (146 upregulated in workers and 159 downregulated), 51 DETs between stylopized and workers (six upregulated in stylopized and 45 downregulated) and 90 DETs between stylopized and gynes (30 upregulated in stylopized and 60 downregulated). There were only three DETs that differed between all three contrasts (figure 2b, centre), including two immune-related genes, Defensin and IRP30. Because of these strong expression differences, these two genes were used for qRT-PCR validation of the RNA-seq results (described below).

Of the 51 DETs that differed between stylopized females and workers, 14 were also related to caste differences. This gene set represents top candidates for parasite manipulation of host caste phenotypic plasticity and is labelled as ‘Parasite manipulation candidate genes’ in figure 2b. These genes include five uncharacterized proteins (electronic supplementary material, table S7), as well as genes putatively related to hormone binding, transcriptional regulation, oxidoreductase activity, cell adhesion and proteolysis (table 1). Several of these functions have been previously associated with caste differential gene expression in Polistes [24,29]. Strikingly, a comparison of read count ratios for the 14 overlapping DETs between workers and gynes (W : G ratio), and workers and stylopized females (W : S ratio) were tightly correlated (Pearson correlation coefficient R² = 0.9831, p < 0.0001; figure 2c).

(c) Gene functions

We performed GO analysis separately on each of the pairwise contrasts in order to identify functional categories of genes associated with caste differences (worker versus gyne) and parasitism (stylopized versus gyne and stylopized versus worker); these DET lists returned 43 clusters, nine clusters and one cluster of GO terms, respectively (electronic supplementary material, table S6). For the caste comparison (worker versus gyne), the top clusters related to largely biosynthetic and metabolic functions (including oxidoreductase and cytochrome P450), stimulus response and post-translational gene silencing. For the stylopization comparisons, stylopized versus gyne clusters primarily consisted of lifespan and metal-binding functions, and the stylopized versus worker cluster was related to immune function (electronic supplementary material, table S6); however, no terms were significantly enriched after correction for multiple testing.

(d) qRT-PCR validation of select RNA-seq results

We selected two genes to validate our RNA-seq results using qRT-PCR. We chose Defensin and IRP30 because both of these genes showed strong expression differences between the three groups (note that they are not among our top candidate genes for host manipulation). We found overall similar gene expression patterns between our RNA-seq results and qRT-PCR validation (using different biological replicates) for both Defensin and IRP30 (figure 3). For Defensin, expression was extremely high in stylopized females, intermediate in gynes and lowest in workers in both RNA-seq and qRT-PCR (figure 3a). For IRP30, expression was lowest in workers and higher in stylopized females (figure 3b), but the pattern in gynes was not consistent across the two types of analysis. The Wilcoxon tests showed significant differences for IRP30 in the qRT-PCR (IRP30: χ² = 6.1091, p = 0.0471) and RNA-seq. For Defensin, the qRT-PCR was not significant (p = 0.20), but a trend in the same direction as the RNA-seq data was apparent.

• Download figure
• Open in new tab
• Download PowerPoint

4. Discussion

In this study, we provide novel evidence supporting the hypothesis that parasites can manipulate host behaviour by inducing changes in gene expression related to innate host phenotypic plasticity. The main finding of this study is that X. vesparum appears to be manipulating its P. dominula wasp host by partially shifting its social caste, and this is associated with changes in the expression of caste-related genes. Three main lines of evidence support this idea. (i) Overall gene expression patterns in stylopized females were shifted towards being gyne-like (figure 2a). (ii) A substantial fraction of the genes that changed in expression because of stylopization were caste-related genes (i.e. also differed in expression between normal gynes and workers, figure 2b). (iii) A subset of genes differentially expressed because of stylopization closely mirrored patterns found in gynes versus workers (figure 2c). These data provide new insights into how parasites can manipulate host brains and behaviour, and suggest that hijacking of plasticity in gene expression could be an effective mechanism by which parasites can shift multiple, coordinated host phenotypes to their advantage.

The most pronounced differences in gene expression (as far as the number of genes) that we uncovered were related to caste differences. There were numerous gene expression differences between workers and gynes, including genes with putative functions in oxidation–reduction processes, carbohydrate and lipid metabolism, and visual stimulus response. Several of these gene functions have been identified in previous studies of caste differential gene expression in Polistes metricus [29,42] and other social insect taxa [51], suggesting that these data reflect typical caste-related expression patterns. We saw fewer differences in gene expression associated with stylopization than related to social caste, suggesting that X. vesparum targets a fairly specific set of transcriptomic and behavioural responses. Some stylopization-related genes also had putative functions associated with immunity, with high expression of both Defensin and IRP30 in stylopized females (figure 3). Upregulation of Defensin in stylopized females agrees with a previous study on X. vesparum effects on larval wasp immune response [52], suggesting that wasps are immunostimulated during the stylopization process.

Although we saw clear evidence of the predicted shift in stylopized female brain gene expression towards being more gyne-like (figure 2), our data do not allow us to conclude whether these gene expression responses are the cause or the consequence of the observed behavioural shifts. Importantly, stylopized wasps were collected before they exhibited gyne-like nest desertion and aggregation behaviour. This suggests that some of the identified genes could have causal roles in some of the subsequent parasitoid-induced shifts in caste-related behaviour. Most promising candidates include genes related to transcriptional regulation (knirps and HMX3-B), oxidation–reduction (hydroxypyruvate reductase) and behavioural response to starvation (Takeout, figure 2c and table 1). However, more studies are needed to understand whether these gene expression changes are a direct cause of parasitoid manipulation, or are downstream from other parasitoid-induced changes in host physiology or hormones.

This study also highlights the potential of using host–parasite systems to elucidate the molecular basis of complex phenotypes. Although hundreds of caste-related gene expression differences have been previously identified in Polistes [53], it has been difficult to pinpoint which genes relate to different aspects of the suite of caste-related behavioural and physiological traits. Because stylopization only appears to target certain modules of gyne-like behaviour (such as nest desertion), rather than the entire phenotype, it can serve as a ‘natural experiment’ [18] to help hone in on key molecular components involved in specific caste endophenotypes.

In summary, we provide novel data on the transcriptomic basis of parasite manipulation of host social behaviour. Our data suggest that X. vesparum shifts the brain expression of genes related to plastic social castes in their host, eliciting an ordinarily quiescent gyne phenotype in the ‘wrong caste’. This is an excellent example of subtle behaviour manipulation, and our results suggest that parasitism can solicit a complex behavioural response in hosts via small but specific shifts in host brain gene expression. We suggest that plastic host behavioural and physiological responses, such as alternate morphs or even personality types [7], may be particularly susceptible to parasite manipulation. By tapping into underlying host phenotypic plasticity, parasites could more efficiently alter multiple, correlated host behavioural and physiological traits via fairly targeted shifts in host regulatory mechanisms. Further research is needed to investigate the generality of such a mechanism across other host–parasite systems, but the plasticity perspective has the potential to provide new insights into our understanding of host–parasite interactions and their coevolution.

Ethics

Wasps in this study were collected and euthanized humanely, according to standard collection protocols.

Data accessibility

Raw data files from the SOLiD RNA-sequencing project described above have been deposited in the NCBI short read archive, BioProject (BioProject: PRJNA312511, Experiments SRX1629425 through SRX1629436). All scripts used to process RNA-Seq data can be found in the following repository (https://github.com/ruolin/WaspsProject).

Authors' contributions

A.C.G.: brain dissection, RNA extraction and submission of samples for sequencing, qRT-PCR validation, data analysis and writing; R.L.: processing of raw sequencing data, preparation of DET dataset and writing; F.M.: sample collection, preparation of samples for RNA-sequencing and writing; L.B.: study conception, sample collection, experimental design and writing; J.K.: study conception, sample collection and writing; C.M.G.: RNA-sequencing/analysis infrastructure and resources; A.L.T.: study conception and experimental design, sample collection and processing, and writing.

Competing interests

We have no competing interests.

Funding

No funding has been received for this article.

Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3726622.