ORIGINAL RESEARCH article
Sec. Chemical Ecology
Variation in the Venom of Parasitic Wasps, Drift, or Selection? Insights From a Multivariate QST Analysis
- 1Université Côte d'Azur, INRA, CNRS, ISA, Sophia Antipolis, France
- 2Université de Rennes 1, UMR-CNRS ECOBIO 6553, Campus de Beaulieu, Rennes, France
- 3Université Lyon 1, Laboratoire de Biométrie et Biologie Évolutive, UMR CNRS 5558, Villeurbanne, France
Differentiation of traits among populations can evolve by drift when gene flow is low relative to drift or selection when there are different local optima in each population (heterogeneous selection), whereas homogeneous selection tends to prevent evolution of such a differentiation. Analyses of geographical variations in venom composition have been done in several taxa such as wasps, spiders, scorpions, cone snails and snakes, but surprisingly never in parasitoid wasps, although their venom should constrain their ability to succeed on locally available hosts. Such a study is now facilitated by the development of an accurate method (quantitative digital analysis) that allows analyzing the quantitative variation of large sets of proteins from several individuals. This method was used here to analyse the venom-based differentiation of four samples of Leptopilina boulardi and five samples of L. heterotoma from populations along a 300 km long south-north gradient in the Rhône-Saône valley (South-East of France). A major result is that the composition of the venom allows to differentiate the populations studied even when separated by few kilometers. We further analyzed these differentiations on the populations (reared under similar conditions to exclude environmental variance) with a QST analysis which compared the variance of a quantitative trait (Q) among the subpopulations (S) to the total variance (T). We also used random forest clustering analyses to detect the venom components the most likely to be adapted locally. The signature of the natural selection was strong for L. heterotoma and L. boulardi. For the latter, the comparison with the differentiation observed at some neutral markers revealed that differentiation was partly due to some local adaptation. The combination of methods used here appears to be a powerful framework for population proteomics and for the study of eco-evolutionary feedbacks between proteomic level and population and ecosystem levels. This is of interest not only for studying field evolution at an intermediate level between the genome and phenotypes, or for understanding the role of evolution in chemical ecology, but also for more applied issues in biological control.
Local adaptation, the fact that populations are more adapted to their own local environment than to another environment, has long been the searched signature of an ongoing natural selection in the field (Futuyma and Moreno, 1988; Kawecki and Ebert, 2004; Poisot et al., 2011). It is considered to be a major cause of the emergence and maintenance of polymorphism (Ravigné et al., 2009). Traits involved in the antagonistic interactions can also be subject to such local adaptations, but with the particularity that the interacting species each represent the environment of the other, which can produce local adaptations induced by coevolution. However, the “geographic mosaic theory of coevolution” predicts that the intensity of such coevolutionary dynamics can vary in space according to the specific composition of local communities (Thompson, 1994, 1997; Gomulkiewicz et al., 2000), with, for example, the presence of some alternative prey or host species.
Local adaptation induced by coevolution is particularly expected between species whose interactions involve resistance-virulence traits, as is the case for host-parasite or host-parasitoid interactions. Parasitoids are insects which develop at the expense of the host, usually leading to its death. For such interactions, Gandon (2002) has shown that the localized coevolution is more likely when the fitness cost to the host and the specificity of the parasite are high. For most parasitoids species, the fitness costs to the host are maximal since most hosts are killed before producing any offspring, and parasitoids virulence has been shown to be specific (Kraaijeveld et al., 1998; Dupas et al., 2009; Rouchet and Vorburger, 2012, 2014). In addition, local adaptation is more likely to occur when organisms are able to choose their environment (Ravigné et al., 2009). Parasitoids may be able to choose host genotypes that maximize their fitness. Indeed, most species of parasitoids were shown to be able to evaluate several traits of a host encountered, such as size, age, parasitism status, or resistance capacity, thanks to numerous types of chemosensory sensilla on their antennae or ovipositor, before deciding to lay an egg or not (van Baaren and Boivin, 1998; Dubuffet et al., 2006; van Baaren et al., 2007). The ability to select its host combined with the specificity of virulence may explain the tendency of similar parasitoid genotypes to parasitize similar host genotypes under certain environmental conditions (Lavandero and Tylianakis, 2013).
If the specificity of virulence and local adaptation exist, they may be supported by the specificity and local adaptation of virulence factors. However, geographic variation of parasitoid virulence factors has only been little studied to date although most adaptive evolutionary constraints are determined at this level of organization (Kawecki and Ebert, 2004; Biron et al., 2006; Wagner, 2012). In many parasitoids, the virulence factors are mainly contained in the venom injected into the host at the time of oviposition, which results in the inhibition of the immune system of the host as well as the regulation of its physiology, or behavior, thus optimizing the development of the parasitoid offspring (Poirié et al., 2009, 2014; Mrinalini et al., 2014; Moreau and Asgari, 2015; Martinson and Werren, 2018;Walker et al., 2018).
Although data for parasitoids are still scarce, geographic variations in venom composition have been investigated in most venomous taxa, e.g., wasps (Perez-Riverol et al., 2017), spiders (Binford, 2001), scorpions (Abdel-Rahman et al., 2009; Rodríguez-ravelo et al., 2013), cone snails (Remigio and Duda, 2008; Abdel-Rahman et al., 2011), and snakes (Wuster et al., 1992; Francischetti et al., 2000; Alape-Giron et al., 2008; Calvete et al., 2011; Holding et al., 2016; Hofmann et al., 2018).
To determine whether a geographical variation is the consequence of local adaptation or drift, tests can be done experimentally or using statistical inferences from field observations. For example, Holding et al. (2016) experimentally evidenced some local adaptations of the venom of rattlesnakes by showing that venom was most effective on squirrels sampled at an altitude similar to that of the rattlesnake from which it was extracted. In addition to these reciprocal-transplant-like experiments, local adaptation can be demonstrated using QST approaches that compare the variance of a quantitative trait (Q) among subpopulations (S) to the total variance (T). They generally help to disentangle the local adaptation from the drift by comparing the differentiation for the traits of interest to that for neutral markers (FST). One of these approaches, the multivariate QST analysis developed by Martin et al. (2008) additionally checks whether traits variances among populations are proportional to their variances within populations, which is expected if the traits are neutral.
In this paper, we have used this multivariate QST approach to assess the effect of drift and selection on the venom of different populations of two parasitoid species of Drosophila, Leptopilina boulardi, and L. heterotoma, along a south-north gradient in which the relative importance of the two species varies. In France, these parasitoids of Drosophila have been the subject of an in-depth follow-up for about 30 years along the Rhône-Saône valley. Drosophila melanogaster and D. simulans, the main Drosophila species in this area, are parasitized by three larval parasitoid species, Asobara tabida, L. boulardi and L. heterotoma. Leptopilina boulardi is generally considered as a specialist of these two host species, while L. heterotoma has been found on eight other Drosophila species as well as flies from the genus Chymomyza and Scaptomyza (Allemand et al., 1999; Fleury et al., 2004; Fleury et al., 2009).
The gradient of temperature along the Rhône valley prevents the establishment of L. boulardi, a Mediterranean species, in its northern part. In the South, L. boulardi outcompetes the other parasitoid species (Allemand et al., 1999; Fleury et al., 2004, 2009). The community is therefore dominated by D. simulans and L. boulardi south of the valley, and D. melanogaster and L. heterotoma in the north (Fleury et al., 2004). However, changes have occurred over the past 30 years in the community with L. boulardi moving northward rapidly (90 km per decade), probably in relation to the increase in temperature associated with global warming (Delava et al., 2014). The presence of this new competitor affects the life history traits of A. tabida (Vayssade et al., 2012; van Baaren et al., 2016; Moiroux et al., 2018) but not L. heterotoma, possibly due to a wider overlap of its ecological niche with that of A. tabida (Vayssade et al., 2012).
Leptopilina parasitoids have been widely studied for various life history traits [fecundity, egg size, egg load at emergence, lifetime, ovigeny index), physiology (lipid content, metabolic rate), mobility (e.g., locomotor activity)] and features associated with parasitism success or competitiveness (Fleury et al., 1995; Vayssade et al., 2012; Vuarin et al., 2012), notably in the Rhône valley. For example, L. heterotoma shows a strong intraspecific genetic variability in the allocation of resources to different life history traits (Vuarin et al., 2012). Genetic differences in locomotor activity were also observed, with parasitoids of the southern populations being more active at the beginning and the end of the photophase (lower temperature), and those of northern populations during the afternoon (higher temperature; Fleury et al., 1995).
The existence of a spatial structuring of populations, a prerequisite for local adaptation, was demonstrated in Iranian populations of L. boulardi using neutral markers (Seyahooei et al., 2011). Moiroux et al. (2013) also showed local adaptation in L. boulardi, with significant differentiation between orchards and forest habitats, possibly related to a different distribution of resources (e.g., host distribution) in these environments. A similar pattern may explain the characteristics of a Congolese line of L. boulardi that is less effective on D. melanogaster than Mediterranean populations but can develop successfully on the tropical species D. yakuba. A genetic analysis using this line and a line from a Mediterranean population suggested that these differences involve two non-linked loci, each necessary for parasitic success on one of the two hosts (Dupas and Carton, 1999). Interestingly, the Congolese line reveals resistance in D. melanogaster and D. yakuba and has highlighted two major resistance genes, one in each species (Hita et al., 1999, 2006). For L. heterotoma, a comparison of the Seattle (USA) and Igé (France) populations revealed a local adaptation to the host D. subobscura, the parasitoid fecundity being higher when parasitizing the sympatric host (Gibert et al., 2010). Overall, the explanations for populations differentiation mainly involved local adaptation, either to host distribution (Moiroux et al., 2013) or to the presence of competitors (Vayssade et al., 2012).
Surprisingly, although frequent local adaptation of traits involved in virulence is expected, the existence of a population structure based on such traits has rarely been tested. For example, Dupas et al. (2003) conducted such analysis worldwide on populations of L. boulardi, showing that the host resistance rate to the Congolese line was broadly similar regardless of their geographic origin, while parasitoids from tropical Africa were less virulent on D. melanogaster compared to those from other areas.
In Leptopilina parasitoids, virulence depends on the venom injected with the egg at oviposition (Poirié et al., 2009, 2014; Moreau and Asgari, 2015; Walker et al., 2018). Although never conclusively demonstrated, there is ample congruent evidence to suggest that differences in venom composition of parasitoids may explain their difference in host range (e.g., Lee et al., 2009). Parasitoid venoms are complex fluids containing many proteins that have been characterized in several species, including L. boulardi, and L. heterotoma (Colinet et al., 2013a; Goecks et al., 2013). Colinet et al. (2013a) investigated the interspecific and intraspecific variation of venom composition by comparing two L. boulardi lines having different virulence properties and a strain of L. heterotoma. L. boulardi and L. heterotoma had no major venom proteins in common while L. boulardi lines shared only about 50% of their most abundant proteins. Finally, inter-individual variations of venom composition were also recently evidenced in strains and populations of species of the Leptopilina and Psyttalia genera (Colinet et al., 2013a,b). This variation in the venom composition may also be at the origin of the rapid evolution of parasitoid virulence observed under experimental evolution (Dion et al., 2011; Rouchet and Vorburger, 2014). Indeed, differential selection of the expression levels of genes encoding venom proteins was observed during the evolution of parasitoid virulence in this last experimental evolution (Dennis et al., 2017). Nevertheless, neutral molecular variation is common, and some of the variations in the venom composition are therefore most likely neutral and do not affect virulence. This raises the question of whether parasitoid populations could be differentiated according to the composition of the venom and which of these differences would be caused by drift or local adaptation.
Here, we have investigated the population structure of L. boulardi and L. heterotoma in the Rhône-Saône valley according to the individual composition of the venom. For L. boulardi, we have studied four populations, two in the south of the Rhône-Saône valley that this species invaded more than 40 years ago (Delava et al., 2014) and two at the northern end of this valley (close to Lyon) that was invaded by this species between 1993 and 2003 (Delava et al., 2014). For L. heterotoma, five populations were studied, one in the south, outside the Rhône-Saône valley (200 km east) where L. boulardi predominates, two in the middle of the Rhône-Saône valley where it appeared between 1993 and 2003, and two in the Saône valley, north of Lyon, where its presence was detected between 2003 and 2011 (Delava et al., 2014). The analysis of the venom composition of individual wasps was performed using a recently developed method based on 1D electrophoresis and further digital and statistical analysis of the intensity of the protein bands (Mathé-Hubert et al., 2015).
We report for the first time a structuration of parasitoid populations based on venom composition, a trait often directly related to the outcome of host-parasitoid interactions. This differentiation among populations in venom composition is partly induced by drift but may also be due to some heterogeneity in the selection pressures experienced by different populations. We used the multivariate QST approach developed by Martin et al. (2008) to detect whether there is selection for the same optimal relative amounts of venom components in all populations (homogeneous selection) or selection for different local optima (heterogeneous selection, or local adaptation), possibly in relation with some variation in the host physiology or resistance. For L. boulardi, data showed that part of the venom-based population structure is induced by a certain heterogeneity in selection within different populations, i.e., some local adaptation. For L. heterotoma, we could only perform part of the analysis, which evidenced a high evolutionary potential of the venom composition. These results are compared with those previously obtained for various parasitoid traits in the Rhône-Saône valley and are discussed in the context of coevolutionary theories.
Materials and Methods
Sampling and Analyzed Individuals
The four L. boulardi populations were sampled in orchards in September 2010. Sampling locations were Ste Foy-Lès-Lyon and St Laurent d'Agny for northern populations (25 and 29 females sampled, respectively) and Éyguières and Avignon for southern populations (20 and 11 females sampled, respectively). Southern and northern L. boulardi populations are distant from 200 km (Figure 1 and Table S1). Each field-collected individual was isolated and allowed to parasitize larvae of a D. melanogaster strain sampled in Ste Foy-Lès-Lyon. The offspring of these individuals were then stored separately at −80°C. The venom of one daughter was analyzed per field-sampled female.
Figure 1. Sampling locations of L. heterotoma (green) and L. boulardi (red). Northward progression of L. boulardi was evidenced by Delava et al. (2014): L. boulardi invaded in locations ➀ between 1993 and 2003 and in locations ➁ between 2003 and 2011.
The five populations of L. heterotoma were sampled in orchards in 2008 (n = 15–20) and then maintained on D. melanogaster (Ste Foy-Lès-Lyon) until the beginning of September 2013, when individuals were stored at −80°C. As discussed later, we cannot exclude that some of the evolution we detected occurred in the laboratory between 2008 and 2013. For this reason, we will refer to these populations as “laboratory-maintained populations.” The five laboratory-maintained populations originated from field samples from Uchizy and Montbellet for northern populations, Sonnay and Épinouze for middle populations, and Vence as a southern population but located outside the Rhône-Saône valley. Northern and middle populations are distant from 130 km whereas Vence is located 350 km away from northern populations (Figure 1). We analyzed 20 individuals for each laboratory maintained population.
Venom Characterization: Sample Preparation and Venom Profile Digital Analysis
The venom apparatus were dissected individually, and venom extracted from the venom reservoir treated as described in Mathé-Hubert et al. (2015), and loaded onto 1D SDS-PAGE gels (Any-kD Mini-PROTEAN® TGXTM, Bio-Rad) for L. heterotoma (eight gels) and L. boulardi (seven gels). For both species, we ensured that all gels contained individuals from all studied populations. Briefly, gels were silver stained after migration (Morrissey, 1981) and photographed several times during the protein revelation step (digital camera EOS-5D-MkII, Canon, Japan). Taking several pictures during protein revelation allowed to account for most of the heterogeneity in protein quantities by ensuring that, for all lanes, at least one picture was available on which all shared bands were revealed without being too saturated (see examples of pictures in Figure 2). The high-resolution pictures (5,626 × 3,745 pixels; 16 bits; TIFF file) were then analyzed semi-automatically with a quantitative digital analysis based on lanes transformation into intensity profiles by Phoretix 1D (TotalLab, UK), which were then analyzed by R functions. These functions use the 2nd derivative of these profiles to semi-automatically identify bands that are common to a large number of individuals (hereafter “reference bands”). The intensity of these reference bands is then quantified in each individual lane and normalized to account for the remaining heterogeneity in protein amounts (see details in Mathé-Hubert et al., 2015). The analysis resulted in the choice of 29 and 32 “reference bands” of identified molecular weight for L. heterotoma and L. boulardi, respectively. Then, to remove the remaining variation in the total lane intensity and to measure the intensity of these reference bands in each lane, we used the following combination of parameters for L. heterotoma (Lh) and L. boulardi (Lb), [(i) “height” for Lh and “volume” for Lb (maximal intensity between borders of reference bands or total intensity between borders, respectively), (ii) background removed in Phoretix-1D with a “rolling ball” of 10,000 pixels of radius for both species, (iii) cyclicloess for Lh and quantiles normalization for Lb]. These two combinations of parameters were selected by comparing the signal to noise ratio and the multicollinearity of the dataset produced by each combination of parameters (details in Mathé-Hubert et al., 2015). These normalized intensities of the reference bands are the variables that describe the venom composition.
Figure 2. Example of SDS-PAGE gels for venom analysis. Example of SDS-PAGE gels used to characterize the variation in the venom composition of L. boulardi (A) and L. heterotoma (B). Each lane contains the venom of a single wasp whose population of origin is indicated at the top of the lane (Éy, Éyguières; Av, Avignon; SFL, Ste Foy-Lès-Lyon; SLA, St Laurent d'Agny; Mo, Montbellet; Ve, Vence; So, Sonnay; Ép, Épinouze; Uch, Uchizy). The reference bands that significantly discriminate between populations after a Bonferroni correction are indicated.
Lastly, the intensities of some reference bands can be highly correlated, either because they partly overlap (see Figure 2) or because of linkage disequilibrium (for more details, see Mathé-Hubert et al., 2015). Such a linkage disequilibrium between genes coding for a higher or a lower expression of different venom proteins may indeed induce some correlation between the bands containing the two proteins. Since the linkage disequilibrium between traits should be low for the multivariate QST analysis (Martin et al., 2008), we merged reference bands that were too highly correlated. These set of correlated reference bands were detected using a UPGMA dendrogram in which the distance between pairs of reference bands was measured as 1 − |r| where r is the Pearson correlation coefficient (Figure S2). Then, we merged the bands whose correlation was higher than 0.6 in a composite reference band estimated as the mean of the merged bands. For L. boulardi, we built four composite reference bands from reference bands (16, 17), (19, 20), (22, 28), and (30, 31, 32) (Figure S2A). For L. heterotoma, we built six composite reference bands from reference bands (1–3), (6, 7), (8, 9), (15, 21), (19, 20), and (28, 29) (Figure S2B).
Statistical Analyses of the Population Structuration
We performed all statistical analyses in R 3.4.4. The workflow of these analyses is described in Figure S1.
Power of the Venom to Discriminate Populations
For each pair of populations, we used a Random Forest (RF) clustering algorithm to predict the population of origin of individuals based on the composition of their venom (the centered and scaled intensities of the reference bands). We conducted the RFs using the package party to account for the multicollinearity and high number of predictors compared to the low number of individuals per population (Hothorn et al., 2006; Strobl et al., 2008). For each RF, we have estimated the discriminating power of the venom composition as the proportion of accurately classified individuals. The significance of this discriminatory power was estimated using 10,000 permutations. To avoid overfitting, this was done using predictions based on the out-of-the-bag data. For all RFs fitted to L. boulardi, we used the weight argument to ensure that all populations have the same weight albeit the headcounts were not balanced.
Effect of Selection on the Variation in Venom Composition Within and Among Populations
We used a QST analysis to disentangle between drift, heterogeneous selection (local adaptation) and homogeneous selection (Figures S1A–C). The QST compares the variance of a quantitative trait (Q) among the subpopulations (S) to the total variance (T).
As we had no information on the genetic relatedness between the analyzed individuals, this QST analysis was not performed on the additive genetic variance (strict-sense QST) but on the variation between individuals raised in a common environment i.e., the total genetic variance (broad-sense QST). This broad-sense QST may be biased downward which is further described in the discussion. The multivariate QST analysis developed by Martin et al. (2008) relies on two complementary tests, the first being more powerful when different traits (here venom components) are subjected to different types of selection regime (i.e., homogeneous, heterogeneous, or neutral), the second when they all undergo the same selection pressures (either homogeneous or heterogeneous). For more details see Figure S1C. The R code used to perform these tests was developed by Martin et al. (2008). To ensure its long-lasting accessibility, it is also included in the Supplementary Material S1.
Test 1.H0: Covariances amon g (D) and within (G) populations are proportional (D = ρ × G)
For both species, we first summarized the total variation in the venom composition using a principal components analysis (PCA; Dray and Dufour, 2007) in order to visualize how this variation is partitioned within and among populations. Then, we plotted the first two principal components (PC) of the PCA using the R function s.class (package ade4) that groups individuals of the same population. This allows to determine whether, as expected under neutrality, the dimension with the greatest variation within each population is also the one with the greatest variation between populations. To improve this visual comparison of the variation within vs. among populations, we projected on these plots the first PC of a PCA performed on the variation within populations and the first PC of a PCA performed on the variation among populations. These two additional PCA were performed using the functions wca and bca (package ade4), respectively. If the variation among populations is only induced by drift, then these axes of variation within and among populations should be mostly collinear. We measured this collinearity with the Pearson correlation coefficient. A function performing this analysis was added to the Supplementary Material S1 (see also Figure S1A).
This expectation of proportionality between the variations within and among populations is the null hypothesis of the first test of the multivariate QST analysis developed by Martin et al. (2008). Indeed, under the null hypothesis that the components of venom are neutral, all the variations observed among populations appeared because the isolated populations drifted independently of each other. In each population, for each trait (here venom components), the amount of drift is proportional to the amount of variance within the population. This gives the statistical expectation that the average matrix of covariances within populations (G) is proportional to the matrix of covariances among populations (D): D = ρ × G, where the two matrices D and G are estimated by a MANOVA and the coefficient ρ of proportionality is estimated by a function developed by Martin et al. (2008). As the sample size per population is small, the p-values were computed using the correction based on the Bartlett adjustment of likelihood ratios (Martin et al., 2008). In this analysis, the confidence interval of the proportionality coefficient ρ is estimated by maximum-likelihood (for more details, see Martin et al., 2008). This confidence interval is used for the second test.
Test 2.H0 The proportionality coefficient
The second test relies on the comparison between the average differentiation observed for the traits of interest (measured by the proportionality coefficient ρ) and the differentiation observed for some neutral markers (measured by the FST). Information on such neutral differentiation was only available for three of the four populations of L. boulardi studied (Avignon, Ste Foy-Lès-Lyon, and St Laurent d'Agny) and it was not available for L. heterotoma. This test was therefore performed only on these three L. boulardi populations. Under the null hypothesis that, on average, the traits used to compute ρ are neutral, then ρ = . This is typically tested by checking for an overlap between the confidence intervals of ρ and that of (Martin et al., 2008). The only information we had on neutral differentiation was a point estimate of the L. boulardi FST (= 0.084). However, we did not need the confidence interval of the FST because it turned out that the point estimate of itself fall within the confidence interval of ρ. This estimation of the FST is based on the RAD sequencing of L. boulardi individuals from the Rhône–Saône valley which led to the genotyping of about 16 individuals per population (the three studied populations) on 484 loci (Delava, in preparation).
Identification of the Venom Components Possibly Under Heterogeneous Selection
The identification of venom components possibly under heterogeneous selection was done in two steps. First, we detected the reference bands most likely to be the ones under heterogeneous selection, and then we identified the components of the venom they contain.
For the first step, we used a clustering RF algorithm to predict, for both L. heterotoma and L. boulardi, the population of origin of the individuals based on the composition of their venom (the centered and scaled intensities of the reference bands). We conducted the RFs and computed the conditional importance of each explanatory variable (i.e., reference bands) using the package party (Hothorn et al., 2006; Strobl et al., 2008). We then tested the significance of each explanatory variable by comparing its conditional importance, its estimated null distribution using 10,000 permutations (Hapfelmeier and Ulm, 2013). The p-values resulting from this analysis were then Bonferroni corrected and we selected for further investigation the reference bands that were still significant.
We tentatively identified the proteins present in these reference bands by matching these bands with the bands on the 1D electrophoresis gels used for L. boulardi and L. heterotoma venom proteomics (Colinet et al., 2013a). Since one band can contain different proteins and a given protein can be found in different bands, the number of peptides matches from the mass spectrometry performed by Colinet et al. (2013a) was used to classify the proteins in the bands as abundant or not using a threshold of 10 peptides matches in the proteomic analysis. In these bands, low abundance proteins are unlikely to drive the observed population structure.
Power of the Venom to Discriminate Between the Populations
In the first step of the analysis, we checked whether the composition of the venom differed among populations, and to what extent these differences were strong enough to predict the population of origin of the wasps based on their venom composition. To the best of our knowledge, this question had never been investigated in a parasitoid wasp. Our field sampling design was particularly well-suited for such analysis since the distance between populations varied wildly (from 14 to 230 km for L. boulardi and from 3.3 to 350 km for L. heterotoma).
The RFs on each L. boulardi population pair indicated that, on average, venom accurately classified 74–89% of individuals (Table 1). These classifications were significant for three combinations: two of the four South-North comparisons (Avignon vs. Ste Foy-Lès-Lyon and Éyguières vs. St Laurent d'Agny) and one of the two comparisons of populations of the same area (Ste Foy-Lès-Lyon vs. St Laurent d'Agny, two populations only 14 km apart). Furthermore, the classification was only a little more accurate when comparing distant populations than when comparing nearby populations (average accuracies: 81 and 74.1% respectively; Table 1). The RFs performed on each pair of L. heterotoma laboratory-maintained populations allowed to classify 85 to 97% of individuals by their venom (Table 2). These classifications were significant in all comparisons. The highest accuracy was achieved when comparing Vence to another population (average accuracies: 95.6%; Table 2), but, excluding Vence, there was almost no difference between comparisons involving distant populations or nearby (average accuracies: 90.6 and 87.5% respectively; Table 2).
Effect of Selection on the Variation in Venom Composition Within and Among Populations
To detect whether selection affected variation in venom composition among and within populations, we used multivariate QST analysis, which performed two tests. The first investigated whether all traits were under the same selection regime (homogeneous, heterogeneous or neutral) by checking whether the (co)variances among populations were proportional to the (co)variances within populations. If such proportionality exists, then the dimensions that have the most variance among and within populations should be collinear. We provided a visual assessment of this part of the analysis based on PCAs. The second test investigated whether, on average, traits were under homogeneous or heterogeneous selection by testing whether the coefficient of proportionality (ρ) estimated in the first test was different from .
For L. boulardi, 34.3% of the total variation in venom composition was summarized by the first two axes of the PCA. The axis summarizing the variation among populations discriminated mainly the populations of the North (Ste Foy-Lès-Lyon and St Laurent d'Agny) of those of the South (Avignon and Éyguières), although it also partly discriminated between these two southern populations. This axis was mainly collinear with the first axis of the global PCA. On the opposite, the axis summarizing the variation within populations was mainly collinear with the second axis of the global PCA (Figure 3A). This reveals that instead of the expected collinearity under neutrality, the two axes were mainly orthogonal (Pearson correlation of 0.09), thus showing that the reference bands did not all undergo the same type of selection (more details on Figure S1C).
Figure 3. PCAs on the venom composition. The variation in the venom composition of L. boulardi (A) and L. heterotoma (B) is summarized using a PCA. Each dot corresponds to a single wasp that is positioned according to its venom composition. Then, the individuals of the same population are grouped and the ellipse summarizes their dispersion by indicating the standard deviation within the population. The percentage of variance explained is given at the end of each axis. The solid blue line and the dotted green line are the projections of the axes maximizing the within and the among populations variance, respectively. Under neutrality, they are expected to be collinear.
A similar result was obtained for L. heterotoma. The first two axes of the PCA explained 39.9% of the total variation in venom composition, and the within and among populations axes were also predominantly orthogonal, with a Pearson correlation of 0.16. On the among populations axis, the population of Vence was completely separated from other populations and that of Montbellet was partially separated (Figure 3B).
In agreement with this visual analysis of the venom-based structure of populations, the first test of the QST analysis yielded a highly significant result (for both L. boulardi and L. heterotoma p < 0.001; H0: D = ρ × G where G and D are the within and among populations covariance matrices, respectively and ρ is the proportionality coefficient; see section Materials and Methods). This analysis estimated the proportionality coefficient and its 95% confidence interval (L. boulardi: ρ = 0.243; ρϵ[0.176; 0.395], which corresponds to a QST of 0.108 [0.081; 0.165]; L. heterotoma: ρ = 0.738; ρϵ[0.567;1.06], which corresponds to a QST of 0.270 [0.221; 0.346]). The coefficient of proportionality of L. heterotoma was three times higher than that obtained for L. boulardi and the confidence intervals did not overlap, meaning that this difference was significant. For L. boulardi we also compared the estimated coefficient of proportionality (ρ) with the FST to perform the 2nd test which allowed to check whether there was some homogeneous or heterogeneous selection (more details on Figure S1C). This test was not significant since the estimate of FST = 0.084 gave an estimate of , which falls within the confidence interval of ρ.
Identification of the Venom Components Possibly Under Heterogeneous Selection
Reasoning that the locally adapted reference bands are the ones that discriminate the most between populations, we used RF to predict the populations of origin of individuals according to the intensity of reference bands, and we selected for further investigation those which remained significant after a Bonferroni correction. For L. boulardi, this led to identification of nine references bands (Figure 4) which, in agreement with the PCA (Figures 2, 3A), mainly discriminated between populations of the north and populations of the south, except for the B12 bands which mainly discriminated Ste Foy-Lès-Lyon from other populations. Also, as observed on the PCA, Avignon was more strongly discriminated than Éyguières (bands B26 and B29), although this tended to be the opposite for bands B2 and B6.
Figure 4. Venom protein bands that discriminate between L. boulardi populations. Each plot shows the distribution intensities of a reference band (centered and scaled). B30_31_32 is a composite reference band. The Bonferroni corrected p-values are shown at the top with the ID of the reference bands. Reference bands are sorted from low to high molecular weight. These bands can be observed in Figure 2A.
For L. heterotoma, this allowed identification of 13 references bands that, in agreement with the PCA (Figure 3B), distinguished mainly between Vence and other laboratory-maintained populations (B4, B6_7, B8_9, B16, B23, and B28_29). However, there were some other differences (Figure 5). Montbellet was separated from other laboratory-maintained populations for bands B5 and B10. The band B12 tended to isolate Épinouze, and the bands B8_9, B14, and B19_20, Uchizy. The band B11 varied along a south-north gradient, with the lowest intensities at Vence (South), then regrouping Sonnay and Épinouze (middle of the sampled area) and finally Uchizy and Montbellet (North). In contrast, the bands B15_21 isolated Montbellet (North) and Uchizy (middle) from the other three laboratory-maintained populations. Most of these populations had different average intensity for band B16.
Figure 5. Venom protein bands that discriminate between L. heterotoma populations. Each plot shows the distribution of the centered and scaled intensities of a reference band which discriminate L. heterotoma populations. The areas under the curves are the same for all populations. Composite reference bands are indicated with an underscore in the ID. Bonferroni corrected p-values are shown next to the ID of the reference bands. Reference bands are sorted from low to high molecular weight. These bands can be observed in Figure 2B.
For both species, a putative function could be assigned to only a few of these reference bands by comparing the 1D protein electrophoreses of this study to that used by Colinet et al. (2013a) for the venom proteomics (Tables 3, 4). For L. boulardi, the band B2 contained two proteins without any predicted function, the band B24 contained a serine protease inhibitor of the serpin superfamily (LbSPNm) and bands B27 and B29 mainly contained RhoGAP domain-containing proteins (LbGAP1, 2, 3, and 5). For L. heterotoma, bands B4 and B20 contained proteins of no predicted function. Bands B16 and B23 contained RhoGAP domain-containing proteins (LhGAP3 and LhGAP2, respectively), an additional protein in the band B23 being a serine protease inhibitor of the serpin superfamily (LhSPN). Bands B28 and B29 contained an aspartylglucosaminidase (AGA) and a protein of unknown function.
Table 3. Correspondence between L. boulardi bands that discriminate between populations and their putative protein content.
Table 4. Correspondence between L. heterotoma bands that discriminate between the laboratory maintained populations and their putative protein content.
The venom is a crucial component of the success of most endoparasitoids and its composition varies among individuals (Colinet et al., 2013b). However, although population biology studies of the geographical variation in venom composition was conducted on several venomous taxa (Wuster et al., 1992; Francischetti et al., 2000; Binford, 2001; Alape-Giron et al., 2008; Remigio and Duda, 2008; Abdel-Rahman et al., 2009, 2011; Calvete et al., 2011; Rodríguez-ravelo et al., 2013; Holding et al., 2016; Perez-Riverol et al., 2017; Hofmann et al., 2018), it was not the case for parasitoids. This is surprising since, in parasitoid wasps, venom should constrain the capacity to succeed on locally available hosts. Performing such a study is now facilitated by the development of a rapid and accurate method to analyse the quantitative variation of large sets of proteins from several individual samples (Mathé-Hubert et al., 2015). This method was used here to analyse the venom-based differentiation of four L. boulardi and five L. heterotoma samples from populations along a 300 km-long south—north gradient in the Rhône-Saône valley (South-East of France).
A major result of these analyses was that the venom composition, characterized by the relative proportions of the different constitutive proteins, allowed a significant discrimination of populations, even if only a few km apart. Such a differentiation among populations can evolve by drift when gene flow is low relative to drift or by selection when there are different local optima in each population (heterogeneous selection), whereas homogeneous selection tends to prevent evolution of differentiation among populations.
We further analyzed these differentiations with a QST analysis that partitions the variance of a quantitative trait (Q) among subpopulations (S) relative to the total population (T), using the multivariate QST analysis developed by Martin et al. (2008). In the strict sense, QST analyses should be based on additive genetic variance, but in practice, they also use QST in the broad-sense that applies to the total genetic variance (additive + epistasis + dominance; Bouétard et al., 2014; Porth et al., 2015), as well as the PST that considers the phenotypic (P) variance, thus also including the non-genetic environmental variance (Brommer, 2011).
In this study, we could exclude the existence of environmental variance since individuals were raised in a common environment, meaning that the phenotypic variance is equal to the total genetic variance. We therefore used broad-sense QST analyses which assume the lack of non-additive genetic variances, an approximation that can bias the QST estimates. This bias is generally considered low because, in most cases, QST is computed for multigenic traits (Leinonen et al., 2013). However, the amount of a given venom component may depend on a low number of loci (although see Albert et al., 2014), so this bias could be more substantial. It generally induces an underestimation of the QST, both when it is due to neglected epistasis variance (Whitlock, 1999) and when due to neglected dominance variance (Goudet and Büchi, 2006; Leinonen et al., 2013) although this has been more controversial (Goudet and Martin, 2007). For a more detailed discussion on this bias, see Leinonen et al. (2013). Downwards biased QST estimates are too conservative for the detection of heterogeneous selection and too sensitive when detecting homogeneous selection, which is reassuring since, in such an analysis, heterogeneous selections are generally more interesting than homogeneous selections.
Variation of the Venom Composition Among Populations: Drift or/and Selection?
The venom-based differentiation between distant populations (South vs. North) was only slightly higher than that between nearby populations. This is quite surprising since some populations were very close to each other. For example, venom significantly discriminated Ste Foy-Lès-Lyon and St Laurent d'Agny, two northern populations only 14 km apart. This differentiation can originate from drift, heterogeneous selection, or a combination of both.
A strong drift could come from the northward extension of distribution of L. boulardi. Indeed, its repartition area has progressed 170 km northwards from 1993 to 2011 (Delava et al., 2014). The wavefront of expanding populations is known to drift significantly because the individuals that drive this expansion are sampled from previously sampled individuals (Travis et al., 2007; Hallatschek and Nelson, 2010). When the wavefront is expanding fast enough, the drift is expected to create a gradient of differentiation from the source to the wavefront (e.g., Moreau et al., 2011). However, only a slight North-South differentiation was observed for L. boulardi, and of the same order as the differentiation between nearby populations. Another explanation would be a metapopulation dynamic with frequent extinctions and recolonizations, due to a high migration rate. Such metapopulation dynamic is expected to create differentiation, even between close populations, since colonization events are associated with founder effects as already demonstrated (Weisser, 2000; Rauch and Weisser, 2007); Nyabuga et al., 2011.
Alternatively, the similar magnitude of differentiation between distant populations (North vs. South) and nearby populations could also be explained if: (i) the effect of drift on differentiations is low, either due to a low drift or a high migration rate and (ii) the heterogeneous selection is strong enough to counter the effect of drift and migration, thus creating some local adaptation. In keeping with this, the QST analysis detected some selection. Indeed, there was a highly significant lack of proportionality between variations within and among populations. Such a proportionality is expected in the absence of natural selection since if the variation among populations comes only from the drift occurring separately within each, the magnitude of drift should be proportional to the magnitude of the variation within populations. This lack of proportionality thus reveals the variation in the selection pressures experienced by the different venom components.
Such variation among traits (here venom components) can have different origins. While the observed variation is most likely neutral for some traits, the non-neutral variation of others may fall anywhere in a continuum described by the following three categories. Firstly, most selection pressures could be homogeneous among populations, which would keep all populations around the same optimal values and give a proportionality coefficient (ρ) lower than . Secondly, most selection pressures could be heterogeneous among populations, which would push them toward different optimal values . Thirdly, there could be a combination of both homogeneous and heterogeneous selection pressures, with no detectable differences between the levels of the two types of selection regimes . This is the situation identified for the venom of L. boulardi. While the presence of homogeneous selection was highly expected for a complex trait related to fitness such as parasitoid wasp venom, that of heterogeneous selection was not, and its level was therefore unknown. The absence of a significant difference between ρ and suggests that both are of comparable strength.
A much higher differentiation of venom was observed for L. heterotoma compared to L. boulardi. This high structuration among populations could have resulted from a low gene flow among populations. However, previous studies in the South-East of France evidenced a low structuration based on neutral markers, indicating thus a consistent gene flow (overall FST = 0.036; Ris, 2003). This leaves two non-mutually exclusive hypotheses for explaining the strength of this venom-based differentiation. Differentiation might be adaptive, resulting from some heterogeneity among populations in the optimal venom composition. Alternatively, all or part of the differentiation may have evolved by drift during the 5 years of laboratory rearing of the populations between field sampling and the venom analysis (see for example Simoes et al., 2010).
As with L. boulardi, there was a very significant lack of proportionality between variations within and among populations that could have been generated by heterogeneous or homogeneous selection or a combination of both. The second test would have allowed discriminating between these hypotheses but was not doable for L. heterotoma. An additional level of uncertainty comes from where the detected selection may have occurred. Indeed, the conditions of rearing in laboratory could have imposed certain modifications in the optimal composition of the venom to which the five laboratory-maintained populations could have adapted by homogeneous selection (van Lenteren, 2003).
In summary, for L. heterotoma, strong differentiation and lack of proportionality may reflect two non-mutually exclusive scenarios. Heterogeneous selection in the field alone could explain both. In addition, the combination of homogeneous laboratory and/or field selection and laboratory drift may explain the lack of proportionality and strong differentiation, respectively. Importantly, whatever the reasons for the observed differentiation and non-proportionality, results highlight the high amount of genetic variation in the venom composition of L. heterotoma.
Local Adaptation of the Venom Composition?
For L. heterotoma, we could not test for local adaptation of the venom composition. However, several ecological features may have created it. Firstly, both spatial and seasonal variations in the relative abundances of Drosophila species have been described in the Rhône valley (Allemand et al., 1999; Ris, 2003). L. heterotoma is a generalist species whose venom may therefore contain either “generalist” virulence factors or combination of factors, each of which being effective for a subset of the host range. In the latter case, the local availability of the hosts could strongly influence the local composition of the venom. This variation in host availability is also likely affected by the competition with L. boulardi, both a specialist and a powerful competitor, which was absent from northern populations in 2003 and present in 2011. Its arrival may have induced changes in the exploitation pattern of Drosophila species by L. heterotoma, resulting in different selection pressures on venom components.
Interestingly, a strong latitudinal genetic differentiation of L. heterotoma was also previously reported, based on fitness traits (Fleury et al., 1995, 2004). However, this could be explained either by differences in abiotic conditions (e.g., temperature) or by a variation in the presence of L. boulardi as a competitor, or both. Another piece of information suggesting some heterogeneous selection could be the putatively adaptive North-South genetic differentiation for the time of preferred locomotor activity, with an ongoing shift between Lyon and Antibes that has been highlighted by Fleury et al. (1995). Indeed, Vence is not far from Antibes, and the middle populations of the Rhône Valley are not far from Lyon. Such a temperature-related difference in the period of activity could have other consequences on both the host and parasitoid physiology, which may affect the optimal venom composition.
For L. boulardi, the QST analysis detected some local adaptation of the composition of the venom, which could be induced by certain biotic or abiotic variations. Local adaptation could indeed be an optimization of something other than virulence if the expression of venom proteins is associated with other physiological traits. For example, if a venom protein is encoded by a gene (i) that is also expressed elsewhere than in the venom gland and (ii) that shows variation in cis-regulation of expression, then a local adaptation could be detected on the venom composition even though the selection takes place due to the protein effect in tissues other than the venom gland. However, such indirect selection on the composition of the venom should be rather rare because the changes it induces are likely not neutral if these venom proteins are costly to produce (Morgenstern and King, 2013).
Alternatively, this might reflect some variability in host availability. Indeed, in the studied area, L. boulardi parasitizes both D. melanogaster, predominant in the North, and D. simulans, predominant in the South (Fleury et al., 2004), and the optimal venom composition differ for success on both species (Cavigliasso et al., in preparation). The North could therefore be a hotspot of co-evolution with D. melanogaster and a coldspot with D. simulans, and vice versa for the South. However, this North-South variation in host availability cannot explain why venom composition differs between close populations. This variation between neighboring populations is the main surprise of this study and deserves to be analyzed with larger samples. As mentioned in the extinction/recolonization scenario, this differentiation could come from drift. If alternatively, it is a local adaptation of virulence, it may be an adaptation to an intra-specific variation in D. melanogaster, the main host of L. boulardi in this region (Fleury et al., 2004).
Such a “very local” adaptation may seem implausible given the low FST between these populations. However, local adaptation induced by coevolution differs from classical local adaptation. For example, since parasitoids and hosts interact antagonistically the local adaptation of one species in a given population is also the local maladaptation of the other species. The interactions between local adaptation and coevolution were studied by Gandon (2002) and Gandon and Nuismer (2009). They reveal that when genetic drift is substantial, higher migration rates actually promote local adaptation, as migrants increase the evolutionary potential (genetic diversity) of the population which is needed to overcome the effect of drift and win the arms race (Gandon, 2002). Also, although drift generally decreases local adaptation, it may sometimes have the opposite effect in the case of local adaptation of traits involved in an antagonistic co-evolutionary dynamic (Gandon and Nuismer, 2009). Indeed, drift occurring independently in each population of two species of a co-evolving “couple” may induce a certain spatial heterogeneity to which both species can adapt locally, which reinforces the spatial heterogeneity initiated by drift. Such simultaneous substantial amounts of drift within each population combined with migration among populations are precisely the predicted situation under the extinction/recolonization scenario we previously considered. Such scenario is not unlikely if we consider (i) the short generation time of these species, (ii) their high fertility, (iii) their ability to fly and disperse over long distances with the wind, (iv) the known seasonal variations in populations size, and finally, (v) the possibility of local extinctions or large fluctuations in the sizes of populations of hosts and parasitoids driven by the parasitic interaction, which would create a strong drift. These conditions are generally likely to create a metapopulation dynamic (Fronhofer et al., 2012).
If we combine this view of host-parasitoid metapopulation dynamic with the notion of informed dispersal (Clobert et al., 2009), then local adaptation between close populations for a trait as important as virulence may be expectable. As parasitoids can accurately estimate the suitability of a host before oviposition (van Baaren and Boivin, 1998; Dubuffet et al., 2006; van Baaren et al., 2007), their decision to settle or continue their migration, when emigrating, may depend on the prevalence of hosts that fit their needs. Such non-random dispersal strongly increases the possibility of local adaptation (Ravigné et al., 2009).
Identifying the Venom Components Likely Under Heterogeneous Selection
To identify the traits under heterogeneous selection, we used the RF clustering algorithm to detect reference bands that discriminate significantly among populations. Indeed, heterogeneously selected venom components provide the greatest discriminating power among populations (high variance among populations relative to the variance within populations). We identified venom components (reference bands) that significantly discriminate between populations. However, even reference bands for which variation is neutral are expected to sometimes discriminate populations. We thus applied a Bonferroni correction to keep only the bands with the highest discriminating power. This selected 9 and 13 reference bands for L. boulardi and L. heterotoma, respectively. These two sets of reference bands should not be interpreted in the same way since we did not formally identify heterogeneous selection for L. heterotoma, which we did for L. boulardi.
The content of these reference bands was tentatively identified by manually matching them with the 1D electrophoresis gels analyzed by proteomics by Colinet et al. (2013a). This proteomics study characterized the venom content of L. boulardi and L. heterotoma lines, which are probably not representative of the whole diversity observed in the field. Nevertheless, there should be an overall correspondence between the content of the reference bands analyzed in this study and that of the bands analyzed by Colinet et al. (2013a). This identified 16 proteins as potentially under heterogeneous selection, six of which have no identified putative function. This highlights that our prior-less method can be used to discover novel active proteins. Indeed, venoms have strong and diverse biological activities but they contain a great diversity of proteins that might not all be worth studying. Our approach allows selecting without a priori proteins likely to contribute to parasitism success on the host. The other identified proteins are members of the RhoGAP and serpin families found in reference bands of both Leptopilina species, as well as an aspartylglucosaminidase (AGA) found in two L. heterotoma reference bands. Among the proteins identified in the discriminating reference bands are LbGAP, a L. boulardi RhoGAP seemingly required for parasitism success on a resistant D. melanogaster strain (Labrosse et al., 2005b; Colinet et al., 2010), and LbSPNm, a serpin of L. boulardi closely related to a member of this family involved in suppressing the Drosophila immune defense (Colinet et al., 2009). This suggests that at least some of the discriminating reference bands contain proteins that may be involved in virulence. The observed variation in their amount in the venom is thus unlikely to be neutral, especially since the production of many venom proteins is probably energetically expensive (Morgenstern and King, 2013).
In this study, we illustrate how the method developed by (Mathé-Hubert et al., 2015) can be used for the venomics of populations or more generally for population proteomics. It shows for the first time that parasitoid wasp venom can be used to discriminate populations, including geographically very close ones. This is accurately demonstrated for L. boulardi whereas, differentiation among populations of L. heterotoma might be partly due to drift in the laboratory.
We also illustrated the use of the multivariate QST analysis developed by Martin et al. (2008) on the data generated by the method of (Mathé-Hubert et al., 2015) to detect the signature of natural selection. This signature is strong for both L. heterotoma and L. boulardi. For the latter, the comparison with the differentiation observed for neutral markers revealed that differentiation was partly due to some local adaptation.
Using proteomic data, we were able to tentatively identify the proteins contained in the bands that discriminate populations and therefore whose quantity might have evolved under heterogeneous selection. This identified several proteins known or suspected to be involved in virulence, but also proteins that have no identified putative function and would therefore deserve further investigation.
The combination of these two methods appears to be a powerful framework for population proteomics (Biron et al., 2006) and for the study of eco-evolutionary feedbacks between the proteomic level and the population and ecosystem levels (Diz et al., 2012). This is of interest not only to study evolution in the field at an intermediate level between the genome and phenotypes, or to understand the role of evolution in chemical ecology, but also for more applied issues. For example, biological control agents, sampled in the field in their area of origin, are often reared in the laboratory on a host species other than the target one (van Lenteren, 2003). The framework presented here could be used to monitor their potential unwanted adaptation to the laboratory host.
The data generated from the gel pictures analysis are available in the Supplementary Material S2. The raw gel pictures will be made available by the authors, upon request, without undue reservation.
DC, HM-H, J-LG, and MP conceived and designed research. JV was responsible for the sampling and rearing or the L. heterotoma populations. LK performed the analysis of the venom protein composition and data acquisition. HM-H designed and performed the statistical analyses. ÉD did the estimate of the FST. HM-H, MP, DC, J-LG, ÉD, and JV performed the writing and editing of the manuscript. All authors read and approved the final manuscript.
This work received support from the French National Research Agency (CLIMEVOL project, ANR-08- BLAN-0231) and the Department of Plant Health and Environment from the French National Institute for Agricultural Research (INRA). It was performed in the context of the Investments for the Future LABEX SIGNALIFE: program reference ANR-11-LABX-0028. HM-H was funded by the Provence Alpes Côte d'Azur (PACA) region and the Department of Plant Health and Environment from the French National Institute for Agricultural Research (INRA).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We are very grateful to Patricia Gibert and Stéphanie Llopis for maintaining and providing the L. boulardi and L. heterotoma populations and for producing and sending us the individuals analyzed in this study. We also thank C. Rebuf for technical assistance.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00156/full#supplementary-material
Table S1. Sampling locations.
Figure S1. Workflow of the statistical analysis.
Figure S2. Analysis of the correlation between the reference bands for L. boulardi and L. heterotoma (panel A and B respectively).
Supplementary Material S1. R functions allowing to reproduce the Qst analysis (tests and graphics).
Supplementary Material S2. Data produced in this work.
Abdel-Rahman, M. A., Abdel-Nabi, I. M., El-Naggar, M. S., Abbas, O. A., and Strong, P. N. (2011). Intraspecific variation in the venom of the vermivorous cone snail Conus vexillum. Comp. Biochem. Physiol. C Toxicol. Pharmacol. 154, 318–325. doi: 10.1016/j.cbpc.2011.06.019
Abdel-Rahman, M. A., Omran, M. A., Abdel-Nabi, I. M., Ueda, H., and McVean, A. (2009). Intraspecific variation in the Egyptian scorpion Scorpio maurus palmatus venom collected from different biotopes. Toxicon 53, 349–359. doi: 10.1016/j.toxicon.2008.12.007
Alape-Giron, A., Sanz, L., Madrigal, M., Sasa, M., and Calvete, J. J. (2008). Snake venomics of the lancehead pitviper bothrops asper: geographic, individual, and ontogenetic variations. J. Proteome Res. 7, 3556–3571. doi: 10.1021/pr800332p
Albert, F. W., Treusch, S., Shockley, A. H., Bloom, J. S., and Kruglyak, L. (2014). Genetics of single-cell protein abundance variation in large yeast populations. Nature 506, 494–497. doi: 10.1038/nature12904
Allemand, R., Fleury, F., Lemaitre, C., and Bouletreau, M. (1999). Population dynamics and competitive interactions in two species of Leptopilina (Hymenoptera: Figitidae) which parasitize Drosophila in the Rhone valley (S-E France). Ann. Soc. Entomol. Fr. 35, 97–103.
Biron, D. G., Loxdale, H. D., Ponton, F., Moura, H., Marché, L., Brugidou, C., et al. (2006). Population proteomics: an emerging discipline to study metapopulation ecology. Proteomics 6, 1712–1715. doi: 10.1002/pmic.200500423
Bouétard, A., Côte, J., Besnard, A.-L., Collinet, M., and Coutellec, M.-A. (2014). Environmental versus anthropogenic effects on population adaptive divergence in the freshwater snail Lymnaea stagnalis. PLoS ONE 9:e106670. doi: 10.1371/journal.pone.0106670
Calvete, J. J., Sanz, L., Pérez, A., Borges, A., Vargas, A. M., Lomonte, B., et al. (2011). Snake population venomics and antivenomics of Bothrops atrox: paedomorphism along its transamazonian dispersal and implications of geographic venom variability on snakebite management. J. Proteomics 74, 510–527. doi: 10.1016/j.jprot.2011.01.003
Clobert, J., Le Galliard, J.-F., Cote, J., Meylan, S., and Massot, M. (2009). Informed dispersal, heterogeneity in animal dispersal syndromes and the dynamics of spatially structured populations. Ecol. Lett. 12, 197–209. doi: 10.1111/j.1461-0248.2008.01267.x
Colinet, D., Deleury, E., Anselme, C., Cazes, D., Poulain, J., Azema-Dossat, C., et al. (2013a). Extensive inter- and intraspecific venom variation in closely related parasites targeting the same host: the case of Leptopilina parasitoids of Drosophila. Insect Biochem. Mol. Biol. 43, 601–611. doi: 10.1016/j.ibmb.2013.03.010
Colinet, D., Dubuffet, A., Cazes, D., Moreau, S., Drezen, J.-M., and Poirié, M. (2009). A serpin from the parasitoid wasp Leptopilina boulardi targets the Drosophila phenoloxidase cascade. Dev. Comp. Immunol. 33, 681–689. doi: 10.1016/j.dci.2008.11.013
Colinet, D., Mathé-Hubert, H., Allemand, R., Gatti, J. L., and Poirié, M. (2013b). Variability of venom components in immune suppressive parasitoid wasps: from a phylogenetic to a population approach. J. Insect Physiol. 59, 205–212. doi: 10.1016/j.jinsphys.2012.10.013
Colinet, D., Schmitz, A., Cazes, D., and Gatti, J. (2010). The origin of intraspecific variation of virulence in an eukaryotic immune suppressive parasite. PLoS Pathog. 6:206. doi: 10.1371/journal.ppat.1001206
Colinet, D., Schmitz, A., Depoix, D., Crochard, D., and Poirié, M. (2007). Convergent use of RhoGAP toxins by eukaryotic parasites and bacterial pathogens. PLoS Pathog. 3:e203. doi: 10.1371/journal.ppat.0030203
Coulette, Q., Lemauf, S., Colinet, D., Prévost, G., Anselme, C., Poirié, M., et al. (2017). Biochemical characterization and comparison of aspartylglucosaminidases secreted in venom of the parasitoid wasps Asobara tabida and Leptopilina heterotoma. PLoS ONE 12:e0181940. doi: 10.1371/journal.pone.0181940
Delava, E., Allemand, R., Léger, L., Fleury, F., and Gibert, P. (2014). The rapid northward shift of the range margin of a Mediterranean parasitoid insect (Hymenoptera) associated with regional climate warming. J. Biogeogr. 41, 1379–1389. doi: 10.1111/jbi.12314
Dion, E., Zélé, F., Simon, J.-C., and Outreman, Y. (2011). Rapid evolution of parasitoids when faced with the symbiont-mediated resistance of their hosts. J. Evol. Biol. 24, 741–750. doi: 10.1111/j.1420-9101.2010.02207.x
Diz, A. P., Martínez-Fernández, M., and Rolán-Alvarez, E. (2012). Proteomics in evolutionary ecology: linking the genotype with the phenotype. Mol. Ecol. 21, 1060–1080. doi: 10.1111/j.1365-294X.2011.05426.x
Dubuffet, A., Alvarez, C. I. R., Drezen, J.-M., van Alphen, J. J. M., and Poirié, M. (2006). Do parasitoid preferences for different host species match virulence? Physiol. Entomol. 31, 170–177. doi: 10.1111/j.1365-3032.2006.00505.x
Dupas, S., Carton, Y., and Poirié, M. (2003). Genetic dimension of the coevolution of virulence-resistance in Drosophila - parasitoid wasp relationships. Heredity 90, 84–89. doi: 10.1038/sj.hdy.6800182
Dupas, S., Dubuffet, A., Carton, Y., and Poirié, M. (2009). Local, geographic and phylogenetic scales of coevolution in Drosophila-parasitoid interactions. Adv. Parasitol. 70, 281–295. doi: 10.1016/S0065-308X(09)70011-9
Fleury, F., Allemand, R., Fouillet, P., and Boulétreau, M. (1995). Genetic variation in locomotor activity rhythm among populations of Leptopilina heterotoma (Hymenoptera: Eucoilidae), a larval parasitoid of Drosophila species. Behav. Genet. 25, 81–89.
Fleury, F., Ris, N., Allemand, R., Fouillet, P., Carton, Y., and Boulétreau, M. (2004). Ecological and genetic interactions in Drosophila-parasitoids communities: a case study with D. melanogaster, D. simulans and their common Leptopilina parasitoids in south-eastern France. Genetica 120, 181–194. doi: 10.1023/B:GENE.0000017640.78087.9e
Francischetti, I. M., Gombarovits, M., Valenzuela, J., Carlini, C., and Guimaraes, J. (2000). Intraspecific variation in the venoms of the South American rattlesnake (Crotalus durissus terrificus). Comp. Biochem. Physiol. Part C 127, 23–36. doi: 10.1016/S0742-8413(00)00129-8
Gibert, P., Allemand, R., Henri, H., and Huey, R. B. (2010). Local adaptation and evolution of parasitoid interactions in an invasive species, Drosophila subobscura. Evol. Ecol. Res. 12, 873–883. Available online at: http://www.evolutionary-ecology.com/issues/v12/n07/hhar2606.pdf
Goecks, J., Mortimer, N. T., Mobley, J. A., Bowersock, G. J., Taylor, J., and Schlenke, T. A. (2013). Integrative approach reveals composition of endoparasitoid wasp venoms. PLoS ONE 8:e64125. doi: 10.1371/journal.pone.0064125
Goudet, J., and Büchi, L. (2006). The effects of dominance, regular inbreeding and sampling design on Q(ST), an estimator of population differentiation for quantitative traits. Genetics 172, 1337–1347. doi: 10.1534/genetics.105.050583
Hita, M., Espagne, E., Lemeunier, F., Pascual, L., Carton, Y., Periquet, G., et al. (2006). Mapping candidate genes for Drosophila melanogaster resistance to the parasitoid wasp Leptopilina boulardi. Genet. Res. 88, 81–91. doi: 10.1017/S001667230600841X
Hita, M. T., Poirié, M., Leblanc, N., Lemeunier, F., Lutcher, F., Frey, F., et al. (1999). Genetic localization of a Drosophila melanogaster resistance gene to a parasitoid wasp and physical mapping of the region. Genome Res. 9, 471–481. doi: 10.1101/gr.9.5.471
Hofmann, E. P., Rautsaw, R. M., Strickland, J. L., Holding, M. L., Hogan, M. P., Mason, A. J., et al. (2018). Comparative venom-gland transcriptomics and venom proteomics of four Sidewinder Rattlesnake (Crotalus cerastes) lineages reveal little differential expression despite individual variation. Sci. Rep. 8, 1–15. doi: 10.1038/s41598-018-33943-5
Holding, M. L., Biardi, J. E., Gibbs, H. L., and Holding, M. L. (2016). Coevolution of venom function and venom resistance in a rattlesnake predator and its squirrel prey. Proc. R. Soc. B Biol. Sci. 283:20152841. doi: 10.1098/rspb.2015.2841
Labrosse, C., Eslin, P., Doury, G., Drezen, J. M., and Poirié, M. (2005a). Haemocyte changes in D. melanogaster in response to long gland components of the parasitoid wasp Leptopilina boulardi: a Rho-GAP protein as an important factor. J. Insect Physiol. 51, 161–170. doi: 10.1016/j.jinsphys.2004.10.004
Labrosse, C., Stasiak, K., Lesobre, J., Grangeia, A., Huguet, E., Drezen, J. M., et al. (2005b). A RhoGAP protein as a main immune suppressive factor in the Leptopilina boulardi (Hymenoptera, Figitidae)-Drosophila melanogaster interaction. Insect Biochem. Mol. Biol. 35, 93–103. doi: 10.1016/j.ibmb.2004.10.004
Lavandero, B., and Tylianakis, J. M. (2013). Genotype matching in a parasitoid-host genotypic food web: an approach for measuring effects of environmental change. Mol. Ecol. 22, 229–238. doi: 10.1111/mec.12100
Lee, M. J., Kalamarz, M. E., Paddibhatla, I., Small, C., Rajwani, R., and Govind, S. (2009). Virulence factors and strategies of Leptopilina spp.: selective responses in Drosophila hosts. Adv. Parasitol. 70, 123–145. doi: 10.1016/S0065-308X(09)70005-3
Leinonen, T., McCairns, R. J. S., O'Hara, R. B., and Merilä, J. (2013). QST–FST comparisons: evolutionary and ecological insights from genomic heterogeneity. Nat. Rev. Genet. 14, 179–190. doi: 10.1038/nrg3395
Martin, G., Chapuis, E., and Goudet, J. (2008). Multivariate QST-FST comparisons: a neutrality test for the evolution of the G matrix in structured populations. Genetics 180, 2135–2149. doi: 10.1534/genetics.107.080820
Mathé-Hubert, H., Gatti, J.-L., Colinet, D., Poirié, M., and Malausa, T. (2015). Statistical analysis of the individual variability of 1D protein profiles as a tool in ecology: an application to parasitoid venom. Mol. Ecol. Resour. 15, 1120–1132. doi: 10.1111/1755-0998.12389
Moiroux, J., Baaren, J., Van Poyet, M., Couty, A., Eslin, P., and Le Roux, V. (2018). Response of life-history traits to artificial and natural selection for virulence and nonvirulence in a Drosophila parastitoid, Asobara tabida. Insect Sci. 25, 317–327. doi: 10.1111/1744-7917.12428
Moiroux, J., Delava, E., Fleury, F., and van Baaren, J. (2013). Local adaptation of a Drosophila parasitoid: habitat-specific differences in thermal reaction norms. J. Evol. Biol. 26, 1108–1116. doi: 10.1111/jeb.12122
Moreau, C., Bhérer, C., Vézina, H., Jomphe, M., Labuda, D., and Excoffier, L. (2011). Deep human genealogies reveal a selective advantage to be on an expanding wave front. Science 334, 1148–1150. doi: 10.1126/science.1212880
Moreau, S. J. M., Cherqui, A., Doury, G., Dubois, F., Fourdrain, Y., Sabatier, L., et al. (2004). Identification of an aspartylglucosaminidase-like protein in the venom of the parasitic wasp Asobara tabida (Hymenoptera: Braconidae). Insect Biochem. Mol. Biol. 34, 485–492. doi: 10.1016/j.ibmb.2004.03.001
Mrinalini, A. L. S., Wright, J., Martinson, E., Wheeler, D., and Werren, J. H. (2014). Parasitoid venom induces metabolic cascades in fly hosts. Metabolomics 11, 350–366. doi: 10.1007/s11306-014-0697-z
Nyabuga, F. N., Loxdale, H. D., Heckel, D. G., and Weisser, W. W. (2011). Temporal genetic structuring of a specialist parasitoid, Lysiphlebus hirticornis Mackauer (Hymenoptera: Braconidae) attacking a specialist aphid on tansy. Biol. J. Linn. Soc. 102, 737–749. doi: 10.1111/j.1095-8312.2011.01620.x
Perez-Riverol, A., Roberto, J., Musacchio, A., Sergio, M., and Brochetto-braga, M. R. (2017). Wasp venomic: unravelling the toxins arsenal of Polybia paulista venom and its potential pharmaceutical applications. J. Proteomics 161, 88–103. doi: 10.1016/j.jprot.2017.04.016
Poisot, T., Bever, J. D., Nemri, A., Thrall, P. H., and Hochberg, M. E. (2011). A conceptual framework for the evolution of ecological specialisation. Ecol. Lett. 14, 841–851. doi: 10.1111/j.1461-0248.2011.01645.x
Porth, I., Klápste, J., McKown, A. D., La Mantia, J., Guy, R. D., Ingvarsson, P. K., et al. (2015). Evolutionary quantitative genomics of Populus trichocarpa. PLoS ONE 10:e0142864. doi: 10.1371/journal.pone.0142864
Ravigné, V., Dieckmann, U., and Olivieri, I. (2009). Live where you thrive: joint evolution of habitat choice and local adaptation facilitates specialization and promotes diversity. Am. Nat. 174, E141–E169. doi: 10.1086/605369
Ris, N. (2003). Hétérogénéité spatiale, plasticité phénotypique et trade-off environnementaux: rôle de l'espèce hôte et de la température dans la différenciation génétique des populations du parasitoïde Leptopilina heterotoma (Hymenoptera).
Rodríguez-ravelo, R., Coronas, F. I. V., Zamudio, F. Z., González-morales, L., López, G. E., Urquiola, A. R., et al. (2013). The Cuban scorpion Rhopalurus junceus (Scorpiones, Buthidae): component variations in venom samples collected in different geographical areas. J. Venom. Anim. Toxins Incl. Trop. Dis. 19, 1–10. doi: 10.1186/1678-9199-19-13
Rouchet, R., and Vorburger, C. (2014). Experimental evolution of parasitoid infectivity on symbiont-protected hosts leads to the emergence of genotype specificity. Evolution 68, 1607–1616. doi: 10.1111/evo.12377
Seyahooei, M. A., van Alphen, J. J. M., and Kraaijeveld, K. (2011). Genetic structure of Leptopilina boulardi populations from different climatic zones of Iran. BMC Ecol. 11:4. doi: 10.1186/1472-6785-11-4
Simoes, P., Pascual, M., Coelho, M. M., and Matos, M. (2010). Divergent evolution of molecular markers during laboratory adaptation in Drosophila subobscura. Genetica 138, 999–1009. doi: 10.1007/s10709-010-9486-4
Travis, J. M. J., Münkemüller, T., Burton, O. J., Best, A., Dytham, C., and Johst, K. (2007). Deleterious mutations can surf to high densities on the wave front of an expanding population. Mol. Biol. Evol. 24, 2334–2343. doi: 10.1093/molbev/msm167
van Baaren, J., Boivin, G., Bourdais, D., and Roux, O. (2007). “Antennal sensilla of hymenopteran parasitic wasps: variations linked to host exploitation behavior,” in Modern Research and Educational Topics in Microscopy, eds. A. Méndez-Vilas and J. Díaz (Badajoz: Elsevier), 345–352.
van Baaren, J., Dufour, C. M. S., Pierre, J.-S., Martel, V., and Louapre, P. (2016). Evolution of life-history traits and mating strategy in males: a case study on two populations of a Drosophila parasitoid. Biol. J. Linn. Soc. 117, 231–240. doi: 10.1111/bij.12644
van Lenteren, J. C. (2003). “Need for quality control of mass- produced biological control agents,” in Quality Control and Production of Biological Control Agents Theory and Testing Procedures, ed. J. C. van Lenteren (Wageningen: CABI Publishing), 1–18.
Vayssade, C., Martel, V., Moiroux, J., Fauvergue, X., Van Alphen, J. J. M., and van Baaren, J. (2012). The response of life-history traits to a new species in the community: a story of Drosophila parasitoids from the Rhône and Saône valleys. Biol. J. Linn. Soc. 107, 153–165. doi: 10.1111/j.1095-8312.2012.01918.x
Vinchon, S., Moreau, S. J. M., Drezen, J. M., Prévost, G., and Cherqui, A. (2010). Molecular and biochemical analysis of an aspartylglucosaminidase from the venom of the parasitoid wasp Asobara tabida (Hymenoptera: Braconidae). Insect Biochem. Mol. Biol. 40, 38–48. doi: 10.1016/j.ibmb.2009.12.007
Vuarin, P., Allemand, R., Moiroux, J., van Baaren, J., and Gibert, P. (2012). Geographic variations of life history traits and potential trade-offs in different populations of the parasitoid Leptopilina heterotoma. Naturwissenschaften 99, 903–912. doi: 10.1007/s00114-012-0972-7
Walker, A. A., Robinson, S. D., Yeates, D. K., Jin, J., Baumann, K., Dobson, J., et al. (2018). Entomo-venomics: the evolution, biology and biochemistry of insect venoms. Toxicon 154, 15–27. doi: 10.1016/j.toxicon.2018.09.004
Keywords: adaptive divergence, local adaptation, multivariate QST, antagonistic coevolution, individual 1D SDS-PAGE, population proteomics, Leptopilina wasps, Drosophila
Citation: Mathé-Hubert H, Kremmer L, Colinet D, Gatti J-L, Van Baaren J, Delava É and Poirié M (2019) Variation in the Venom of Parasitic Wasps, Drift, or Selection? Insights From a Multivariate QST Analysis. Front. Ecol. Evol. 7:156. doi: 10.3389/fevo.2019.00156
Received: 01 March 2019; Accepted: 23 April 2019;
Published: 10 May 2019.
Edited by:Kartik Sunagar, Indian Institute of Science (IISc), India
Reviewed by:Mrinalini, National University of Singapore, Singapore
Yehu Moran, Hebrew University of Jerusalem, Israel
Darin Rokyta, State College of Florida, United States
Micaiah Ward, Florida State University, United States, in collaboration with reviewer DR
Copyright © 2019 Mathé-Hubert, Kremmer, Colinet, Gatti, Van Baaren, Delava and Poirié. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hugo Mathé-Hubert, email@example.com
†Present Address: Hugo Mathé-Hubert, LIEC UMR, Université de Lorraine, CNRS, Metz, France
In Memoriam: This paper is dedicated to the memory of Prof. Roland Allemand