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, interindividual variations of venom composition were also recently evidenced in strains and populations of species of the Leptopilina and Psyttalia genera (Colinet et al., 2013a,b). This variation in the venom composition may also be at the origin of the rapid evolution of parasitoid virulence observed under experimental evolution (Dion et al., 2011; Rouchet and Vorburger, 2014). Indeed, differential selection of the expression levels of genes encoding venom proteins was observed during the evolution of parasitoid virulence in this last experimental evolution (Dennis et al., 2017). Nevertheless, neutral molecular variation is common, and some of the variations in the venom composition are therefore most likely neutral and do not affect virulence. This raises the question of whether parasitoid populations could be differentiated according to the composition of the venom and which of these differences would be caused by drift or local adaptation.

Here, we have investigated the population structure of L. boulardi and L. heterotoma in the Rhône-Saône valley according to the individual composition of the venom. For L. boulardi, we have studied four populations, two in the south of the Rhône-Saône valley that this species invaded more than 40 years ago (Delava et al., 2014) and two at the northern end of this valley (close to Lyon) that was invaded by this species between 1993 and 2003 (Delava et al., 2014). For L. heterotoma, five populations were studied, one in the south, outside the Rhône-Saône valley (200 km east) where L. boulardi predominates, two in the middle of the Rhône-Saône valley where it appeared between 1993 and 2003, and two in the Saône valley, north of Lyon, where its presence was detected between 2003 and 2011 (Delava et al., 2014).The analysis of the venom composition of individual wasps was performed using a recently developed method based on 1D electrophoresis and further digital and statistical analysis of the intensity of the protein bands (Mathé-Hubert et al., 2015).

We report for the first time a structuration of parasitoid populations based on venom composition, a trait often directly related to the outcome of host-parasitoid interactions. This differentiation among populations in venom composition is partly induced by drift but may also be due to some heterogeneity in the selection pressures experienced by different populations. We used the multivariate QST approach developed by Martin et al. (2008) to detect whether there is selection for the same optimal relative amounts of venom components in all populations (homogeneous selection) or selection for different local optima (heterogeneous selection, or local adaptation), possibly in relation with some variation in the host physiology or resistance. For L. boulardi, data showed that part of the venom-based population structure is induced by a certain heterogeneity in selection within different populations, i.e., some local adaptation. For L. heterotoma, we could only perform part of the analysis, which evidenced a high evolutionary Mathé-Hubert et al. Parasitic Wasp Venom Locally Adapts

potential of the venom composition. These results are compared with those previously obtained for various parasitoid traits in the Rhône-Saône valley and are discussed in the context of coevolutionary theories.

# MATERIALS AND METHODS

# Sampling and Analyzed Individuals

The four L. boulardi populations were sampled in orchards in September 2010. Sampling locations were Ste Foy-Lès-Lyon and St Laurent d'Agny for northern populations (25 and 29 females sampled, respectively) and Éyguières and Avignon for southern populations (20 and 11 females sampled, respectively). Southern and northern L. boulardi populations are distant from 200 km (**Figure 1** and **Table S1**). Each field-collected individual was isolated and allowed to parasitize larvae of a D. melanogaster strain sampled in Ste Foy-Lès-Lyon. The offspring of these individuals were then stored separately at −80◦C. The venom of one daughter was analyzed per field-sampled female.

The five populations of L. heterotoma were sampled in orchards in 2008 (n = 15–20) and then maintained on D. melanogaster (Ste Foy-Lès-Lyon) until the beginning of September 2013, when individuals were stored at −80◦C. As discussed later, we cannot exclude that some of the evolution we detected occurred in the laboratory between 2008 and 2013. For this reason, we will refer to these populations as "laboratory-maintained populations." The five laboratorymaintained 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 R 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

between 2003 and 2011.

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

FIGURE 2 | Example of SDS-PAGE gels for venom analysis. Example of SDS-PAGE gels used to characterize the variation in the venom composition of L. boulardi (A) and L. heterotoma (B). Each lane contains the venom of a single wasp whose population of origin is indicated at the top of the lane (Éy, Éyguières; Av, Avignon; SFL, Ste Foy-Lès-Lyon; SLA, St Laurent d'Agny; Mo, Montbellet; Ve, Vence; So, Sonnay; Ép, Épinouze; Uch, Uchizy). The reference bands that significantly discriminate between populations after a Bonferroni correction are indicated.

coding for a higher or a lower expression of different venom proteins may indeed induce some correlation between the bands containing the two proteins. Since the linkage disequilibrium between traits should be low for the multivariate QST analysis (Martin et al., 2008), we merged reference bands that were too highly correlated. These set of correlated reference bands were detected using a UPGMA dendrogram in which the distance between pairs of reference bands was measured as 1 − |r| where r is the Pearson correlation coefficient (**Figure S2**). Then, we merged the bands whose correlation was higher than 0.6 in a composite reference band estimated as the mean of the merged bands. For L. boulardi, we built four composite reference bands from reference bands (16, 17), (19, 20), (22, 28), and (30, 31, 32) (**Figure S2A**). For L. heterotoma, we built six composite reference bands from reference bands (1–3), (6, 7), (8, 9), (15, 21), (19, 20), and (28, 29) (**Figure S2B**).

# Statistical Analyses of the Population Structuration

We performed all statistical analyses in R 3.4.4. The workflow of these analyses is described in **Figure S1**.

### Power of the Venom to Discriminate Populations

For each pair of populations, we used a Random Forest (RF) clustering algorithm to predict the population of origin of individuals based on the composition of their venom (the centered and scaled intensities of the reference bands). We conducted the RFs using the package party to account for the multicollinearity and high number of predictors compared to the low number of individuals per population (Hothorn et al., 2006; Strobl et al., 2008). For each RF, we have estimated the discriminating power of the venom composition as the proportion of accurately classified individuals. The significance of this discriminatory power was estimated using 10,000 permutations. To avoid overfitting, this was done using predictions based on the out-of-the-bag data. For all RFs fitted to L. boulardi, we used the weight argument to ensure that all populations have the same weight albeit the headcounts were not balanced.

### Effect of Selection on the Variation in Venom Composition Within and Among Populations

We used a QST analysis to disentangle between drift, heterogeneous selection (local adaptation) and homogeneous selection (**Figures S1A–C**). The QST compares the variance of a quantitative trait (Q) among the subpopulations (S) to the total variance (T).

As we had no information on the genetic relatedness between the analyzed individuals, this QST analysis was not performed on the additive genetic variance (strictsense QST) but on the variation between individuals raised in a common environment i.e., the total genetic variance (broad-sense QST). This broad-sense QST may be biased downward which is further described in the discussion. The multivariate QST analysis developed by Martin et al. (2008) relies on two complementary tests, the first being more powerful when different traits (here venom components) are subjected to different types of selection regime (i.e., homogeneous, heterogeneous, or neutral), the second when they all undergo the same selection pressures (either homogeneous or heterogeneous). For more details see **Figure S1C**. The R code used to perform these tests was developed by Martin et al. (2008). To ensure its long-lasting accessibility, it is also included in the **Supplementary Material S1**.

## **Test 1.H0: Covariances among (D) and within (G) populations are proportional (D** =ρ×**G)**

For both species, we first summarized the total variation in the venom composition using a principal components analysis (PCA; Dray and Dufour, 2007) in order to visualize how this variation is partitioned within and among populations. Then, we plotted the first two principal components (PC) of the PCA using the R function s.class (package ade4) that groups individuals of the same population. This allows to determine whether, as expected under neutrality, the dimension with the greatest variation within each population is also the one with the greatest variation between populations. To improve this visual comparison of the variation within vs. among populations, we projected on these plots the first PC of a PCA performed on the variation within populations and the first PC of a PCA performed on the variation among populations. These two additional PCA were performed using the functions wca and bca (package ade4), respectively. If the variation among populations is only induced by drift, then these axes of variation within and among populations should be mostly collinear. We measured this collinearity with the Pearson correlation coefficient. A function performing this analysis was added to the **Supplementary Material S1** (see also **Figure S1A**).

This expectation of proportionality between the variations within and among populations is the null hypothesis of the first test of the multivariate QST analysis developed by Martin et al. (2008). Indeed, under the null hypothesis that the components of venom are neutral, all the variations observed among populations appeared because the isolated populations drifted independently of each other. In each population, for each trait (here venom components), the amount of drift is proportional to the amount of variance within the population. This gives the statistical expectation that the average matrix of covariances within populations (G) is proportional to the matrix of covariances among populations (D): D = ρ × G, where the two matrices D and G are estimated by a MANOVA and the coefficient ρ of proportionality is estimated by a function developed by Martin et al. (2008). As the sample size per population is small, the pvalues were computed using the correction based on the Bartlett adjustment of likelihood ratios (Martin et al., 2008). In this analysis, the confidence interval of the proportionality coefficient ρ is estimated by maximum-likelihood (for more details, see Martin et al., 2008). This confidence interval is used for the second test.

### **Test 2.H0 The proportionality coefficient** ρ = **2FST 1**−**FST**

The second test relies on the comparison between the average differentiation observed for the traits of interest (measured by the proportionality coefficient ρ) and the differentiation observed for some neutral markers (measured by the FST). Information on such neutral differentiation was only available for three of the four populations of L. boulardi studied (Avignon, Ste Foy-Lès-Lyon, and St Laurent d'Agny) and it was not available for L. heterotoma. This test was therefore performed only on these three L. boulardi populations. Under the null hypothesis that, on average, the traits used to compute ρ are neutral, then ρ = 2FST 1−FST . This is typically tested by checking for an overlap between the confidence intervals of ρ and that of 2FST 1−FST (Martin et al., 2008). The only information we had on neutral differentiation was a point estimate of the L. boulardi FST (= 0.084). However, we did not need the confidence interval of the FST because it turned out that the point estimate of 2FST 1−FST itself fall within the confidence interval of ρ. This estimation of the FST is based on the RAD sequencing of L. boulardi individuals from the Rhône–Saône valley which led to the genotyping of about 16 individuals per population (the three studied populations) on 484 loci (Delava, in preparation).

## Identification of the Venom Components Possibly Under Heterogeneous Selection

The identification of venom components possibly under heterogeneous selection was done in two steps. First, we detected the reference bands most likely to be the ones under heterogeneous selection, and then we identified the components of the venom they contain.

For the first step, we used a clustering RF algorithm to predict, for both L. heterotoma and L. boulardi, the population of origin of the individuals based on the composition of their venom (the centered and scaled intensities of the reference bands). We conducted the RFs and computed the conditional importance of each explanatory variable (i.e., reference bands) using the package party (Hothorn et al., 2006; Strobl et al., 2008). We then tested the significance of each explanatory variable by comparing its conditional importance, its estimated null distribution using 10,000 permutations (Hapfelmeier and Ulm, 2013). The p-values resulting from this analysis were then Bonferroni corrected and we selected for further investigation the reference bands that were still significant.

We tentatively identified the proteins present in these reference bands by matching these bands with the bands on the 1D electrophoresis gels used for L. boulardi and L. heterotoma venom proteomics (Colinet et al., 2013a). Since one band can contain different proteins and a given protein can be found in different bands, the number of peptides matches from the mass spectrometry performed by Colinet et al. (2013a) was used to classify the proteins in the bands as abundant or not using a threshold of 10 peptides matches in the proteomic analysis. In these bands, low abundance proteins are unlikely to drive the observed population structure.

# 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. Ste Foy-Lès-Lyon and Éyguières vs. St Laurent d'Agny) and one of the two comparisons of populations of the same area (Ste Foy-Lès-Lyon vs. St Laurent d'Agny, two populations only 14 km apart). Furthermore, the classification was only a little more accurate when comparing distant

### TABLE 1 | Venom-based discrimination between L. boulardi populations.

For each pair of populations, the proportion of individuals accurately classified by the RF is given with the permutation-based p-value of the null hypothesis (in brackets). Sample sizes are indicated close to the location name.

For each pair of populations, the proportion of individuals accurately classified by the RF is given with the permutation-based p-value of the null hypothesis (in brackets). Sample sizes are indicated close to the location name.

populations than when comparing nearby populations (average accuracies: 81 and 74.1% respectively; **Table 1**). The RFs performed on each pair of L. heterotoma laboratory-maintained populations allowed to classify 85 to 97% of individuals by their venom (**Table 2**). These classifications were significant in all comparisons. The highest accuracy was achieved when comparing Vence to another population (average accuracies: 95.6%; **Table 2**), but, excluding Vence, there was almost no difference between comparisons involving distant populations or nearby (average accuracies: 90.6 and 87.5% respectively; **Table 2**).

# Effect of Selection on the Variation in Venom Composition Within and Among Populations

To detect whether selection affected variation in venom composition among and within populations, we used multivariate QST analysis, which performed two tests. The first investigated whether all traits were under the same selection regime (homogeneous, heterogeneous or neutral) by checking whether the (co)variances among populations were proportional to the (co)variances within populations. If such proportionality exists, then the dimensions that have the most variance among and within populations should be collinear. We provided a visual assessment of this part of the analysis based on PCAs. The second test investigated whether, on average, traits were under homogeneous or heterogeneous selection by testing whether the coefficient of proportionality (ρ) estimated in the first test was different from 2FST 1− FST .

For L. boulardi, 34.3% of the total variation in venom composition was summarized by the first two axes of the PCA. The axis summarizing the variation among populations discriminated mainly the populations of the North (Ste Foy-Lès-Lyon and St Laurent d'Agny) of those of the South (Avignon and Éyguières), although it also partly discriminated between these two southern populations. This axis was mainly collinear with the first axis of the global PCA. On the opposite, the axis summarizing the variation within populations was mainly collinear with the second axis of the global PCA (**Figure 3A**). This reveals that instead of the expected collinearity under neutrality, the two axes were mainly orthogonal (Pearson correlation of 0.09), thus showing that the reference bands did not all undergo the same type of selection (more details on **Figure S1C**).

A similar result was obtained for L. heterotoma. The first two axes of the PCA explained 39.9% of the total variation in venom composition, and the within and among populations axes were also predominantly orthogonal, with a Pearson correlation of 0.16. On the among populations axis, the population of Vence was completely separated from other populations and that of Montbellet was partially separated (**Figure 3B**).

In agreement with this visual analysis of the venom-based structure of populations, the first test of the QST analysis yielded a highly significant result (for both L. boulardi and L. heterotoma p < 0.001; H0: D = ρ × G where G and D are the within and among populations covariance matrices, respectively and ρ is the proportionality coefficient; see section Materials and Methods). This analysis estimated the proportionality coefficient and its 95% confidence interval (L. boulardi: ρ = 0.243; ρǫ [0.176; 0.395], which corresponds to a QST of 0.108 [0.081; 0.165]; L. heterotoma: ρ = 0.738; ρǫ [0.567; 1.06], which corresponds to a QST of 0.270 [0.221; 0.346]). The coefficient of proportionality of L. heterotoma was three times higher than that obtained for L. boulardi and the confidence intervals did not overlap, meaning that this difference was significant. For L. boulardi we also compared the estimated coefficient of proportionality (ρ) with the FST to perform the 2nd test H0 : ρ = 2FST 1−FST which allowed to check whether there was some homogeneous or heterogeneous selection (more details on **Figure S1C**). This test was not significant since the estimate of FST = 0.084 gave an estimate of 2FST 1−FST = 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

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.

bands which mainly discriminated Ste Foy-Lès-Lyon from other populations. Also, as observed on the PCA, Avignon was more strongly discriminated than Éyguières (bands B26 and B29), although this tended to be the opposite for bands B2 and B6.

For L. heterotoma, this allowed identification of 13 references bands that, in agreement with the PCA (**Figure 3B**), distinguished mainly between Vence and other laboratorymaintained 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.

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.

(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

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

Band contents were obtained from the 1D-SDS-PAGE proteomic study by Colinet et al. (2013a). Only proteins for which at least 10 peptide matches in Mascot searches against unisequences identified by transcriptomics of the venom apparatus were considered as abundant and used. The number of proteins in the band, their predicted function and the number of peptides matches is provided. NA, data not available.

studies of the geographical variation in venom composition was conducted on several venomous taxa (Wuster et al., 1992; Francischetti et al., 2000; Binford, 2001; Alape-Giron et al., 2008; Remigio and Duda, 2008; Abdel-Rahman et al., 2009, 2011; Calvete et al., 2011; Rodríguez-ravelo et al., 2013; Holding et al., 2016; Perez-Riverol et al., 2017; Hofmann et al., 2018), it was not the case for parasitoids. This is surprising since, in parasitoid wasps, venom should constrain the capacity to succeed on locally available hosts. Performing such a study is now facilitated by the development of a rapid and accurate method to analyse the quantitative variation of large sets of proteins from several individual samples (Mathé-Hubert et al., 2015). This method was used here to analyse the venom-based differentiation of four L. boulardi and five L. heterotoma samples from populations along a 300 km-long south—north gradient in the Rhône-Saône valley (South-East of France).

A major result of these analyses was that the venom composition, characterized by the relative proportions of the different constitutive proteins, allowed a significant discrimination of populations, even if only a few km apart. Such a differentiation among populations can evolve by drift when gene flow is low relative to drift or by selection when there are different local optima in each population (heterogeneous selection), whereas homogeneous selection tends to prevent evolution of differentiation among populations.

We further analyzed these differentiations with a QST analysis that partitions the variance of a quantitative trait (Q) among subpopulations (S) relative to the total population (T), using the multivariate QST analysis developed by Martin et al. (2008). In the strict sense, QST analyses should be based on additive genetic variance, but in practice, they also use QST in the broad-sense that applies to the total genetic variance (additive + epistasis + dominance; Bouétard et al., 2014; Porth et al., 2015), as well as the PST that considers the phenotypic (P) variance, thus also including the non-genetic environmental variance (Brommer, 2011).

In this study, we could exclude the existence of environmental variance since individuals were raised in a common environment, meaning that the phenotypic variance is equal to the total genetic variance. We therefore used broad-sense QST analyses which assume the lack of non-additive genetic variances, an approximation that can bias the QST estimates. This bias is generally considered low because, in most cases, QST is

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

Band contents were obtained from the 1D-SDS-PAGE proteomic study by Colinet et al. (2013a). Only proteins for which at least 10 peptide matches in Mascot searches against unisequences identified by transcriptomics of the venom apparatus were considered as abundant and used. The number of proteins in the band, their predicted function and the number of peptides matches is provided. NA, data not available.

computed for multigenic traits (Leinonen et al., 2013). However, the amount of a given venom component may depend on a low number of loci (although see Albert et al., 2014), so this bias could be more substantial. It generally induces an underestimation of the QST, both when it is due to neglected epistasis variance (Whitlock, 1999) and when due to neglected dominance variance (Goudet and Büchi, 2006; Leinonen et al., 2013) although this has been more controversial (Goudet and Martin, 2007). For a more detailed discussion on this bias, see Leinonen et al. (2013). Downwards biased QST estimates are too conservative for the detection of heterogeneous selection and too sensitive when detecting homogeneous selection, which is reassuring since, in such an analysis, heterogeneous selections are generally more interesting than homogeneous selections.

# Variation of the Venom Composition Among Populations: Drift or/and Selection?

### 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 Ste Foy-Lès-Lyon and St Laurent d'Agny, two northern populations only 14 km apart. This differentiation can originate from drift, heterogeneous selection, or a combination of both.

A strong drift could come from the northward extension of distribution of L. boulardi. Indeed, its repartition area has progressed 170 km northwards from 1993 to 2011 (Delava et al., 2014). The wavefront of expanding populations is known to drift significantly because the individuals that drive this expansion are sampled from previously sampled individuals (Travis et al., 2007; Hallatschek and Nelson, 2010). When the wavefront is expanding fast enough, the drift is expected to create a gradient of differentiation from the source to the wavefront (e.g., Moreau et al., 2011). However, only a slight North-South differentiation was observed for L. boulardi, and of the same order as the differentiation between nearby populations. Another explanation would be a metapopulation dynamic with frequent extinctions and recolonizations, due to a high migration rate. Such metapopulation dynamic is expected to create differentiation, even between close populations, since colonization events are associated with founder effects as already demonstrated (Weisser, 2000; Rauch and Weisser, 2007); Nyabuga et al., 2011.

Alternatively, the similar magnitude of differentiation between distant populations (North vs. South) and nearby populations could also be explained if: (i) the effect of drift on differentiations is low, either due to a low drift or a high migration rate and (ii) the heterogeneous selection is strong enough to counter the effect of drift and migration, thus creating some local adaptation. In keeping with this, the QST analysis detected some selection. Indeed, there was a highly significant lack of proportionality between variations within and among populations. Such a proportionality is expected in the absence of natural selection since if the variation among populations comes only from the drift occurring separately within each, the magnitude of drift should be proportional to the magnitude of the variation within populations. This lack of proportionality thus reveals the variation in the selection pressures experienced by the different venom components.

Such variation among traits (here venom components) can have different origins. While the observed variation is most likely neutral for some traits, the non-neutral variation of others may fall anywhere in a continuum described by the following three categories. Firstly, most selection pressures could be homogeneous among populations, which would keep all populations around the same optimal values and give a proportionality coefficient (ρ) lower than 2FST 1−FST . Secondly, most selection pressures could be heterogeneous among populations, which would push them toward different optimal values ρ > 2FST 1−FST . 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 ρ ≈ 2FST 1−FST . 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 2FST 1−FST 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 FST = 0.036; Ris, 2003). This leaves two non-mutually exclusive hypotheses for explaining the strength of this venom-based differentiation. Differentiation might be adaptive, resulting from some heterogeneity among populations in the optimal venom composition. Alternatively, all or part of the differentiation may have evolved by drift during the 5 years of laboratory rearing of the populations between field sampling and the venom analysis (see for example Simoes et al., 2010).

As with L. boulardi, there was a very significant lack of proportionality between variations within and among populations that could have been generated by heterogeneous or homogeneous selection or a combination of both. The second test would have allowed discriminating between these hypotheses but was not doable for L. heterotoma. An additional level of uncertainty comes from where the detected selection may have occurred. Indeed, the conditions of rearing in laboratory could have imposed certain modifications in the optimal composition of the venom to which the five laboratorymaintained populations could have adapted by homogeneous selection (van Lenteren, 2003).

In summary, for L. heterotoma, strong differentiation and lack of proportionality may reflect two non-mutually exclusive scenarios. Heterogeneous selection in the field alone could explain both. In addition, the combination of homogeneous laboratory and/or field selection and laboratory drift may explain the lack of proportionality and strong differentiation, respectively. Importantly, whatever the reasons for the observed differentiation and non-proportionality, results highlight the high amount of genetic variation in the venom composition of L. heterotoma.

## Local Adaptation of the Venom Composition?

For L. heterotoma, we could not test for local adaptation of the venom composition. However, several ecological features may have created it. Firstly, both spatial and seasonal variations in the relative abundances of Drosophila species have been described in the Rhône valley (Allemand et al., 1999; Ris, 2003). L. heterotoma is a generalist species whose venom may therefore contain either "generalist" virulence factors or combination of factors, each of which being effective for a subset of the host range. In the latter case, the local availability of the hosts could strongly influence the local composition of the venom. This variation in host availability is also likely affected by the competition with L. boulardi, both a specialist and a powerful competitor, which was absent from northern populations in 2003 and present in 2011. Its arrival may have induced changes in the exploitation pattern of Drosophila species by L. heterotoma, resulting in different selection pressures on venom components.

Interestingly, a strong latitudinal genetic differentiation of L. heterotoma was also previously reported, based on fitness traits (Fleury et al., 1995, 2004). However, this could be explained either by differences in abiotic conditions (e.g., temperature) or by a variation in the presence of L. boulardi as a competitor, or both. Another piece of information suggesting some heterogeneous selection could be the putatively adaptive North-South genetic differentiation for the time of preferred locomotor activity, with an ongoing shift between Lyon and Antibes that has been highlighted by Fleury et al. (1995). Indeed, Vence is not far from Antibes, and the middle populations of the Rhône Valley are not far from Lyon. Such a temperature-related difference in the period of activity could have other consequences on both the host and parasitoid physiology, which may affect the optimal venom composition.

For L. boulardi, the QST analysis detected some local adaptation of the composition of the venom, which could be induced by certain biotic or abiotic variations. Local adaptation could indeed be an optimization of something other than virulence if the expression of venom proteins is associated with other physiological traits. For example, if a venom protein is encoded by a gene (i) that is also expressed elsewhere than in the venom gland and (ii) that shows variation in cis-regulation of expression, then a local adaptation could be detected on the venom composition even though the selection takes place due to the protein effect in tissues other than the venom gland. However, such indirect selection on the composition of the venom should be rather rare because the changes it induces are likely not neutral if these venom proteins are costly to produce (Morgenstern and King, 2013).

Alternatively, this might reflect some variability in host availability. Indeed, in the studied area, L. boulardi parasitizes both D. melanogaster, predominant in the North, and D. simulans, predominant in the South (Fleury et al., 2004), and the optimal venom composition differ for success on both species (Cavigliasso et al., in preparation). The North could therefore be a hotspot of co-evolution with D. melanogaster and a coldspot with D. simulans, and vice versa for the South. However, this North-South variation in host availability cannot explain why venom composition differs between close populations. This variation between neighboring populations is the main surprise of this study and deserves to be analyzed with larger samples. As mentioned in the extinction/recolonization scenario, this differentiation could come from drift. If alternatively, it is a local adaptation of virulence, it may be an adaptation to an intra-specific variation in D. melanogaster, the main host of L. boulardi in this region (Fleury et al., 2004).

Such a "very local" adaptation may seem implausible given the low FST between these populations. However, local adaptation induced by coevolution differs from classical local adaptation. For example, since parasitoids and hosts interact antagonistically the local adaptation of one species in a given population is also the local maladaptation of the other species. The interactions between local adaptation and coevolution were studied by Gandon (2002) and Gandon and Nuismer (2009). They reveal that when genetic drift is substantial, higher migration rates actually promote local adaptation, as migrants increase the evolutionary potential (genetic diversity) of the population which is needed to overcome the effect of drift and win the arms race (Gandon, 2002). Also, although drift generally decreases local adaptation, it may sometimes have the opposite effect in the case of local adaptation of traits involved in an antagonistic co-evolutionary dynamic (Gandon and Nuismer, 2009). Indeed, drift occurring independently in each population of two species of a co-evolving "couple" may induce a certain spatial heterogeneity to which both species can adapt locally, which reinforces the spatial heterogeneity initiated by drift. Such simultaneous substantial amounts of drift within each population combined with migration among populations are precisely the predicted situation under the extinction/recolonization scenario we previously considered. Such scenario is not unlikely if we consider (i) the short generation time of these species, (ii) their high fertility, (iii) their ability to fly and disperse over long distances with the wind, (iv) the known seasonal variations in populations size, and finally, (v) the possibility of local extinctions or large fluctuations in the sizes of populations of hosts and parasitoids driven by the parasitic interaction, which would create a strong drift. These conditions are generally likely to create a metapopulation dynamic (Fronhofer et al., 2012).

If we combine this view of host-parasitoid metapopulation dynamic with the notion of informed dispersal (Clobert et al., 2009), then local adaptation between close populations for a trait as important as virulence may be expectable. As parasitoids can accurately estimate the suitability of a host before oviposition (van Baaren and Boivin, 1998; Dubuffet et al., 2006; van Baaren et al., 2007), their decision to settle or continue their migration, when emigrating, may depend on the prevalence of hosts that fit their needs. Such non-random dispersal strongly increases the possibility of local adaptation (Ravigné et al., 2009).

# Identifying the Venom Components Likely Under Heterogeneous Selection

To identify the traits under heterogeneous selection, we used the RF clustering algorithm to detect reference bands that discriminate significantly among populations. Indeed, heterogeneously selected venom components provide the greatest discriminating power among populations (high variance among populations relative to the variance within populations). We identified venom components (reference bands) that significantly discriminate between populations. However, even reference bands for which variation is neutral are expected to sometimes discriminate populations. We thus applied a Bonferroni correction to keep only the bands with the highest discriminating power. This selected 9 and 13 reference bands for L. boulardi and L. heterotoma, respectively. These two sets of reference bands should not be interpreted in the same way since we did not formally identify heterogeneous selection for L. heterotoma, which we did for L. boulardi.

The content of these reference bands was tentatively identified by manually matching them with the 1D electrophoresis gels analyzed by proteomics by Colinet et al. (2013a). This proteomics study characterized the venom content of L. boulardi and L. heterotoma lines, which are probably not representative of the whole diversity observed in the field. Nevertheless, there should be an overall correspondence between the content of the reference bands analyzed in this study and that of the bands analyzed by Colinet et al. (2013a). This identified 16 proteins as potentially under heterogeneous selection, six of which have no identified putative function. This highlights that our prior-less method can be used to discover novel active proteins. Indeed, venoms have strong and diverse biological activities but they contain a great diversity of proteins that might not all be worth studying. Our approach allows selecting without a priori proteins likely to contribute to parasitism success on the host. The other identified proteins are members of the RhoGAP and serpin families found in reference bands of both Leptopilina species, as well as an aspartylglucosaminidase (AGA) found in two L. heterotoma reference bands. Among the proteins identified in the discriminating reference bands are LbGAP, a L. boulardi RhoGAP seemingly required for parasitism success on a resistant D. melanogaster strain (Labrosse et al., 2005b; Colinet et al., 2010), and LbSPNm, a serpin of L. boulardi closely related to a member of this family involved in suppressing the Drosophila immune defense (Colinet et al., 2009). This suggests that at least some of the discriminating reference bands contain proteins that may be involved in virulence. The observed variation in their amount in the venom is thus unlikely to be neutral, especially since the production of many venom proteins is probably energetically expensive (Morgenstern and King, 2013).

# 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 QST analysis developed by Martin et al. (2008) on the data generated by the method of (Mathé-Hubert et al., 2015) to detect the signature of natural selection. This signature is strong for both L. heterotoma and L. boulardi. For the latter, the comparison with the differentiation observed for neutral markers revealed that differentiation was partly due to some local adaptation.

Using proteomic data, we were able to tentatively identify the proteins contained in the bands that discriminate populations and therefore whose quantity might have evolved under heterogeneous selection. This identified several proteins known or suspected to be involved in virulence, but also proteins that have no identified putative function and would therefore deserve further investigation.

The combination of these two methods appears to be a powerful framework for population proteomics (Biron et al., 2006) and for the study of eco-evolutionary feedbacks between the proteomic level and the population and ecosystem levels (Diz et al., 2012). This is of interest not only to study evolution in the field at an intermediate level between the genome and phenotypes, or to understand the role of evolution in chemical ecology, but also for more applied issues. For example, biological control agents, sampled in the field in their area of origin, are often reared in the laboratory on a host species other than the target one (van Lenteren, 2003). The framework presented here could be used to monitor their potential unwanted adaptation to the laboratory host.

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

# REFERENCES

# 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 FST. HM-H, MP, DC, J-LG, ÉD, and JV performed the writing and editing of the manuscript. All authors read and approved the final manuscript.

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

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

(Hymenoptera: Figitidae) which parasitize Drosophila in the Rhone valley (S-E France). Ann. Soc. Entomol. Fr. 35, 97–103.

virulence? Physiol.  Sci. Rep. 8, 1–15. doi: 10.1038/s41598-018-33943-5

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

two populations of a Drosophila parasitoid. Biol. J. Linn. Soc. 117, 231–240. doi: 10.1111/bij.12644

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

Copyright © 2019 Mathé-Hubert, Kremmer, Colinet, Gatti, Van Baaren, Delava and Poirié. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Venom Costs and Optimization in Scorpions

### Edward R. J. Evans 1 \*, Tobin D. Northfield2,3, Norelle L. Daly 1 and David T. Wilson1 \*

*1 Centre for Molecular Therapeutics, Australian Institute of Tropical Health and Medicine, James Cook University, Cairns, QLD, Australia, 2 Department of Entomology, Tree Fruit Research and Extension Center, Washington State University, Wenatchee, WA, United States, 3 Centre for Tropical Environmental and Sustainability Science, College of Science and Engineering, James Cook University, Cairns, QLD, Australia*

Scorpions use venoms as weapons to improve prey capture and predator defense, and these benefits must be balanced against costs associated with its use. Venom costs involve direct energetic costs associated with the production and storage of toxins, and indirect fitness costs arising from reduced venom availability. In order to reduce these costs, scorpions optimize their venom use via evolutionary responses, phenotypic plasticity, and behavioral mechanisms. Over long timescales, evolutionary adaptation to environments with different selection pressures appears to have contributed to interspecific variation in venom composition and stinger morphology. Furthermore, plastic responses may allow scorpions to modify and optimize their venom composition as pressures change. Optimal venom use can vary when facing each prey item and potential predator encountered, and therefore scorpions display a range of behaviors to optimize their venom use to the particular situation. These behaviors include varying sting rates, employing dry stings, and further altering the volume and composition of venom injected. Whilst these cost-reducing mechanisms are recognized in scorpions, relatively little is understood about the factors that influence them. Here, we review evidence of the costs associated with venom use in scorpions and discuss the mechanisms that have evolved to minimize them.

Keywords: scorpion, venom, optimization, evolution, plasticity, behavior

# INTRODUCTION

Venomous organisms inject chemical cocktails into their predators and prey in order to disrupt normal biological functioning in their target (Fry et al., 2009; Casewell et al., 2013). These chemical weapons are often rich in proteins, peptides, and small molecules (Inceoglu et al., 2003; Escoubas et al., 2008; Calvete et al., 2009; Fry et al., 2009; Villar-Briones and Aird, 2018). Whilst venom provides survival benefits by aiding in prey capture and predator defense, the benefits come with costs. These costs are two-fold, involving direct energetic costs associated with production and storage of toxins (McCue, 2006; Nisani et al., 2007, 2012), and further indirect costs associated with a reduced capacity to capture prey or defend when supplies are depleted. Whilst these costs have different types of impacts on venomous animals, the methods to reduce these costs can overlap. It has been proposed that due to costs associated with venom use, organisms will meter/optimize the volume of venom they inject in order to use their venom as economically as possible in different situations (Hayes et al., 2002; Wigger et al., 2002; Hayes, 2008; Morgenstern and King, 2013). Research into optimal venom use has primarily focussed on snakes (Hayes, 1995, 2008; Hayes et al., 2002; Young et al., 2002), spiders (Wigger et al., 2002; Wullschleger and Nentwig, 2002; Hostettler and Nentwig, 2006; Nelsen et al., 2014; Cooper et al., 2015), and scorpions (Edmunds and Sibly, 2010; Nisani and Hayes, 2011; Lira et al., 2017), with the latter being the focus of this review.

### Edited by:

*Maria Vittoria Modica, Stazione Zoologica Anton Dohrn, Italy*

### Reviewed by:

*Timothy Jackson, The University of Melbourne, Australia Cassandra Marie Modahl, National University of Singapore, Singapore*

### \*Correspondence:

*Edward R. J. Evans edwardrobertjonathan.evans@ my.jcu.edu.au David T. Wilson<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Chemical Ecology, a section of the journal Frontiers in Ecology and Evolution*

> Received: *27 March 2019* Accepted: *13 May 2019* Published: *06 June 2019*

### Citation:

*Evans ERJ, Northfield TD, Daly NL and Wilson DT (2019) Venom Costs and Optimization in Scorpions. Front. Ecol. Evol. 7:196. doi: 10.3389/fevo.2019.00196*

Scorpions utilize venom to both capture prey and defend against predators, which can include organisms with very different physiologies and susceptibility to venom components (Gangur et al., 2017; van der Meijden et al., 2017). Over 400 million years of evolution (Dunlop and Selden, 2009) have led scorpions to develop a wide range of mechanisms that help minimize the costs of venom use. Successful prey capture and predator defense will ultimately affect a scorpion's evolutionary fitness, and therefore selection on venom composition and concentration is generally influenced by both prey and/or predators (Tian et al., 2008; Weinberger et al., 2009; Gangur et al., 2018). As selection pressures vary between environments, so will optimal investment in venom (Gangur et al., 2018). Scorpions adapted to different ecological niches often show large differences in venom composition (de la Vega et al., 2010) and stinger morphology (van der Meijden et al., 2013; van der Meijden and Kleinteich, 2017), and this likely reflects responses to different selection pressures.

# THE COSTS OF VENOM USE IN SCORPIONS

# Direct

Research on the energetic demands of venom use have focussed on costs of production, rather than maintenance, as it is difficult to measure the energy used in maintaining and storing toxins experimentally. Nisani et al. (2007) showed that after depleting the venom glands of scorpions (Parabuthus transvaalicus), metabolic rates increased by 39% during the first 3 days of regeneration. A later study (Nisani et al., 2012) found milked P. transvaalicus had on average a 21% higher metabolic rate than un-milked scorpions during the first 8 days of regeneration, but in this second study the rate did not rise as high during the first 3 days (Nisani et al., 2007). Metabolic rates fluctuated throughout the experiment, and the authors suggested this likely reflected the asynchronous regeneration of toxins (Nisani et al., 2012). Whilst differences were observed between the studies, both identify a large increase in metabolic rate above baseline levels, indicating that in scorpions venom production is an important energetic expense (Nisani et al., 2007, 2012). These studies also likely underestimate total energetic costs, as venom regeneration can take longer than 8 days to complete (Carcamo-Noriega et al., 2019). Furthermore, scorpion venom varies between species in complexity, toxins utilized, and volume stored and injected (de la Vega et al., 2010; Sunagar et al., 2013; van der Meijden et al., 2015); all of which alter energy requirements. Energetic costs are also dependent on the metabolic rate, which is likely to vary between species adapted to different ecological niches. Furthermore, as scorpions are ectotherms their metabolic rate will vary with environmental conditions, as Nime et al. (2013) showed scorpion activity is positively correlated with temperature.

# Indirect

Indirect costs are associated with the ecological limitations arising from depleted venom supplies, such as increased predation risk or reduced ability to capture prey. Scorpions can store a limited volume of venom (van der Meijden et al., 2015), and regeneration of toxins can take at least 2 weeks to be complete (Carcamo-Noriega et al., 2019), reducing the ability to use venom for prey capture and predator defense until venom supplies are restored. Ecological costs of venom depletion cannot easily be quantified, as they will vary widely between species and environments with fluctuating selection pressures. Nonetheless, evaluation of behavioral changes that arise when venom stores are depleted may serve as evidence of ecological costs. Such a behavioral response has yet to be reported in scorpions, but some spiders with depleted venom supplies will adapt their hunting behavior to target easily caught prey, as has been found in the wandering spider Cupiennius salei (Wullschleger and Nentwig, 2002). Scorpions and spiders often share many of their natural predators and prey items, therefore the ecological costs of venom depletion in both organisms may bear some similarity. Comparable behaviors in scorpions might include a switch to smaller prey that can be successfully captured with only the pedipalps while venom is replenishing.

# THE EVOLUTION OF OPTIMAL VENOM USE

Compared with other venomous taxa, scorpions are unusual in that they possess two main weapons: their stinging apparatus and their pedipalps. Species vary in their relative investment in these two weapons depending on their ecological niche, leading to the great morphological diversity in scorpion stinging apparatus and pedipalps seen between species (**Figure 1**) (van der Meijden et al., 2013). Burrowing species, such as members of the family Scorpionidae, are often sit-and-wait predators (Hadley and Williams, 1968; Bub and Bowerman, 1979; Shachak and Brand, 1983; Shivashankar, 1994), and possess large pedipalps that can be used to dig, grab passing prey, and block predators from entering their burrow (van der Meijden et al., 2010). Large pedipalps are often accompanied by a small stinging apparatus. Small tail size may be due to a trade-off in investment between pedipalps and stinging apparatus (van der Meijden et al., 2013), or a reduced tail may simply improve mobility in the confines of a burrow. Vagrant scorpion species, such as many members of the family Buthidae, have no permanent residence and forage more actively (Hadley and Williams, 1968). This group of scorpions generally rely more heavily on their sting to capture prey and defend themselves, and have evolved powerful and highly mobile tails (Warburg, 1998; Coelho et al., 2017). Many of these scorpions have developed potent venom, and the family Buthidae contains all species with medically significant stings (Santos et al., 2016). Buthids often possess small pedipalps, suggesting an evolutionary trade-off between pedipalps and stinging apparatus may be present (van der Meijden et al., 2013), and/or small pedipalps may improve mobility and energetic efficiency while actively foraging. The evolution of potent venom in buthids may further have reduced the advantages that large pedipalps provide, making them energetically unfavorable.

# BEHAVIORAL MECHANISMS TO OPTIMIZE VENOM USE

Optimal use of venom will vary between each interaction with prey and predators, influenced by the size and identity of the prey/predator. This has led scorpions to evolve a range of behavioral mechanisms allowing them to optimize their venom use when facing specific prey (Edmunds and Sibly, 2010) and under particular levels of threat (Nisani and Hayes, 2011).

# The "Decision" to Sting

To reduce the costs associated with unnecessary venom use, scorpions adapt their hunting strategies to the particular prey items targeted (Simone et al., 2018), and are less likely to use venom when capturing small or easily subdued prey (Rein, 1993; Edmunds and Sibly, 2010). Spiders appear to target the injection of their venom toward the thorax or head of prey items to maximize venom efficiency (Wigger et al., 2002; Carlson et al., 2014). However, to our knowledge it is not known if scorpions seek to apply stings to an optimal location, or if an optimal sting location even exists, as one study on Bothriurus bonariensis found that sting location did not affect the time taken to subdue prey (Simone et al., 2018).

Both the size and activity of prey items can influence a scorpion's choice to sting, as Parabuthus liosoma, Parabuthus pallidus and Hadrurus spadix sting larger and more active prey items more frequently (Rein, 1993; Edmunds and Sibly, 2010). Rein (1993) observed that the scorpions did not use their sting immediately when encountering prey, but would rather grab with their pedipalps and apply stings if the prey continued to struggle, presumably to minimize venom use whilst ensuring predation success. Ontogenetic changes in stinging behavior can also be used to optimize venom use. Older Paruroctonus boreus and Pandinus imperator use their larger pedipalps to overpower prey and sting less often, avoiding using venom (Cushing and Matherne, 1980; Casper, 1985).

In addition to trade-offs between the use of venom and pedipalps for prey capture, there may be trade-offs in venom use and mobility when avoiding predation, as faster scorpions appear less likely to sting predators (Carlson et al., 2014; Miller et al., 2016). For example, female Centuroides vittatus scorpions that are heavier and less mobile, are more likely to sting a potential predator than males, which are more likely to sprint to safety (Carlson et al., 2014; Miller et al., 2016). Furthermore, within each sex, sprint speed decreases and sting rate increases with mass, indicating that higher rates of aggression are associated with reduced mobility (Carlson et al., 2014; Miller et al., 2016). Through fleeing, the males and smaller scorpions are not only able to avoid being eaten, but they also save their venom supply for future encounters and do not need to expend energy regenerating toxins, thus reducing ecological and energetic venom costs.

Unlike other scorpions, seven species of Parabuthus spray venom defensively, which may cause irritation to the sensitive tissues, such as eyes, of predators (Newlands, 1974; Nisani and Hayes, 2015). P. transvaalicus defensively spray when presented with both air-flow and touch stimuli simultaneously, suggesting the behavior may be optimized toward high-threat scenarios where defensive tactics must be implemented before a predator gets close enough to sting (Nisani and Hayes, 2015). These scorpions are often attacked by predators such as grasshopper mice, which disarm scorpions by biting off their tails, and therefore place their face in close proximity to the telson (Nisani and Hayes, 2015). Compared with injection, sprayed venom has a higher risk of missing its target, likely increasing the costs of venom necessary to deter predators. However, these costs are likely offset by the advantage of deterring predators while the predator is still at a distance.

# "Dry" Stings

Scorpions can still avoid venom use when stinging in defensive encounters by employing "dry" stings, where no venom is injected. Scorpions readily utilize dry stings in defensive situations (Nisani and Hayes, 2011; Lira et al., 2017; Rasko et al., 2018). The factors that determine whether venom is injected in a sting, however, are currently unclear. In P. transvaalicus, dry stinging behavior is correlated with threat level, with the scorpions employing dry stings more frequently when the threat level is low (Nisani and Hayes, 2011). Additionally, this study suggests P. transvaalicus, when induced to sting multiple times in succession, use dry stings more frequently early in the stinging sequence and are more likely to inject venom as the threat persists (Nisani and Hayes, 2011). The combination of these results provides evidence that scorpions use dry stings in low threat situations to optimize their venom use. In contrast to the findings by Nisani and Hayes (2011) studies on dry stinging behavior in other scorpion species have produced differing results. Lira et al. (2017) presented Tityus stigmurus with the same "low-threat" and "high-threat" stimuli described by Nisani and Hayes (2011), and found no correlation between dry sting rate and threat level. Furthermore, it has been shown that repeated simulated attacks against Hadrurus arizonensis lead to an increase in dry sting rate, despite venom remaining in the gland (Rasko et al., 2018). This latter result seems counter-intuitive, as scorpions might be expected to increase their defensive investment as a threat persists. Whilst the studies investigating the factors influencing dry stinging behavior in scorpions are limited, the evidence supports the idea that at least some scorpion species utilize dry stings as a means to optimize their venom in defensive contexts (Nisani and Hayes, 2011), while others may not (Lira et al., 2017; Rasko et al., 2018). Further research should aim to identify whether interspecific differences are truly occurring, or if methodological differences between the studies are responsible for the observed differences. Furthermore, it is not currently known if scorpions utilize dry stings as a tactic to save venom when capturing prey, as spiders use dry bites (Malli et al., 1998; Wigger et al., 2002). Spiders, however, have the ability to masticate their prey with their fangs, but the dry sting of a scorpion provides comparatively little aid in the incapacitation of prey. It is therefore unlikely that scorpions use dry stings to save venom when targeting prey.

# Volume Injected

While scorpion stinging behavior involves a dichotomy between dry stings vs. stings with venom injected, scorpions also have the ability to vary the volume of venom they inject, both within each sting and through the application of multiple stings (Nisani and Hayes, 2011; van der Meijden et al., 2015). P. transvaalicus injected twice as much venom per single sting in highthreat situations compared with low-threat situations, indicating scorpions may use this tactic to vary their defensive investment in response to perceived threat level. The defensive sprays of P. transvaalicus display variable duration and flow rate suggesting the volume expelled could be controlled by contraction of the venom gland, but it is not currently known if the volume sprayed is influenced by threat level (Nisani and Hayes, 2015).

The volume injected in single stings may be limited by morphological constraints, or the time that the aculeus is pierced into the target (van der Meijden et al., 2015). When scorpions are faced with repeated attacks from predators, they will continue to defensively sting as the attacks continue. Experiments into the defensive investment of scorpions in response to predation threat suggest that scorpions will repeatedly sting predators as the threat persists (Lira et al., 2017), but it is unclear whether the investment per attack increases or decreases with sting number (Nisani and Hayes, 2011; Rasko et al., 2018). Targeting prey, scorpions often hold on with their pedipalps and judiciously apply stings as the prey continues to struggle (Casper, 1985; Rein, 1993). The rate of stings increases with both prey size and activity (Edmunds and Sibly, 2010), suggesting that scorpions are being frugal with their venom, and only apply extra stings as necessary.

# Composition Injected

In addition to reducing venom costs by metering the volume of venom they inject, scorpions are able to alter the composition of venom injected into their target and avoid unnecessarily injecting costly venom components. Scorpion venom is heterogeneous and changes in composition as it is expelled from the aculeus

(**Figure 2**). As the venom is secreted from the aculeus tip, the initial expulsion is a clear liquid, followed by an opalescent liquid, and finally turns milky colored and viscous (Yahel-Niv and Zlotkin, 1979). These different secretions also vary in toxicity (Yahel-Niv and Zlotkin, 1979). Inceoglu et al. (2003) found that the initial clear secretion in P. transvaalicus constituted around 5% of the total venom volume within the gland, and they termed this "prevenom". Prevenom and main milky venom have distinct compositions and modes of action in both invertebrate and vertebrate targets (Inceoglu et al., 2003). The different mode of actions of prevenom and main venom may act together to induce greater toxicity, in a similar way to the toxin cabals employed by cone snails which target different pathways simultaneously (Olivera et al., 1999; Inceoglu et al., 2003). P. transvaalicus prevenom was found to contain six-times less peptide and protein concentration compared to the main venom, but a 16-fold higher potassium (K+) salt concentration (Inceoglu et al., 2003). The authors suggested this extremely high K+ concentration causes large and rapid depolarization of nerves in the target, causing quick paralysis in insects and pain in vertebrates (Inceoglu et al., 2003). Prevenom not only contains a much lower peptide/protein concentration, but a comparatively simplistic composition (Inceoglu et al., 2003). It is therefore expected that prevenom is metabolically cheaper to produce than the main venom, as K+ salt likely requires less energy to be replenished than peptides requiring production and folding (Inceoglu et al., 2003). The relative costs of prevenom vs. main venom have not been calculated experimentally, but in P. transvaalicus prevenom components appear to regenerate quickly and at little metabolic cost compared with other toxins (Nisani et al., 2012). Prevenom may therefore have evolved as a mechanism to avoid injecting larger volumes of peptide rich mixtures, thereby minimizing both metabolic and ecological

costs of depletion, although this connection is difficult to test experimentally (Nisani et al., 2012).

By using prevenom first, scorpions can save their main peptide-rich venom for high-threat situations or when initially low-threat situations escalate. This hypothesis is supported by evidence that the composition of venom (prevenom vs. main venom) scorpions inject is context dependent. In lowthreat situations, scorpions are able to avoid injecting their metabolically "expensive" mixtures of toxins, by injecting only prevenom (Nisani and Hayes, 2011). P. transvaalicus inject their main venom more frequently in high-threat situations, and in later stings when induced to sting repeatedly at both low and high threat levels (Nisani and Hayes, 2011). Furthermore, in lowthreat situations, T. stigmurus injected prevenom in all trials, but when faced with the high-threat treatments most of the scorpions injected their main milky venom secretion (Lira et al., 2017). The use of prevenom in low threat situations not only minimizes metabolic costs but also reduces ecological costs, as prevenom appears to regenerate faster than the main venom components (Nisani et al., 2012).

# ADAPTIVE PLASTICITY

Recent evidence suggests scorpions can modify their venom composition in response to predator exposure (Gangur et al., 2017). Repeated periodical encounters with a surrogate vertebrate predator (a taxidermied mouse) over a 6 week period led Hormurus waigiensis to appear to produce a higher relative abundance of some vertebrate specific toxins used in defensive situations, and a lower relative abundance of certain toxins specific to their invertebrate prey (Gangur et al., 2017). This study provided the first evidence for adaptive plasticity in venom compositions, and suggested it has evolved as a mechanism to allow for the optimization of venom use (Gangur et al., 2017). Modification of venom composition in response to environmental pressures could allow scorpions to further optimize venom use in different environments. In environments with few predators, scorpions may not require large quantities of defensive toxins, but as predator abundance increases so does the need to defend themselves. Therefore, ability to plastically change venom composition can allow scorpions to prioritize their resources and minimize the costs of venom use. It is currently unclear what environmental cues (e.g., olfactory) led to the plastic response observed by Gangur et al. (2017). Furthermore, it is unclear if the response was targeted specifically at the presence of the mouse, or if it was a uniform response to increased predation pressure.

Unlike the response from simulated top-down predation pressure, Gangur et al. (2017) did not identify changes in venom composition in response to a scavenging vs. predacious diet, where venom is not required for prey capture. This may be due to the unpredictable nature of scavenging, and the potential need to kill prey in the future, regardless of current carrion availability. Alternatively, venom may need to be maintained for its defensive function. In contrast, it may also be that experimental conditions do not represent the bottom-up pressures experienced in the wild, as crickets may be more easily subdued than natural prey items. H. waigiensis is a burrowing species that has evolved large pedipalps and a small stinging apparatus, and further studies should evaluate whether more active species that rely more heavily on their sting to capture prey respond differently to a changing diet.

# CONCLUSIONS AND FUTURE DIRECTIONS

Scorpions experience direct costs associated with the production and storage of toxins, and indirect costs associated with impaired ecological function when their venom is depleted. Optimal venom use minimizes these costs, maximizing the survival benefit venom provides. On the broadest scales, optimal venom investment has contributed to the divergence of stinger morphology and venom compositions between species adapted to different environments (Tian et al., 2008; Sunagar et al., 2013; van der Meijden et al., 2013). Optimal venom use can be influenced by factors such as prey/predator identity, and scorpions therefore utilize a suite of behavioral tactics to minimize waste. These include varying sting frequency, employing dry stings, and further controlling the volume and composition of venom injected (Nisani and Hayes, 2011). Scorpions may also plastically adapt their venom composition (Gangur et al., 2017), allowing them to optimize venom use as selection pressures change. Whilst the presence of these mechanisms and behaviors are well-documented, the factors influencing them are poorly understood. Current knowledge of venom optimization has generally relied upon correlative research, where the selective forces driving the correlations are inferred, rather than directly measured. There is evidence of venom costs, benefits for prey capture and predator defense, and behavioral and trait phenotypes that appear to reduce these costs and maximize benefits. However, there is little direct evidence tying changes in phenotypes to changes in costs or benefits to describe a mechanistic link. Controlled selection experiments or phylogenetic studies that consider species interactions can help describe links between selection and evolutionary response in arms races (Pimentel, 1968; Kursar et al., 2009; Toju et al., 2011; Betts et al., 2018), and may help better describe how observed venom optimization mechanisms have evolved. Future work is needed to investigate whether observed changes are due to adaptive responses or physiological limitations, the extent that these mechanisms are influenced by the environment, and how widespread they are across different scorpion species.

# AUTHOR CONTRIBUTIONS

The manuscript was written by EE and edited by TN, DW, and ND.

# FUNDING

The Northcote Trust supports EE with a Northcote Graduate Scholarship.

### Evans et al. Venom Optimization in Scorpions

# REFERENCES

M. Hoggren, M. E. Douglas, H. W. Greene (Eagle Mountain: Eagle Mountain Publishing), 207–233.

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

Copyright © 2019 Evans, Northfield, Daly and Wilson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Multifunctional Toxins in Snake Venoms and Therapeutic Implications: From Pain to Hemorrhage and Necrosis

Camila R. Ferraz 1,2, Arif Arrahman3 , Chunfang Xie3 , Nicholas R. Casewell 4 , Richard J. Lewis 1 , Jeroen Kool 3 and Fernanda C. Cardoso1 \*

1 Centre for Pain Research, Institute for Molecular Bioscience, The University of Queensland, Brisbane, QLD, Australia, 2 Departamento de Ciências Patológicas, Centro de Ciências Biológicas, Universidade Estadual de Londrina, Londrina, Brazil, 3 Division of BioAnalytical Chemistry, Amsterdam Institute for Molecules Medicines and Systems, Vrije Universiteit Amsterdam, Amsterdam, Netherlands, 4 Centre for Snakebite Research & Interventions, Liverpool School of Tropical Medicine, Liverpool, United Kingdom

### Edited by:

Kartik Sunagar, Indian Institute of Science (IISc), India

### Reviewed by:

Kempaiah Kemparaju, University of Mysore, India Stephen P. Mackessy, University of Northern Colorado, United States

> \*Correspondence: Fernanda C. Cardoso<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Chemical Ecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 31 March 2019 Accepted: 24 May 2019 Published: 19 June 2019

### Citation:

Ferraz CR, Arrahman A, Xie C, Casewell NR, Lewis RJ, Kool J and Cardoso FC (2019) Multifunctional Toxins in Snake Venoms and Therapeutic Implications: From Pain to Hemorrhage and Necrosis. Front. Ecol. Evol. 7:218. doi: 10.3389/fevo.2019.00218 Animal venoms have evolved over millions of years for prey capture and defense from predators and rivals. Snake venoms, in particular, have evolved a wide diversity of peptides and proteins that induce harmful inflammatory and neurotoxic effects including severe pain and paralysis, hemotoxic effects, such as hemorrhage and coagulopathy, and cytotoxic/myotoxic effects, such as inflammation and necrosis. If untreated, many envenomings result in death or severe morbidity in humans and, despite advances in management, snakebite remains a major public health problem, particularly in developing countries. Consequently, the World Health Organization recently recognized snakebite as a neglected tropical disease that affects ∼2.7 million p.a. The major protein classes found in snake venoms are phospholipases, metalloproteases, serine proteases, and three-finger peptides. The mechanisms of action and pharmacological properties of many snake venom toxins have been elucidated, revealing a complex multifunctional cocktail that can act synergistically to rapidly immobilize prey and deter predators. However, despite these advances many snake toxins remain to be structurally and pharmacologically characterized. In this review, the multifunctional features of the peptides and proteins found in snake venoms, as well as their evolutionary histories, are discussed with the view to identifying novel modes of action and improving snakebite treatments.

Keywords: snake venoms, multifunctional toxins, pathological mechanisms, evolution, snakebite treatment, pain, hemotoxicity, myotoxicity

# INTRODUCTION

The composition and evolutionary histories of animal venoms have fascinated the scientific community for centuries. Venoms have evolved over millions of years to facilitate prey capture and/or defense from predators and rivals. Snake venoms, in particular, likely originated in the Cenozoic Era (Fry, 2005; Fry et al., 2006), and they are amongst the most well-characterized of animal venoms, comprising a complex mixture of toxic, and pharmacologically-active proteins and peptides (Casewell et al., 2013; Chan et al., 2016). Tragically, snake envenomation is a significant health and economic burden worldwide. It is estimated 1.8–2.7 million snakebites and 81,410– 137,880 deaths occur annually worldwide (Kasturiratne et al., 2008; Gutierrez et al., 2017), a problem which is mostly associated with agricultural work, especially in South and Southeast Asia, sub-Saharan Africa, and Central and South America (Kasturiratne et al., 2008; Harrison et al., 2019). In 2017, the World Health Organization (WHO) finally recognized the snakebite as a priority neglected tropical disease in which morbidity and mortality affects mostly individuals under 30 years old, who are often the most economically productive members of a community (WHO, 2018).

Snake venoms have a distinct complexity when compared to venoms from other animals such as spiders, scorpions, and cone snails (Zelanis and Tashima, 2014). In these animal venoms, the pharmacological effects are primarily caused by disulfide bridged peptides, whilst snake venoms consist of a more diverse array of larger proteins and peptides which results in a wider variety of pharmacological and toxicological effects (Zhang, 2015). These venoms comprise 50–200 components distributed in dominant and secondary families which can be presented in multiple proteins and peptides isoforms (Vonk et al., 2011; Slagboom et al., 2017; Tasoulis and Isbister, 2017). The dominant families are secreted phospholipases A2 (PLA2s), snake venom metalloproteinases (SVMP), snake venom serine proteases (SVSP), and three-finger peptides (3FTX), while the secondary families comprise cysteine-rich secretory proteins, Lamino acid oxidases, kunitz peptides, C-type lectins, disintegrins, and natriuretic peptides (Slagboom et al., 2017; Tasoulis and Isbister, 2017; Munawar et al., 2018). Interestingly, snake venom composition varies interspecifically (Fry et al., 2008; Tasoulis and Isbister, 2017), as well as intraspecifically, with many factors influencing this diversity including age (Dias et al., 2013), gender (Menezes et al., 2006; Zelanis et al., 2016), location (Durban et al., 2011; Goncalves-Machado et al., 2016), diet (Barlow et al., 2009), and season (Gubensek et al., 1974). This variability phenomenon underpins toxin diversity and multifunctionality, and is of great importance to be considered in antivenom production and envenomation treatment (Gutierrez et al., 2017).

The pharmacological effects of snake venoms are classified into three main types, hemotoxic, neurotoxic, and cytotoxic (WHO, 2010). The major toxins involved in these effects are the PLA2s, SVMPs, SVSPs, and 3FTXs, that alone or in combination, are responsible for the multiple pharmacological effects occurring in snakebite victims. For example, some PLA2s and 3FTX are able to act on pre- or post-synaptic junctions as antagonist of ion channels and nicotinic or muscarinic receptors to induce severe neurotoxicity such as paralysis and respiratory failure (Fry et al., 2003; Lynch, 2007; Casewell et al., 2013; Harris and Scott-Davey, 2013; Tsetlin, 2015). In addition, other PLA2s and 3FTXs, along with SVMPs, cause local tissue damage resulting in swelling, blistering, bruising, and necrosis, and systemic effects such as hypovolemic shock (Gutierrez and Rucavado, 2000; Gutierrez et al., 2005; Harris and Scott-Davey, 2013; Rivel et al., 2016). Furthermore, SVSPs and SVMPs induce hemostatic and cardiovascular effects as coagulopathy, hypotension and hemorrhage (Slagboom et al., 2017). Interestingly, some PLA2s, SVSPs, and SVMPs are also capable of triggering severe pain by modulating pain pathways through activation of ion channels, such as transient receptor potential vanilloid type 1 (TRPV1) and acid-sensing ion channel (ASIC) (Bohlen et al., 2011; Zhang et al., 2017); and/or by pain sensitization through inflammatory mediators (Zychar et al., 2010; Menaldo et al., 2013; Ferraz et al., 2015; Mamede et al., 2016; Zambelli et al., 2017b). The inflammation induced by the elapid and viper venoms is widely reported to produce pain or hyperalgesia in human and in experimental models (Hifumi et al., 2015; Bucaretchi et al., 2016; Kleggetveit et al., 2016; Mamede et al., 2016). Unfortunately, these are not completely reversed by antivenom and anti-inflammatory therapies (Picolo et al., 2002; Ferraz et al., 2015; Hifumi et al., 2015).

The toxicological effects induced by snakebite are currently treated with intravenous administration of antivenom in combination to analgesics, fluid therapy, hemodialysis and/or antibiotics (Gutierrez et al., 2017). Although sufficient in most cases, snakebite treatments have been challenged by the continuous high numbers of clinical illness and mortality associated with snakebites worldwide (WHO, 2018). Furthermore, chronic morbidity following snakebites have been underestimated, with many victims reporting chronic symptoms in the bitten region, including complex regional pain syndrome (CRPA) (Seo et al., 2014; Kleggetveit et al., 2016) and musculoskeletal disabilities (Jayawardana et al., 2016).  Therefore, current research efforts are directed to the development of more effective snakebite therapies able to generically fully inhibit the major toxic components of snake venoms in order to better overcome severe acute and chronic effects caused by snakebite.

In light of the public health importance and the complexity of snake venoms, in this review we highlight the multifunctionality, structure-activity relationships and evolution of proteins and peptides in snake venoms. We aim to provide a better understanding of their action mechanisms and effects, and to bring attention to their undetermined targets and a host of potential novel therapeutic targets that might have implications for improving the treatments of snakebites.

# PHOSPHOLIPASES (PLA2s)

Phospholipases A2 play an important role in the neurotoxic and myotoxic effects of snakebites (Harris and Scott-Davey, 2013). These proteins have molecular masses of 13–15 kDa and are classified into groups I and II, which are found as major components in the venoms of Elapidae and Viperidae, respectively (Six and Dennis, 2000; Harris and Scott-Davey, 2013). In addition, a third group of PLA2s, termed IIE, have been predominately recovered from the venoms of nonfront fanged snakes, although their importance in the venom arsenal remains unclear (Fry et al., 2012; Perry et al., 2018). Studies reconstructing the evolutionary history of this multilocus gene family have demonstrated that each of these PLA2s types (I, II, and IIE) have been independently recruited into snake venom systems (Fry et al., 2012; Junqueira-De-Azevedo et al., 2015), suggesting they have evolved their toxic properties by convergent evolution (Lomonte and Rangel, 2012). Although PLA2s from vipers and elapids share similar enzymatic properties, both types have undergone extensive gene duplication over evolutionary time, seemingly facilitating the evolution of new toxic functions, and resulting in different patterns of residue conservation (Lynch, 2007; Vonk et al., 2013; Dowell et al., 2016; **Figure 1A**). In addition, the venoms of viper snakes contain isoforms of group II PLA2s that are catalyticallyactive (e.g., Asp49) and catalytically-less active (e.g., Lys49) (Lomonte and Rangel, 2012).

A number of PLA2s exert strong myotoxic effects which often lead to severe necrosis (Harris and Maltin, 1982; Gutierrez and Ownby, 2003), and many of these toxins also promote inflammation, including edema formation, cytokine production and leukocyte recruitment, pain by inducing thermal allodynia and mechanical hyperalgesia, paralysis through block of neuromuscular transmission and intensify hemorrhage by inhibiting coagulation (**Table 1**) (Camara et al., 2003; Chacur et al., 2003; Camargo et al., 2008; Teixeira et al., 2011; Lomonte and Rangel, 2012; Harris and Scott-Davey, 2013; Casais-E-Silva et al., 2016; Costa et al., 2017; Zambelli et al., 2017b; Zhang et al., 2017). Neurotoxic effects caused by these toxins, as well as some of their proinflammatory effects, occurs via the modulation of pre-synaptic terminals as well as sensory nerveendings (Camara et al., 2003; Harris and Scott-Davey, 2013; Sribar et al., 2014; Zhang et al., 2017). The PLA2s pre-synaptic effects are characteristic of β-neurotoxins and target the motor nerve terminals at the neuromuscular junction (Sribar et al., 2014; Gutierrez et al., 2017). The mechanisms in which certain PLA2s exert their pre-synaptic effects are not fully understood, and the primary targets remain unidentified, although the PLA2 β-bungarotoxin is known to bind to K+ channels at the pre-synaptic terminals via an accessory kunitz subunit (Benishin, 1990; Sribar et al., 2014; **Figure 1B**). Overall, these presynaptic effects induce robust exocytosis of the neurotransmitters vesicles reserves which consequently lead to the depletion of neurotransmitter release in the neuromuscular junction to promote muscle paralysis (Harris et al., 2000; Sribar et al., 2014; Gutierrez et al., 2017).

The inflammation induced by PLA2s has non-neurogenic and neurogenic (substance-P dependent) components (Camara et al., 2003; Camargo et al., 2008; Costa et al., 2017; Zhang et al., 2017). The non-neurogenic component is mostly caused by the hydrolysis of membrane lipids that generate potent pro-inflammatory lipid mediators (Costa et al., 2017). Additional non-neurogenic and neurogenic inflammations induced by PLA2s use more complex mechanisms still not fully understood. For example, leukocyte recruitment (De Castro et al., 2000), mastocytes degranulation (Menaldo et al., 2017), and macrophage activation (Triggiani et al., 2000; Granata et al., 2006; Giannotti et al., 2013) were demonstrated to occur independently of the generic PLA2s lipid hydrolysis catalytic activity. Furthermore, substance-P mediated neurogenic inflammation has been described to be induced by PLA2s from Crotalus durissus cascavella (Camara et al., 2003) and from Naja mossambica (Camargo et al., 2008). Interestingly, the C-terminal of Myotoxin-II (a Lys49-PLA2) isolated from Bothrops asper was able to activate macrophages, showing this region maybe be crucial for the observed enzymatic-independent inflammation (Giannotti et al., 2013) (**Figure 1A**).

The pain induced by PLA2s is driven by inflammatory processes and sensory neuronal activation. Bradykinin is an important mediator of the inflammatory pain induced by PLA2s (Moreira et al., 2014; Urs et al., 2014; Mamede et al., 2016; Zambelli et al., 2017b). It induces mechanical hyperalgesia dependent on the production of TNF-α, IL-1β, and prostaglandins (Cunha et al., 1992). This suggests that PLA2s contribute to an increase in arachidonic acid release from cell membranes and its availability to be processed by cyclooxygenase resulting in prostaglandin production (Verri et al., 2006). Corroborating this hypothesis, studies performed in rodents have demonstrated that PLA2s isolated from different snake venoms induced hyperalgesia mediated by biogenic amines, cytokines, prostaglandins, sympathomimetic amines, ATP, K+ release, purinergic receptor activation, and glial cell activation (Nunez et al., 2001; Chacur et al., 2003, 2004; Zhang et al., 2017). Interestingly, the presence of a strong catalytic activity in the PLA2s is not essential for its nociceptive activity as observed by the nociceptive effects of PLA2s "Lys49" variants (Lomonte et al., 1994; Rong et al., 2016; Zhang et al., 2017). Direct activation of sensory neurons was demonstrated by MitTx from Micrurus tener tener, a heteromeric complex between a PLA2 and a kunitz peptide (Bohlen et al., 2011), and by the Lys49 PLA2 BomoTx from the Brazilian viper Bothrops moojeni (Zhang et al., 2017). MitTx activates somatosensory neurons and was found to be a potent and selective agonist of ASIC channels (**Figure 1C**). This agonistic effect induces robust pain behavior in mice via activation of ASIC1 channels on capsaicin-sensitive nerve fibers (Bohlen et al., 2011). BomoTx also activated a cohort of sensory neurons to induce ATP release followed by activation of purinergic receptors (Zhang et al., 2017). Unfortunately, the primary target of this neuronal activation is still unknown. This same toxin induced non-neurogenic inflammatory pain, thermal hyperalgesia dependent on TRPV1 channels-expressing nerve fibers, and mechanical allodynia dependent on P2X2/P2X3 expressing fibers (Zhang et al., 2017).

The multifunctionality of PLA2s is evidenced by their myotoxic, neurotoxic and enzymatic functions, as well as by their inflammatory properties. There is evidence that separate domains and regions of the PLA2s structure participate in these various activities (**Figures 1A,B**). For example, for the Lys49-PLA2 from Bothrops asper and Agkistrodon piscivorus piscivorus, the Cterminal region of these toxins (residues 115–129) were identified as the active sites responsible for their myotoxic effects (Lomonte et al., 1994; Nunez et al., 2001) (**Figure 1A**). Interestingly, the same C-terminal region in BpirPLA2-I isolated from Bothrops

pirajai had anticoagulant activity through inhibition of platelet aggregation (Teixeira et al., 2011). Crotoxin B, an Asp49-PLA2, and a major component of the venom of Crotalus durissus terrificus, has toxic active sites fully independent of its enzymatic activity (Soares et al., 2001; **Figure 1D**), while Lemnitoxin, a basic class I PLA2 isolated from Micrurus lemniscatus, has strong myotoxic and proinflammatory effects but no neurotoxic activity (Casais-E-Silva et al., 2016). A detailed mutational study using the PLA2 OS2 from the Australian Taipan snake Oxyuranus scutellatus scutellatus revealed that a 500-fold loss in enzymatic activity had only a minor effect on its neurotoxicity (Rouault et al., 2006). Furthermore, the enzymatic activity of OS2 was dependent of the N- and C-terminal regions, and the Nterminal region had a major role in the central nervous system

and in the nociceptive and/or endematogenic properties in BThrTx1 (Zambelli et al., 2017a). The central histidine in the catalytic site of Crotoxin B is highlighted in red.

neurotoxicity. An alanine scan of the Lys49-PLA2 from Bothrops jararacussu (BThTx-I) demonstrated distinct regions involved in the hyperalgesia and edema (Zambelli et al., 2017a). In this study, the mutant Arg118Ala lost both nociceptive and edematogenic properties, Lys115Ala and Lys116Ala lost the nociceptive effects without interfering with the edema formation and Lys122Ala lost the nociceptive properties and had weak inflammatory effects (**Figure 1E**). Similarly, an independent study showed the Lys122Ala substitution led to reduced membrane damaging and myotoxic activities (Ward et al., 2002). This C-terminal region is characterized by the presence of basic and hydrophobic residues which have been strongly associated with the ability of PLA2s to interact and penetrate the lipid bilayer (Delatorre et al., 2011; Gutierrez and Lomonte, 2013).

The variability of PLA2 isoforms observed in the venom of different snakes may reflect a variety of factors, including the evolutionary history, phylogeography, diet, and/or environmental conditions relating to a species or populations within a species (Zancolli et al., 2019). Many snake venom toxins are known to be encoded by multi-locus gene families (Casewell et al., 2013), and the resulting genes within those families have been found clustered together in arrays on microchromosomes, likely as the result of tandem gene duplication events (Vonk et al., 2013). The process of gene duplication and loss underpins the evolution of many snake venom toxin families, including the PLA2s (Lynch, 2007; Vonk et al., 2013; Casewell et al., 2014; Dowell et al., 2016), with duplications likely initially stimulating a gene dosing effect while also freeing duplicates from evolutionary constraints, and thus enabling a scenario that may facilitate protein sub- and/or neo-functionalization (Lynch and Conery, 2000). Indeed, studies have demonstrated that extremely divergent venom phenotypes (e.g., neurotoxic vs. haemorrhagic) observed within populations of the same snake species, or between closely related species, are at least partially the result of major genomic differences in PLA2 toxin loci, with variation at different gene complexes resulting in markedly different haplotypes (Dowell et al., 2018; Zancolli et al., 2019). It remains unclear as to the specific processes that underpin such diversity, although natural selection driven by environmental factors and hybridization events have both been proposed (Dowell et al., 2018).

# SNAKE VENOM METALLOPROTEINASES (SVMPs)

Snake Venom MetalloProteinases (SVMPs) are zinc-dependent proteinases ranging from 20 to 110 kDa in size and are categorized into P-I, P-II, and P-III classes according to their structural domains (Hite et al., 1994; Jia et al., 1996; Fox and Serrano, 2005) (**Figure 2**). These toxins are major components of viper venoms and play a key role in the toxicity of these snake venoms (**Table 1**; Tasoulis and Isbister, 2017). Venom SVMPs have evolved from ADAM (a disintegrin and metalloproteinase) proteins, specifically ADAM28 (Casewell, 2012), with the P-III being the most basal structural variant consisting of metalloproteinase, disintegrin-like, and cysteine-rich domains (Moura-Da-Silva et al., 1996; Casewell et al., 2011; **Figure 2B**). Subsequently, P-II SVMPs diverged from P-IIIs and consist of a metalloproteinase and disintegrin domain, with the latter typically detected in venom as a proteolytically processed product (Fox and Serrano, 2008; Casewell et al., 2011). The final class, P-I SVMPs which consist only of the metalloproteinase domain, appeared to have evolved on multiple independent occasions in specific lineages as a result of loss of the P-II disintegrinencoding domain (Casewell et al., 2011). Throughout this diverse evolutionary history, SVMPs show evidence of extensive gene duplication events, coupled with bursts of accelerated molecular evolution (Casewell et al., 2011; Vonk et al., 2013). However, while all three classes of SVMPs have been described from viper venoms, only P-III SVMPs have been detected in elapid and "colubrid" venoms (Casewell et al., 2015). While these P-III SVMPs are typically relatively lowly abundant venom components in elapid snakes (e.g., <10% of venom toxins), they can be major components in "colubrids" (Mackessy and Saviola, 2016; Pla et al., 2017; Tasoulis and Isbister, 2017; Modahl et al., 2018a). These abundance differences likely underpin the distinct pathologies observed following envenomings by snakes found in these families. SVMPs contribute extensively to the hemorrhagic and coagulopathic venom activities following bites by viperid snakes, and the diversity of SVMPs isoforms often present in their venom likely facilitate synergistic effects, such as simultaneous action on multiple steps of the blood clotting cascade (Kini and Koh, 2016; Slagboom et al., 2017). Certain "colubrid" snakes, whether medically important or not, also show evidence of having multiple SVMP toxins in their venom, even if they do not show the sub-class diversity observed in the vipers (Mackessy and Saviola, 2016; Pla et al., 2017; Modahl et al., 2018a; Perry et al., 2018). However, it is relatively uncommon for elapid snakebites to cause systemic hemotoxicity (Slagboom et al., 2017) and this is likely a consequence of those venoms exhibiting little diversity or abundance of SVMPs, and instead usually being dominated by neurotoxic toxin families such as the 3FTXs and PLA2s (Tasoulis and Isbister, 2017).

As mentioned above, the SVMPs are known for their hemorrhagic activity as well as for their ability to influence multiple steps of the blood clotting cascade, resulting in a lethal combination of systemic hemorrhage and incoagulable blood in prey and/or victims (Markland and Swenson, 2013). Research has revealed that the effects of SVMP-induced hemorrhage relies on a mechanism that occurs in two steps (Gutierrez et al., 2005; Escalante et al., 2011). First, SVMPs cleave the basement membrane and adhesion proteins of endothelial cells-matrix complex to weaken the capillary vessels. During the second stage, the endothelial cells detach from the basement membrane and become extremely thin, resulting in disruption of the capillary walls and effusion of blood from the fragile capillary walls. In addition to the proteinase activity, SVMPs impact on homeostasis by altering coagulation, which contributes to their toxic hemorrhagic effects (Markland, 1998; Takeda et al., 2012; Slagboom et al., 2017). This occurs through modulation of factors such as fibrinogenase and fibrolase that mediate the coagulation cascade, depletion of pro-coagulation

Ferraz et al.

Their functional properties associated to pain, inflammation, hemorrhage, necrosis and paralysis are described.

the disintegrin-like domain is highlighted in green and the cysteine-rich domain is highlighted in blue. (B) Cartoon representation of the three-dimensional structure of the class P-III metalloproteinase VAPB2 from Crotalus atrox (PDB 2DW0). The metalloproteinase domain is colored in orange, the disintegrin-like domain (D-like) is colored in green and the cysteine-rich domain (Cys-rich) is colored in blue. The disulphide bridges are colored in yellow.

factors through consumption processes (e.g., Factor X, prothrombin and fibrinogen), platelet aggregation inhibition and inflammatory activities (Kamiguti, 2005; Kini, 2005; Takeda et al., 2012; Kini and Koh, 2016; Slagboom et al., 2017; Ainsworth et al., 2018).

Some SVMPs also induce inflammation, including edema, and pain by triggering hyperalgesia (Dale et al., 2004; da Silva et al., 2012; Bernardes et al., 2015). The inflammation induced by jararhagin, a class P-IIIb metalloproteinase with potent hemorrhagic and dermonecrotic activity isolated from Bothrops jararaca, produced TNF-α and IL-1β in vivo (Laing et al., 2003; Clissa et al., 2006), while BJ-PI2, isolated from the same snake species, lacked hemorrhagic and necrotic activities but still induced vascular permeability and inflammatory cell migration (da Silva et al., 2012). Curiously, the edema formation induced by jararhagin was independent of pro-inflammatory mediators such as TNF, IL-1β, and IL-6 (Laing et al., 2003). Neurogenic inflammation was also implicated in the local hemorrhage induced by Bothrops jararaca which was shown to be dependent on serotonin and other neuronal factors (Goncalves and Mariano, 2000). The mechanisms on how neurogenic inflammation is triggered by the snake venom components and how it participates in the hemorrhagic process are still not understood. Pain induced by SVMPs is characterized by hyperalgesia and inflammatory pain, which is dependent on the production of cytokines, nitric oxide, prostaglandins, histamine, leukotrienes, and migration of leukocytes, mast cell degranulation and NFkB activation (Fernandes et al., 2007; Bernardes et al., 2015; De Toni et al., 2015; Ferraz et al., 2015). However, the mechanisms underlying SVMP-induced pain are still poorly understood, with neurogenic inflammation and neuronal excitatory properties still underexplored. The research reported to date indicates that SVMPs induce the production of inflammatory mediators and activate cytokine/chemokine cascades and the release of prostaglandins and sympathetic amines to engender nociceptor sensitization (Verri et al., 2006; De Toni et al., 2015; Ferraz et al., 2015).

The multifunctional properties of SMVPs are also welldescribed. Class P-III SVMPs tend to display stronger hemorrhagic activity compared to P-I and P-II SVMPs, possibly due to the disintegrin-like and cysteine-rich domains enabling binding to relevant targets in the extracellular matrix of capillary vessels. The functions of these domains have been investigated in inflammation, revealing that these domains are sufficient to induce pro-inflammatory responses through production of TNF-α, IL-1β and IL-6, and leukocyte migration, in which mechanisms and primary targets are still unknown (Clissa et al., 2006; Ferreira et al., 2018). These observations suggest that these domains are involved in the inflammatory hyperalgesia induced by SVMPs. Furthermore, the pronounced hemorrhagic and necrotic activities are strongly dependent on biological effects driven by the disintegrin-like and cysteine-rich domains, as observed for BJ-PI2 (da Silva et al., 2012). The hemorrhagic activity of Bothrops jararaca venom was also shown dependent on neurogenic inflammation (Goncalves and Mariano, 2000). These findings implicate the disintegrin-like and/or cysteine-rich domains as key player(s) in these neurogenic mechanisms, and possibly in non-inflammatory pain, from which the neuronal targets are still to be identified.

# SNAKE VENOM SERINE PROTEINASES (SVSPs)

Snake Venom Serine Proteinases (SVSPs) belong to the S1 family of serine proteinases and display molecular masses ranging from 26 to 67 kDa and two distinct structural domains (**Figure 3**). These venom toxins have evolved from kallikreinlike serine proteases and, following their recruitment for use in the venom gland, have undergone gene duplication events giving rise to multiple isoforms (Fry et al., 2008; Vaiyapuri et al., 2012). SVSPs catalyze the cleavage of polypeptide chains on the C-terminal side of positively charged or hydrophobic amino acid residues (Page and Di Cera, 2008; Serrano, 2013). Similar to SVMPs, SVSPs have been described in the venom of a wide variety of snake families, although they are typically only abundant in viper venoms, and much less common in the venoms of elapid and "colubrid" snakes (Tasoulis and Isbister, 2017; Modahl et al., 2018a). Whilst the SVMPs are well-known for their ability to rupture capillary vessels, SVSPs execute their primary toxicity by altering the hemostatic system of their victims, and by inducing edema and hyperalgesia through mechanisms still poorly understood (**Table 1**). Hemotoxic effects caused by SVSPs include perturbations of blood coagulation (pro-coagulant or anti-coagulant), fibrinolysis, platelet aggregation and blood pressure, with potential deadly consequences for snakebite victims (Murakami and Arni, 2005; Kang et al., 2011; Serrano, 2013; Slagboom et al., 2017).

Pro-coagulant SVSPs have been described to activate multiple coagulation factors, including prothrombin and factors V, VII, and X (Kini, 2005; Serrano, 2013). For example, the activation of prothrombin produces thrombin which in turn produces fibrin polymers that are cross-linked. Thrombin also activates aggregation of platelets which, together with the formation of fibrin clots, results in coagulation (Murakami and Arni, 2005). In addition, platelet-aggregating SVSPs will activate the platelet-receptors to promote binding to fibrinogen and clot formation (Yip et al., 2005). These procoagulant and plateletaggregating activities will lead to the rapid consumption of key factors in the coagulation cascade and clot formation. On the other hand, anticoagulant SVSPs effects involve the activation of Protein C, which subsequently inactivates the coagulant factors Va and VIIIa (Kini, 2006). Furthermore, fibrinolytic SVSPs play an important role in the elimination of blood clots by acting as thrombin-like enzymes or plasminogen activators, which eliminates the fibrin in the clots and contributes significantly to the establishment of the coagulopathy (Kang et al., 2011; Serrano, 2013). Through the activation/depletion and inactivation of these coagulation factors, the clotting of blood is prevented, leading to uncoagulable blood, and external and internal bleeding.

Little is known about inflammatory responses and hyperalgesia induced by SVSPs. Studies suggest SVMPs and PLA2s have a pivotal role in the inflammatory responses and pain induced by snake venoms, while SVSPs have an important role in inflammation and a minor role in pain (Zychar et al., 2010; Menaldo et al., 2013). SVSPs in the venoms of Bothrops jararaca and Bothrops pirajai induce inflammation through edema formation, leucocyte migration (mainly neutrophils) and mild mechanical hyperalgesia, however, the mediators involved in these effects are still unknown (Zychar et al., 2010; Menaldo et al., 2013).

# THREE-FINGER TOXINS (3FTXs)

Three-fingers toxins (3FTXs) are non-enzymatic neurotoxins ranging from 58 to 81 residues that contain a three-finger fold structure stabilized by disulfide bridges (Osipov and Utki, 2015; Kessler et al., 2017; **Figure 4A**). They are present mostly in the venoms of elapid and colubrid snakes, and exert their neurotoxic effects by binding postsynaptically at the neuromuscular junctions to induce flaccid paralysis in snakebite victims (Barber et al., 2013). Three-finger toxins differ in length, with short-chain 3FTXs including α-neurotoxins, β-cardiotoxins, cytotoxins, fasciculins and mambalgins, which comprise of 57– 62 residues and 4 disulfide bridges, and long-chain 3FTXs which include α-neurotoxins and γ-neurotoxins, hannalgesin and κneurotoxins, and comprise 66–74 residues and five disulfide bridges. Furthermore, they can exist as monomers and as covalent or non-covalent homo or heterodimers. The diversity of 3FTX isoforms described above are a direct result of a diverse evolutionary history, whereby ancestral 3FTXs have diversified by frequent gene duplication and accelerated rates of molecular evolution. These processes, which are broadly similar to those underpinning the evolution of the other toxin families described above, are particularly associated with the evolution of a high-pressure hollow-fanged venom delivery system observed in elapid snakes (Sunagar et al., 2013). For example, gene duplication events have resulted in the expansion of 3FTX loci from one in non-venomous snakes like pythons, to 19 in the elapid Ophiophagus hannah (king cobra) (Vonk et al., 2013), and selection appears to have acted extensively on surface exposed amino acid residues in these resulting paralogous elapid 3FTX genes (Sunagar et al., 2013). The consequences of this evolutionary history are the differential production of numerous 3FTX isoforms that often exhibit considerable structural differences and distinct biological functions (**Figures 4B–E**). Although many elapid snakes exhibit broad diversity of these functionally varied toxins in their venom (e.g., multiple short- and long-chain 3FTX isoforms), it remains unclear why particular functional variants are enriched in the venoms of certain elapid lineages, such as the cytotoxin-rich venoms of cobras (genus Naja) or the neurotoxin-rich venoms of mambas (genus Dendroaspis) (Tan et al., 2015; Lauridsen et al., 2017; Ainsworth et al., 2018). However, evidence of taxon-specific 3FTXs in the venoms of certain "colubrid" snakes (Pawlak et al., 2006; Mackessy and Saviola, 2016; Modahl et al., 2018a), coupled with pseudogenization of 3FTXs in species that no longer rely on their venom (Li et al., 2005), suggests that

cytotoxin 1 from Naja atra (Uniprot P60304), and calliotoxin from Calliophis bivirgatus (Uniprot P0DL82). (B–E) Cartoon representation of the three-dimensional structure of the 3FTXs short-chain mambalgin-1 from Dendroaspis polylepis polylepis (PDB 5DU1) (B), long-chain α-bungarotoxin from Bungarus multicinctus (PDB 1ABT) (C), non-covalent homodimer α-cobratoxin from Naja kaouthia (PDB 4AEA) (D), and covalent heterodimer irditoxin from Boiga irregularis (PDB 2H7Z) (E). (F) Cartoon representation of the three-dimensional structure of Fascilulin-2 bound to the AChE (PDB 4BDT). The residues Arg24, Lys25, Pro31 and Leu35 which form hydrogen bonds with AChE are shown in orange. (G–J) Cartoon representation of the three-dimensional structure of the muscarinotoxin 1 (MT1, PDB 4DO8) (G) and muscarinotoxin 7 (MT7, PDB 2VLW) (H), and respective analogs displaying the modified loop 1 (PDB 3FEV) (I) and loop 3 (PDB 3NEQ) (J) in light orange color. (K) Neurotoxin II from N. oxiana (PDB 1NOR). The residues Ser29, His31, Gly33 and Thr34 which form hydrogen bonds with the α-subunit of nAChR are shown in orange, and the residue Arg32 forming ionic interactions with the α-subunit of nAChR is shown in yellow. (L) Neurotoxin b (NTb) from O. Hannah (PDB 1TXA). The residues Lys24, Trp26, and Asp28 that form hydrogen bonds with the α-subunit nAChR are shown in orange.

selection for prey capture may be at least partially responsible for influencing differential 3FTX representation.

Despite the shared three-finger fold, the 3FTXs have diverse targets and biological activities.  Their toxic biological effects include flaccid or spastic paralysis due to the inhibition of AChE and ACh receptors (Grant and Chiappinelli, 1985; Changeux, 1990; Marchot et al., 1998; Marquer et al., 2011), and activation of NaV1.4 (Yang et al., 2016) and L-type calcium channels (Garcia et al., 2001) in the periphery, necrosis through the action of cardiotoxins (cytotoxins) (Konshina et al., 2017), alteration of the cardiac rate through modulation of αand β-adrenoreceptors (Rajagopalan et al., 2007; Koivula et al., 2010), and altered homeostasis through inhibition of platelet aggregation (McDowell et al., 1992) and Factor X (Girish and Kini, 2016). In addition to their multitude of bio-activities, 3FTXs can remarkably display toxicities that target distinct classes of organisms as demonstrated in non-front fanged snake venoms that produce 3FTX isoforms which are non-toxic to mice but highly toxic to lizards, and vice-versa (Modahl et al., 2018b).

Some 3FTXs are able to induce analgesia through inhibition of ASIC channels (Salinas et al., 2014), while no 3FTXs are known to be involved in inflammation and hyperalgesia as commonly reported for other snake toxins. Furthermore, 3FTXs are relatively small compared to the other snake toxins discussed herein, and do not exhibit multiple domains to produce their multiple toxic functions. Nevertheless, the number of receptors, ion channels, and enzymes targeted by snake 3FTXs highlights the unique capacity of this fold to modulate diverse biological functions and the arsenal of toxic effects that are induced by 3FTXs. The unique multifunctionality of the 3FTX scafold occurs because of their resistance to degradation and tolerance to mutations and large deletions (Kini and Doley, 2010). Therefore, the structure-activity relationship of the 3FTXs is complex and yet to be fully understood. Their functional sites are located on various segments of the molecule surface. Conserved regions determine structural integrity and correct folding of 3FTXs to form the three loops, including eight conserved cysteine residues found in the core region. Aromatic residues (Tyr25 or Phe27) are conserved in most 3FTXs and essential to their folding. Another conserved features are the antiparallel β-sheet structure and charged amino acid residues also essential to stabilize the native conformation of the protein by forming hydrogen and ionic interactions, respectively (Torres et al., 2001). Additional disulfide bonds can be observed either in the loop I or loop II which can potentially change the activity of the 3FTX in some cases.

Specific amino acid residues in critical segments of the 3FTXs have been identified to be important for binding to their targets. For example, the interactions between fasciculin and AChE enzyme has been studied. The first loop or finger of fasciculin reaches down the outer surface of the enzyme, while the second loop inserts into the active site and exhibit hydrogen bonds and hydrophobic interaction (Harel et al., 1995; **Figure 4F**). Several basic residues in fasciculin make key contacts with AChE. From docking studies, hydrogen bonds, and hydrophobic interactions where shown to establish receptor-toxin assembly. Six amino acid residues (Lys25, Arg24, Asn47, Pro31, Leu35 and Ala12) of fasciculin interact with the AChE residues by forming hydrogen bonds at its active site. Hydrophobic interactions are also observed between eight amino acid residues (Lys32, Cys59, Val34, Leu48, Ser26, Gly36, Thr15, Asn20) from fasciculin and the enzyme active site (Waqar and Batool, 2015). These interactions involve charged residues but lacks intermolecular salt linkages.

Muscarinic toxins from mamba venoms, such as MT1 and MT7 (**Figures 4G,H**), act as highly potent and selective antagonists of M1 receptor subtype through allosteric interactions with the M1 receptor. Fruchart-Gaillard et al. (2012) synthesized seven chimeric 3FTXs based on MT1 and MT7 proteins that have remarkable affinity for α1A-adrenoceptor receptor subtypes but low affinity for M1 receptor. In this study, substitution within loop 1 and loop 3 weaken the toxin interactions with the M1 receptor, resulting in a 2-fold decrease in affinity (**Figures 4I,J**). Furthermore, modifications in loop 2 of the MT1 and MT7 significantly reduce the affinity for the M1 receptor. Interestingly, a significant increase in affinity was achieved on the α1A-adrenoceptor by combined modifications in loops 1 and 3, where loop 1 forms a critical interaction with the receptor (Fruchart-Gaillard et al., 2012). Another muscarinic toxin named MTβ was designed based on ρ-Da1a protein from Black Mamba which is known to have affinity for the muscarinic receptor. The Ser/Ile38 and Ala/Val43 substitutions in ρ-Da1a could be responsible for the increased affinity of MTβ for α1Badrenoreceptor and α1D-adrenoreceptor (Blanchet et al., 2013). These two residues were not located at the tip of the toxin loop, however, they played a critical role in the interactions with their molecular targets (Bourne et al., 2005). These mutations seem to induce a significant change in the structure of the toxin which is mostly due to additional hydrophobic interactions between Ile38 and the aliphatic side chains of the toxin that may induce a slight movement of the side chain surrounding Ser/Ile38.

Neurotoxin II (NTII) from Naja oxiana is a potent blocker of nAChR. NTII is a short chain α-neurotoxin which consists of 61 amino acid residues and four disulfide bridges (**Figure 4K**). A computational model for examining the interactions of NTII with the Torpedo californica nAChR has been studied (Mordvintsev et al., 2005). The model showed that the binding of the short α-neurotoxin occurs by rearranging the aromatic residues in the binding pocket. The insertion of the loop II into the binding pocket of a nAChR induces the neurotoxin activity and significantly determines the toxin-receptor interactions, while loop I and III contact the receptor residues by their tips only and determine the immunogenicity of the short neurotoxins. In the model, the Arg32 NTII residue forms an ionic pair with Trp149 from the nAChR and is observed as the strongest interaction. Hydrogen bond interactions such as Asp30 from loop II with Tyr198 of the nAChR and Lys46 from loop III with Thr191 of the nAChR complement the ionic interaction between NTII and its target receptor. Aside from hydrogen bonds, van-der-Waals interactions were also observed at the fingertip of loop II amino acid residues (Lys25, Trp27, Trp28, Ser29, His31, Gly33, Thr34 and Arg38) (Mordvintsev et al., 2005).

The structure of neurotoxin b (NTb), a long neurotoxin from Ophiophagus hannah, has been elucidated (Peng et al., 1997; **Figure 4L**). Conserved residues in loop II also play an important role in the toxicity of the long neurotoxins by making ionic interactions between toxin and receptor. Positively charged residues Trp27, Lys24 and Asp28 are highly conserved residues in the long neurotoxins. Furthermore, a modification of the Trp27 in the long neurotoxin analog of NTb from king cobra venom led to a significant loss in neurotoxicity. The additional disulphide bridge in loop II of long neurotoxins does not affect the toxin activity. Nevertheless, cleavage of the additional disulphide bridge in loop II can disrupt the positively charged cluster at the tip of loop II. Changes in loop II conformation will affect the binding of the long neurotoxin to the target receptor resulting the loss of neurotoxicity (Peng et al., 1997).

Long and short neurotoxins show sequence homology and similar structure. Previous studies show that many residues located at the tip of loop II are conserved in both short and long neurotoxins. It is consisting of the long central β-sheets forming three loops and globular core. From the studies of α-bungarotoxin and α-cobratoxin, the least conserved regions of the long neurotoxin are the C-terminal and the first loop (Walkinshaw et al., 1980; Juan et al., 1999; Dutertre et al., 2017). However, significant differences between long-chain neurotoxin and short chain neurotoxin are indicated by the immunological reactivity. Many of the residues involved in the antibody-long neurotoxins binding are located in loop II, loop III, and in the Cterminal, while in short neurotoxins the antibody's epitope makes interactions with the loop I and loop II (Engmark et al., 2016).

# THERAPEUTIC IMPLICATIONS

# Treating Snake Envenomation

Animal-derived antivenoms are considered the only specific therapy available for treating snakebite envenoming (Maduwage and Isbister, 2014; Slagboom et al., 2017; Ainsworth et al., 2018). These consist of polyclonal immunoglobulins, such as intact IgGs or F(ab')2, or Fab fragments (Ouyang et al., 1990; Maduwage and Isbister, 2014; Roncolato et al., 2015), derived from hyperimmune animal serum/plasma (typically horse or sheep). When used rapidly and appropriately, they are capable of neutralizing life-threatening systemic envenoming, for example pathologies such as venom-induced coagulopathy, hemorrhage, neurotoxic effects, and/or hypotensive shock (Warrell, 1992; Calvete et al., 2009; Maduwage and Isbister, 2014; Slagboom et al., 2017; Ainsworth et al., 2018).