Variation in the Venom of Parasitic Wasps, Drift, or Selection? Insights From a Multivariate QST Analysis

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 Q_ST 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.

Introduction

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 Q_ST 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 (F_ST). One of these approaches, the multivariate Q_ST 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 Q_ST 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 Q_ST 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 S^te Foy-Lès-Lyon and S^t 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 S^te 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

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 (S^te 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

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, S^te Foy-Lès-Lyon; SLA, S^t 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 Q_ST 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 Q_ST analysis to disentangle between drift, heterogeneous selection (local adaptation) and homogeneous selection (Figures S1A–C). The Q_ST 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 Q_ST analysis was not performed on the additive genetic variance (strict-sense Q_ST) but on the variation between individuals raised in a common environment i.e., the total genetic variance (broad-sense Q_ST). This broad-sense Q_ST may be biased downward which is further described in the discussion. The multivariate Q_ST 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.H₀: 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 Q_ST 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.H₀ The proportionality coefficient $\rho =\frac{2{\text{F}}_{\text{ST}}}{1-{\text{F}}_{\text{ST}}}$

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 F_ST). Information on such neutral differentiation was only available for three of the four populations of L. boulardi studied (Avignon, S^te Foy-Lès-Lyon, and S^t 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 ρ = $\frac{2F_{ST}}{1-F_{ST}}$. This is typically tested by checking for an overlap between the confidence intervals of ρ and that of $\frac{2F_{ST}}{1-F_{ST}}$ (Martin et al., 2008). The only information we had on neutral differentiation was a point estimate of the L. boulardi F_ST (= 0.084). However, we did not need the confidence interval of the F_ST because it turned out that the point estimate of $\frac{2F_{ST}}{1-F_{ST}}$ itself fall within the confidence interval of ρ. This estimation of the F_ST 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.

Results

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. S^te Foy-Lès-Lyon and Éyguières vs. S^t Laurent d'Agny) and one of the two comparisons of populations of the same area (S^te Foy-Lès-Lyon vs. S^t 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).

TABLE 1

Table 1. Venom-based discrimination between L. boulardi populations.

TABLE 2

Table 2. Venom-based discrimination between L. heterotoma populations.

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 Q_ST 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 $\frac{2F_{ST}}{1-F_{ST}}$.

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 (S^te Foy-Lès-Lyon and S^t 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

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 Q_ST analysis yielded a highly significant result (for both L. boulardi and L. heterotoma p < 0.001; H₀: 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 Q_ST of 0.108 [0.081; 0.165]; L. heterotoma: ρ = 0.738; ρϵ[0.567;1.06], which corresponds to a Q_ST 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 F_ST to perform the 2nd test $\left(H_0:\rho =\frac{2F_{ST}}{1-F_{ST}}\right)$ 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 F_ST = 0.084 gave an estimate of $\frac{2F_{ST}}{1-F_{ST}}=0.183$, 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 S^te 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

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

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

Table 3. Correspondence between L. boulardi bands that discriminate between populations and their putative protein content.

TABLE 4

Table 4. Correspondence between L. heterotoma bands that discriminate between the laboratory maintained populations and their putative protein content.

Discussion

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 Q_ST analysis that partitions the variance of a quantitative trait (Q) among subpopulations (S) relative to the total population (T), using the multivariate Q_ST analysis developed by Martin et al. (2008). In the strict sense, Q_ST analyses should be based on additive genetic variance, but in practice, they also use Q_ST 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 P_ST 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 Q_ST analyses which assume the lack of non-additive genetic variances, an approximation that can bias the Q_ST estimates. This bias is generally considered low because, in most cases, Q_ST 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 Q_ST, 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 Q_ST 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?

Leptopilina boulardi

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 S^te Foy-Lès-Lyon and S^t 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 Q_ST 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 $\frac{2F_{ST}}{1-F_{ST}}$. Secondly, most selection pressures could be heterogeneous among populations, which would push them toward different optimal values $\left(\rho >\frac{2F_{ST}}{1-F_{ST}}\right)$. 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 $\left(\rho \approx \frac{2F_{ST}}{1-F_{ST}}\right)$. 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 $\frac{2F_{ST}}{1-F_{ST}}$ suggests that both are of comparable strength.

Leptopilina heterotoma

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 F_ST = 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 Q_ST 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 F_ST 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).

Conclusion

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 Q_ST 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.

Data Availability

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.

Author Contributions

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 F_ST. 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.

Funding

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.

Acknowledgments

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.

Supplementary Material

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.

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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.

Google Scholar

Binford, G. J. (2001). An analysis of geographic and intersexual chemical variation in venoms of the spider Tegenaria agrestis (Agelenidae). Toxicon 39, 955–968. doi: 10.1016/S0041-0101(00)00234-8

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Brommer, J. E. (2011). Whither Pst? The approximation of Qst by Pst in evolutionary and conservation biology. J. Evol. Biol. 24, 1160–1168. doi: 10.1111/j.1420-9101.2011.02268.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Dennis, A. B., Patel, V., Oliver, K. M., and Vorburger, C. (2017). Parasitoid gene expression changes after adaptation to symbiont-protected hosts. Evolution 71, 2599–2617. doi: 10.1111/evo.13333

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Dray, S., and Dufour, A. B. (2007). The ade4 package: implementing the duality diagram for ecologists. J. Stat. Softw. 22, 1–18. doi: 10.18637/jss.v022.i04

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Dupas, S., and Carton, Y. (1999). Two non-linked genes for specific virulence of Leptopilina boulardi against Drosophila melanogaster and D. yakuba. Evol. Ecol. 13, 211–220.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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.

PubMed Abstract | Google Scholar

Fleury, F., Gibert, P., Ris, N., and Allemand, R. (2009). Ecology and life history evolution of frugivorous Drosophila parasitoids. Adv. Parasitol. 70, 3–44. doi: 10.1016/S0065-308X(09)70001-6

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Fronhofer, E. A., Kubisch, A., Hilker, F. M., Hovestadt, T., and Poethke, H. J. (2012). Why are metapopulations so rare? Ecology 93, 1967–1978. doi: 10.2307/23225197.

PubMed Abstract | CrossRef Full Text | Google Scholar

Futuyma, D., and Moreno, G. (1988). The evolution of ecological specialization. Annu. Rev. Ecol. Syst. 19, 207–233.

Google Scholar

Gandon, S. (2002). Local adaptation and the geometry of host – parasite coevolution. Ecol. Lett. 5, 246–256. doi: 10.1046/j.1461-0248.2002.00305.x

CrossRef Full Text | Google Scholar

Gandon, S., and Nuismer, S. L. (2009). Interactions between genetic drift, gene flow, and selection mosaics drive parasite local adaptation. Am. Nat. 173, 212–224. doi: 10.1086/593706

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Gomulkiewicz, R., Thompson, J., Holt, R., Nuismer, S., and Hochberg, M. (2000). Hot spots, cold spots, and the geographic mosaic theory of coevolution. Am. Nat. 156, 156–174. doi: 10.1086/303382

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Goudet, J., and Martin, G. (2007). Under neutrality, QST ≤ FST when there is dominance in an island model. Genetics 176, 1371–1374. doi: 10.1534/genetics.106.067173

CrossRef Full Text | Google Scholar

Hallatschek, O., and Nelson, D. R. (2010). Life at the front of an expanding population. Evolution 64, 193–206. doi: 10.1111/j.1558-5646.2009.00809.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hapfelmeier, A., and Ulm, K. (2013). A new variable selection approach using random forests. Comput. Stat. Data Anal. 60, 50–69. doi: 10.1016/j.csda.2012.09.020

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Hothorn, T., Hornik, K., and Zeileis, A. (2006). Unbiased recursive partitioning: a conditional inference framework. J. Comput. Graph. Stat. 15, 651–674. doi: 10.1198/106186006X133933

CrossRef Full Text | Google Scholar

Kawecki, T. J., and Ebert, D. (2004). Conceptual issues in local adaptation. Ecol. Lett. 7, 1225–1241. doi: 10.1111/j.1461-0248.2004.00684.x

CrossRef Full Text | Google Scholar

Kraaijeveld, A. R., Van Alphen, J., and Godfray, H. (1998). The coevolution of host resistance and parasitoid virulence. Parasitology 116, 29–45.

PubMed Abstract | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinson, E. O., and Werren, J. H. (2018). Venom is beneficial but not essential for development and survival of Nasonia. Ecol. Entomol. 43, 146–153. doi: 10.1111/een.12480

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreau, S. J., and Asgari, S. (2015). Venom proteins from parasitoid wasps and their biological functions. Toxins 7, 2385–2412. doi: 10.3390/toxins7072385

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Morgenstern, D., and King, G. F. (2013). The venom optimization hypothesis revisited. Toxicon 63, 120–128. doi: 10.1016/j.toxicon.2012.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrissey, J. H. (1981). Silver stain for proteins in polyacrylamide gels: a modified procedure with enhanced uniform sensitivity. Anal. Biochem. 117, 307–10.

PubMed Abstract | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Poirié, M., Carton, Y., and Dubuffet, A. (2009). Virulence strategies in parasitoid Hymenoptera as an example of adaptive diversity. C. R. Biol. 332, 311–320. doi: 10.1016/j.crvi.2008.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Poirié, M., Colinet, D., and Gatti, J. L. (2014). Insights into function and evolution of parasitoid wasp venoms. Curr. Opin. Insect Sci. 6, 52–60. doi: 10.1016/j.cois.2014.10.004

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rauch, G., and Weisser, W. W. (2007). Local and spatial dynamics of a host–parasitoid system in a field experiment. Basic Appl. Ecol. 8, 89–95. doi: 10.1016/j.baae.2006.03.007

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Remigio, E. A., and Duda, J. R. (2008). Evolution of ecological specialization and venom of a predatory marine gastropod. Mol. Ecol. 17, 1156–1162. doi: 10.1111/j.1365-294X.2007.03627.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rouchet, R., and Vorburger, C. (2012). Strong specificity in the interaction between parasitoids and symbiont-protected hosts. J. Evol. Biol. 25, 2369–2375. doi: 10.1111/j.1420-9101.2012.02608.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Strobl, C., Boulesteix, A.-L., Kneib, T., Augustin, T., and Zeileis, A. (2008). Conditional variable importance for random forests. BMC Bioinform. 9:307. doi: 10.1186/1471-2105-9-307

PubMed Abstract | CrossRef Full Text | Google Scholar

Thompson, J. (1997). Evaluating the dynamics of coevolution among geographically structured populations. Ecology 78, 1619–1623.

Google Scholar

Thompson, J. N. (1994). The Coevolutionary Process. Chicago, IL: University of Chicago Press.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

van Baaren, J., and Boivin, G. (1998). Genotypic and kin discrimination in a solitary Hymenopterous parasitoid: implications for speciation. Evol. Ecol. 12, 523–534.

Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, A. (2012). The role of robustness in phenotypic adaptation and innovation. Proc. Biol. Sci. 279, 1249–1258. doi: 10.1098/rspb.2011.2293

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Weisser, W. (2000). Metapopulation dynamics in an aphid-parasitoid system. Entomol. Exp. Appl. 97, 83–92. doi: 10.1046/j.1570-7458.2000.00719.x

CrossRef Full Text | Google Scholar

Whitlock, M. (1999). Neutral additive genetic variance in a metapopulation. Genet. Res. 74, 215–221.

PubMed Abstract | Google Scholar

Wuster, W., Otsuka, S., Malhotra, A., and Thorpe, R. S. (1992). Population systematics of Russell's viper: a multivariate study. Biol. J. Linn. Soc. 47, 97–113.

Google Scholar