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

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.

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.

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

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

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.

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.

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.

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