Citation: Avella, I.; Damm, M.;
Freitas, I.; Wüster, W.; Lucchini, N.;
Zuazo, Ó.; Süssmuth, R.D.;
Martínez-Freiría, F. One Size Fits
All—Venomics of the Iberian Adder
(Vipera seoanei, Lataste 1878) Reveals
Low Levels of Venom Variation
across Its Distributional Range.
Toxins 2023,15, 371. https://
doi.org/10.3390/toxins15060371
Received: 29 April 2023
Revised: 18 May 2023
Accepted: 30 May 2023
Published: 1 June 2023
Copyright: © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).
toxins
Article
One Size Fits All—Venomics of the Iberian Adder (Vipera
seoanei, Lataste 1878) Reveals Low Levels of Venom Variation
across Its Distributional Range
Ignazio Avella 1,2,3,*,†, Maik Damm 4,† , Inês Freitas 1,2,3, Wolfgang Wüster 5, Nahla Lucchini 1,2,3,
Óscar Zuazo 6, Roderich D. Süssmuth 4and Fernando Martínez-Freiría1,3,*
1CIBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos, InBIO Laboratório Associado,
[email protected] (N.L.)
2Departamento de Biologia, Faculdade de Ciências, Universidade do Porto, 4099-002 Porto, Portugal
3BIOPOLIS Program in Genomics, Biodiversity and Land Planning, CIBIO, Campus de Vairão,
4485-661 Vairão, Portugal
4Institut für Chemie, Technische Universität Berlin, Straße des 17. Juni 124, 10623 Berlin, Germany;
[email protected] (M.D.)
5Molecular Ecology and Evolution at Bangor, School of Natural Sciences, Bangor University,
Bangor LL57 2UW, UK; w.wuster@bangor.ac.uk
6Calle La Puebla 1, 26250 Santo Domingo de la Calzada, Spain
*Correspondence: [email protected] (I.A.); [email protected] (F.M.-F.)
† These authors contributed equally to this work.
Abstract:
European vipers (genus Vipera) are medically important snakes displaying considerable venom
variation, occurring at different levels in this group. The presence of intraspecific venom variation, however,
remains understudied in several Vipera species. Vipera seoanei is a venomous snake endemic to the northern
Iberian Peninsula and south-western France, presenting notable phenotypic variation and inhabiting
several diverse habitats across its range. We analysed the venoms of 49 adult specimens of V. seoanei
from 20 localities across the species’ Iberian distribution. We used a pool of all individual venoms to
generate a V. seoanei venom reference proteome, produced SDS-PAGE profiles of all venom samples, and
visualised patterns of variation using NMDS. By applying linear regression, we then assessed presence
and nature of venom variation between localities, and investigated the effect of 14 predictors (biological,
eco-geographic, genetic) on its occurrence. The venom comprised at least 12 different toxin families, of
which five (
i.e., PLA2
, svSP, DI, snaclec, svMP) accounted for about 75% of the whole proteome. The
comparative analyses of the SDS-PAGE venom profiles showed them to be remarkably similar across the
sampled localities, suggesting low geographic variability. The regression analyses suggested significant
effects of biological and habitat predictors on the little variation we detected across the analysed V. seoanei
venoms. Other factors were also significantly associated with the presence/absence of individual bands in
the SDS-PAGE profiles. The low levels of venom variability we detected within V. seoanei might be the
result of a recent population expansion, or of processes other than directional positive selection.
Keywords: snake venom; viper; bottom-up proteomics; regional variation; Iberian Peninsula
Key Contribution:
The current study provides the first proteomic analysis of the venom of the
medically important Iberian adder (Vipera seoanei). Bottom-up venomics showed that phospholipases
A
2
(PLA
2
s), snake venom serine proteases (svSPs), disintegrins (DIs), snake venom C-type lectins
and C-type lectin-related proteins (snaclecs), and snake venom metalloproteinases (svMPs) account
for roughly 75% of the adult V. seoanei venom proteome. Comparative analyses of the SDS-PAGE
profiles of whole venoms collected from 49 adult individuals from 20 different localities across the
species’ Iberian range revealed them to be remarkably similar. Regression analyses found significant
correlations between the low levels of venom variation detected and a number of biological and
eco-geographic factors.
Toxins 2023,15, 371. https://doi.org/10.3390/toxins15060371 https://www.mdpi.com/journal/toxins
Toxins 2023,15, 371 2 of 22
1. Introduction
Snake venoms comprise a mixture of proteins, peptides, organic molecules
(
e.g., carbohydrates,
lipids), and salts in an aqueous medium [
1
,
2
], and are able to dis-
rupt the homeostatic processes of the envenomated organism [
3
–
5
]. To date, more than
60 protein
families that compose snake venoms have been identified [
6
]. This complexity
of components harbours the potential for extreme variation, which has been found to occur
frequently and at all taxonomic levels (see [
7
–
9
]). Indeed, broad compositional differences
have been found between venoms of snakes belonging to different families (e.g., Elapidae
and Viperidae [
6
,
10
]), between different genera belonging to the same family (e.g., Aus-
tralian elapids [
11
,
12
]; Old World vipers [
13
]), and between congeneric species (e.g., within
the genera Bothriechis [
14
] and Bothrops [
15
]). Moreover, variation in venom composition
and activity has also been found to occur at the intra-specific level, between juveniles and
adults (e.g., [
16
–
18
]), between sexes (e.g., [
19
–
21
]), and between individuals from different
geographic localities (e.g., [
22
–
24
]), where it may persist even in the face of extensive gene
flow [25].
Extensive evidence suggests that the primary function of snake venoms is prey sub-
jugation, and that the main driver behind the dynamism of its variation is adaptation
to diet [
26
–
28
]. In this scenario, the occurrence of venom variation between conspecific
snakes originating from different areas (i.e., geographic/regional variation) has generally
been associated with differences in their feeding ecology (e.g., [
28
,
29
]. Indeed, specimens
living in environmentally distinct areas are likely under different ecological pressures
(
e.g., different
prey availability), which likely determines changes in their diet. Considering
the critical adaptive value and fast evolutionary rates of snake venom [
30
,
31
], it stands to
reason that these changes might drive the emergence of compositional and/or functional
venom variation [8].
Multiple studies have then supported the adaptive value of geographic variation in
the venom composition of different snake taxa in the face of prey type, availability, and/or
susceptibility to venom (e.g., [
32
,
33
]), and a number of them found correlations between
the occurrence of snake venom variation and different environments. For example, Sousa
et al. (2017) found that Bothrops atrox venoms of specimens from different habitats within
the Amazon presented significantly different activities [
34
]. Similarly, Zancolli et al. (2019)
found a strong correlation between the dramatic venom variation occurring in Crotalus
scutulatus [
35
,
36
] and environmental heterogeneity [
25
]. A correlation between venom
phenotype and climatic variables in the same species was also found by Strickland et al.
(2018) [
37
]. In these studies, the authors hypothesised that climatic/environmental variables
could be determining changes in factors such as snake prey distribution and/or physiol-
ogy, likely affecting the feeding ecology of the snake species, thus ultimately supporting
adaptation to diet as a possible cause of the detected venom variation.
European vipers (genus Vipera) are among the snakes notable for extensive venom vari-
ation, including at the interspecific level (see [
13
]), as well as within some species. In partic-
ular, such variation corresponds to the occurrence of profound clinical differences between
geographically close populations of species such as Vipera aspis and
Vipera berus [38,39]
.
However, several Vipera species remain inadequately studied, especially in terms of the
presence of intraspecific venom variation, which represents a significant knowledge gap
given the medical importance of the genus [40–42].
One of the hitherto neglected species of Vipera is the Iberian adder, Vipera seoanei,
Lataste 1879, a venomous snake belonging to the family Viperidae (subfamily Viperinae)
and a member of the clade of vipers with Euro-Siberian affinity (i.e., Pelias [
43
]). Endemic
to the northern Iberian Peninsula and south-western France, this species mainly inhabits
areas with an Atlantic climate, typically occupying humid ecotones between meadows
and forests and zones with abundant basal vegetation from sea level to about 1900 m
of altitude (
i.e., Cantabrian
Mountains [
44
]). Sister species of the common adder V. berus
(see [
43
]), V. seoanei displays low intraspecific genetic variability, likely as result of a
late Pleistocene expansion from north-western Iberian refugia [
45
]. Despite its shallow
Toxins 2023,15, 371 3 of 22
genetic structure, V. seoanei shows considerable variation in biometric and pholidotic
traits across its range [
46
,
47
]. The species is notable for high levels of polymorphism in
body colouration, with five colour phenotypes currently recognised (i.e., bilineata,cantabrica,
classic,melanistic,uniform). These appear to be geographically structured and not concordant
with mitochondrial DNA haplotypes [45,48].
The diet of V. seoanei mainly comprises small mammals and, less frequently, reptiles,
amphibians, arthropods, and birds [
49
–
51
]. Like many other species of the genus Vipera
(
e.g., V. aspis [52]
;Vipera latastei [
53
]), V. seoanei exhibits ontogenetic shift in diet composition,
with ectotherms (e.g., reptiles, amphibians) constituting more than 70% of the diet of
juvenile vipers, whereas small mammals (i.e., shrews and rodents) account for roughly 90%
of the diet of the adults [
50
]. Interestingly, significant correlations have been found between
climatic and habitat conditions and differences in the frequency of consumption of the prey
items V. seoanei mainly feeds on (i.e., amphibians, reptiles, and small mammals [
50
]). Such
environment-related variation in diet composition, also observed in other Vipera species
(e.g., Vipera ammodytes [
54
]; V. latastei [
55
]), suggests the presence of dietary differences
across the species range.
The Iberian adder is recognised as medically important by the World Health Organi-
zation [
42
] and is considered one of the venomous snake species of major clinical relevance
in Europe [
40
,
41
]. An early study [
56
] aiming at investigating the toxicity of V. seoanei
venom across its Spanish distribution suggested the presence of a West–East gradient,
with venoms from the western populations presenting higher lethal potencies than those
collected from eastern populations (e.g., Galicia (W Spain): LD
50
= 6.9–9.9
µ
g per 20 g
mouse; Basque Country (E Spain): LD
50
= 23.1–23.6
µ
g per 20 g mouse). A more recent
study [
57
] reported toxicity values for V. seoanei of Portuguese origin (i.e., LD
50
= 9.7
µ
g
per 18–20 g mouse) comparable to those previously obtained by Detrai et al. (1990) [
56
]
for individuals from Galicia and north of León, potentially supporting the higher toxic-
ity of the westernmost populations. Interestingly, although V. seoanei venoms collected
from different localities across the species’ Spanish distribution presented differences in
LD
50
values, the corresponding proteinograms appeared to be very similar, suggesting
low levels of interpopulational divergence in venom compositions [
56
]. In spite of this
intriguing scenario, and of the recognised medical relevance of V. seoanei, no other studies
regarding this species’ venom have been conducted, and further information about it is
currently unavailable.
In the present work, we (i) provide the first proteomic characterisation of the venom of
V. seoanei; (ii) assess the level of venom variation across the viper’s range in the Iberian Penin-
sula; (iii) test for associations between biological, geographic, genetic, and eco-geographic
factors and venom variation within this species across its ecologically and physiographically
diverse range.
2. Results
2.1. Protein Composition of V. seoanei Venom
A total of 55 peaks were identified within the chromatographic profile of the V. seoanei
venom pool (Figure 1A). The chromatogram shows several broad peaks with high abun-
dance (e.g., peaks 6, 26, 32, 33, 36, 40, 54, and 55) and complex regions with smaller peaks,
e.g., eluting between 35 and 50 min, as well as between 80 and 95 min. The subsequent
reducing SDS-PAGE analysis revealed that the highest HPLC peaks are often dominated
by few, highly abundant bands (see Figure 1B).
The major toxin families in the venom proteome included phospholipases of type A
2
(PLA
2
) and snake venom serine proteases (svSP), followed by the slightly less abundant dis-
integrins (DI), snake venom C-type lectins and C-type lectin-related proteins (snaclec), and
snake venom metalloproteinases (svMP). Other toxin families such as vascular endothelial
growth factors (VEGF), cysteine-rich secretory proteins (CRISP), L-amino-acid oxidases
(LAAO), Kunitz-type inhibitors (KUN), and venom nerve growth factors (NGF) were also
detected with relative abundances of 6% or less. The peptidic part, made of inhibitors of
Toxins 2023,15, 371 4 of 22
svMPs (svMP-i), natriuretic peptides (NP), and other peptides composed about 11% of the
venom proteome. Non-annotated proteins (NA) accounted for less than 0.2% of the venom
proteome. For further details, see Figure 2and Tables S1 and S2.
Among the PLA
2
s, we detected a variety of different proteoforms, some highly similar
to the non-enzymatic PLA
2
homologue S49 ammodytin L(2) [Q6A394], and two enzymatic
active D49 forms, namely ammodytin I1 and I2 [Q910A1;P34180]. Therefore, V. seoanei
venom includes basic, neutral, and acidic PLA
2
s. The main svSP we identified (6%, peak
44–45) shows high similarity with nikobin [E5AJX2] and an RVV-V
γ
homolog [P18965]
(4%, peak 43–45). Similarities were also found with other svSPs, such as SP-4/5/8 isoforms
[A0A1I9KNR8; A0A1I9KNR5; A0A6B7FPJ0] of Vipera ammodytes ammodytes.
Toxins 2023, 15, x FOR PEER REVIEW 4 of 22
Figure 1. Fractionation of the V. seoanei venom pool. The figure shows RP-HPLC profile (A), with
peak 0 corresponding to the injection peak, and SDS-PAGE Coomassie-stained profile (B) of the
venom pool under reducing conditions. PAGE line nomenclature is based on RP-HPLC fractions.
Labelled bands were cut, subjected to tryptic digestion, and analysed with LC-MS.
The major toxin families in the venom proteome included phospholipases of type A
2
(PLA
2
) and snake venom serine proteases (svSP), followed by the slightly less abundant
disintegrins (DI), snake venom C-type lectins and C-type lectin-related proteins (snaclec),
and snake venom metalloproteinases (svMP). Other toxin families such as vascular endo-
thelial growth factors (VEGF), cysteine-rich secretory proteins (CRISP), L-amino-acid ox-
idases (LAAO), Kunitz-type inhibitors (KUN), and venom nerve growth factors (NGF)
were also detected with relative abundances of 6% or less. The peptidic part, made of in-
hibitors of svMPs (svMP-i), natriuretic peptides (NP), and other peptides composed about
11% of the venom proteome. Non-annotated proteins (NA) accounted for less than 0.2%
of the venom proteome. For further details, see Figure 2 and Tables S1 and S2.
Among the PLA
2
s, we detected a variety of different proteoforms, some highly simi-
lar to the non-enzymatic PLA
2
homologue S49 ammodytin L(2) [Q6A394], and two enzy-
matic active D49 forms, namely ammodytin I1 and I2 [Q910A1;P34180]. Therefore, V.
seoanei venom includes basic, neutral, and acidic PLA
2
s. The main svSP we identified (6%,
peak 44–45) shows high similarity with nikobin [E5AJX2] and an RVV-Vγ homolog
[P18965] (4%, peak 43–45). Similarities were also found with other svSPs, such as SP-4/5/8
isoforms [A0A1I9KNR8; A0A1I9KNR5; A0A6B7FPJ0] of Vipera ammodytes ammodytes.
All svMP fragments were annotated using MS as members of the PIII-svMP subfam-
ily, including ions from SDS bands with molecular weights of 30 kDa or less. This range
Figure 1.
Fractionation of the V. seoanei venom pool. The figure shows RP-HPLC profile (
A
), with
peak 0 corresponding to the injection peak, and SDS-PAGE Coomassie-stained profile (
B
) of the
venom pool under reducing conditions. PAGE line nomenclature is based on RP-HPLC fractions.
Labelled bands were cut, subjected to tryptic digestion, and analysed with LC-MS.
Toxins 2023,15, 371 5 of 22
Toxins 2023, 15, x FOR PEER REVIEW 5 of 22
of lower molecular svMP masses usually originates from mature PI or processed PII frag-
ments but is uncommon for PIII. Sequence analysis showed that peaks 28–30 refer only to
the DC domain of a PIII-svMP, which is strong evidence for PIIIe, since no PII could be
assigned [58,59]. The DIs detected in the venom of V. seoanei were mostly composed of
homologs of dimeric RGD-disintegrins, such as VA6 [P0C6A5], lebein-1α [P83253], and
VB7B [P0C6A7]. Snaclecs were constantly detected in two- or four-band patterns at late
retention times of 85–110 min. These toxins are known to form multimeric structures with
one, two (αβ), or four (αβγδ) subunits in different complexities, such as (αβ)4 heterooc-
tamers, and in combination with svMP P-IIIe [58,60]. We identified several homologs to
different snaclec subunits isolated from Vipera ammodytes ammodytes and from the Eura-
sian viper genera Daboia and Macrovipera (see Tables S1 and S2).
Figure 2. Reference composition of Vipera seoanei venom. The pie chart displays the relative abun-
dances of the toxin families found in the proteome of the V. seoanei venom pool. PLA2, phospho-
lipases A2; svSP, snake venom serine proteases; DI, disintegrins; snaclec, snake venom C-type lectins
and C-type lectin-related proteins; svMP, snake venom metalloproteinases; VEGF, vascular endo-
thelial growth factors; CRISP, cysteine-rich secretory proteins; LAAO, L-amino-acid oxidases; KUN,
Kunitz-type inhibitors; NGF, venom nerve growth factors; svMP-i, svMP inhibitors; NP, natriuretic
peptides; NA, non-annotated proteins.
The less abundant toxins include two classes of growth factors, namely NGF with
two populations of 15 kDa and VEGF, almost exclusively identified as Vammin-1 homo-
logs [APB93447]. A single band population of 25 kDa (peak 36–37) contains CRISP homo-
logs of V. berus CRVP [B7FDI1]. The peptide part (11%) of the produced venom proteome
is dominated by small tripeptidic svMP-is, as is the case for most peptidomes of Old
World viper venoms [13]. Also included in the peptidic part are natriuretic peptides and
fragments of bradykinin-potentiating and C-type natriuretic peptides (see Tables S1 and
S2).
2.2. Assessing Geographic Venom Variation
2.2.1. Analysis of the SDS-PAGE Profiles
In total, 49 SDS-PAGE whole venom profiles were produced from the individual V.
seoanei venoms collected (Figure S1). By visually comparing the obtained profiles, we did
not detect stark differences in terms of the presence/absence of bands, and instead noticed
a considerable overall level of similarity between them. The least diverse profiles had 13
evident bands, whereas the most diverse had 18. Six bands were difficult to identify reli-
ably. Nine bands were present in all venom profiles so were uninformative, and were thus
excluded from further analyses.
Figure 2.
Reference composition of Vipera seoanei venom. The pie chart displays the relative abun-
dances of the toxin families found in the proteome of the V. seoanei venom pool. PLA
2
, phospholipases
A
2
; svSP, snake venom serine proteases; DI, disintegrins; snaclec, snake venom C-type lectins and
C-type lectin-related proteins; svMP, snake venom metalloproteinases; VEGF, vascular endothelial
growth factors; CRISP, cysteine-rich secretory proteins; LAAO, L-amino-acid oxidases; KUN, Kunitz-
type inhibitors; NGF, venom nerve growth factors; svMP-i, svMP inhibitors; NP, natriuretic peptides;
NA, non-annotated proteins. The sum of the percentages does not match 100% because of rounding.
Detailed percentages are reported in Table S1.
All svMP fragments were annotated using MS as members of the PIII-svMP subfamily,
including ions from SDS bands with molecular weights of 30 kDa or less. This range of lower
molecular svMP masses usually originates from mature PI or processed PII fragments but is
uncommon for PIII. Sequence analysis showed that peaks 28–30 refer only to the DC domain
of a PIII-svMP, which is strong evidence for PIIIe, since no PII could be assigned [
58
,
59
].
The DIs detected in the venom of V. seoanei were mostly composed of homologs of dimeric
RGD-disintegrins, such as VA6 [P0C6A5], lebein-1
α
[P83253], and VB7B [P0C6A7]. Snaclecs
were constantly detected in two- or four-band patterns at late retention times of 85–110 min.
These toxins are known to form multimeric structures with one, two (
αβ
), or four (
αβγδ
)
subunits in different complexities, such as (
αβ
)4 heterooctamers, and in combination
with svMP P-IIIe [
58
,
60
]. We identified several homologs to different snaclec subunits
isolated from Vipera ammodytes ammodytes and from the Eurasian viper genera Daboia and
Macrovipera (see Tables S1 and S2).
The less abundant toxins include two classes of growth factors, namely NGF with two
populations of 15 kDa and VEGF, almost exclusively identified as Vammin-1 homologs
[APB93447]. A single band population of 25 kDa (peak 36–37) contains CRISP homologs
of V. berus CRVP [B7FDI1]. The peptide part (11%) of the produced venom proteome is
dominated by small tripeptidic svMP-is, as is the case for most peptidomes of Old World
viper venoms [
13
]. Also included in the peptidic part are natriuretic peptides and fragments
of bradykinin-potentiating and C-type natriuretic peptides (see Tables S1 and S2).
2.2. Assessing Geographic Venom Variation
2.2.1. Analysis of the SDS-PAGE Profiles
In total, 49 SDS-PAGE whole venom profiles were produced from the individual
V. seoanei
venoms collected (Figure S1). By visually comparing the obtained profiles, we
did not detect stark differences in terms of the presence/absence of bands, and instead
noticed a considerable overall level of similarity between them. The least diverse profiles
had 13 evident bands, whereas the most diverse had 18. Six bands were difficult to identify
reliably. Nine bands were present in all venom profiles so were uninformative, and were
thus excluded from further analyses.
Toxins 2023,15, 371 6 of 22
2.2.2. Band Analysis and NMDS
The final binary matrix (Table S3) used to assess the diversity between the 49 SDS-
PAGE whole venom profiles included 10 polymorphic (i.e., not appearing in every pro-
file) bands. Their number varied between four (samples 19VS076, 19VS185, 21VS005,
and 21VS006) and nine (samples 19VS450, 19VS451, and 21VS022) across the profiles
(see Table S3).
The three-dimensional NMDS analysis achieved an overall stress value of 0.08, indi-
cating that the method could suitably represent the differences between the venom profiles.
The NMDS ordination plot generated using the individual NMDS scores on the first two
axes (i.e., NMDS1 and NMDS2; Figure 3) showed low dispersion of the samples across the
defined space and a considerable level of overlap.
Figure 3.
Non-metric multidimensional scaling (NMDS) ordination plot of the 49 V. seoanei venom
SDS-PAGE profiles. Teal circles represent the venom profiles. Notice that some circles correspond to
overlapping profiles (up to six).
2.2.3. Regression Analysis
Of the 14 predictors tested using the NMDS1 scores as the response variable, only four
were significant in univariate linear regression analyses, namely SVL (p= 0.016), COLOUR
(p= 0.028), POPULATION (p= 0.001), and FOREST (p= 0.035). In the univariate linear
regression analyses testing the effect of each of the 14 predictors on the NMDS2 scores, only
the three predictors SEX, COLOUR, and FOREST were significant (p= 0.022, p= 0.041, and
p= 0.033, respectively). Details of the univariate linear regression models performed are
reported in Table 1.
In the multiple regression model comprising individual NMDS1 scores as the response
variable and SVL, COLOUR, and FOREST as predictors, no significant effects were detected.
In the multiple regression model built including individual NMDS2 scores as the response
variable and the three predictors SEX, COLOUR, and FOREST, only SEX had a significant
effect (p= 0.045). The p-value relative to Moran’s I calculated for all variables included
in the two multiple regression models was >0.1, indicating the absence of statistically
significant spatial autocorrelation. Details of the multiple regression models used are
reported in Table 2.
Toxins 2023,15, 371 7 of 22
Table 1.
Results of the univariate linear regression models using the individual NMDS1 and NMDS2
scores. For each predictor tested, sum of squares (Sum Sq), degrees of freedom (Df), F value (F), and
associated p-value (p) are reported. Significant predictors and corresponding p-values are in bold.
NMDS1 NMDS2
Predictor Sum Sq Df F p Sum Sq Df F p
SVL 0.263 1 6.351 0.016 <0.001 1 0.013 0.911
SEX 0.081 1 1.794 0.187 0.205 1 5.626 0.022
COLOUR 0.397 3 3.312 0.028 0.324 3 2.983 0.041
POPULATION 1.538 19 3.508 0.001 1.032 19 1.779 0.079
GEN1 0.029 1 0.635 0.429 <0.001 1 <0.001 0.997
GEN2 0.015 1 0.322 0.573 0.016 1 0.403 0.529
BIO1 0.061 1 1.334 0.254 0.001 1 0.031 0.861
BIO5 0.145 1 3.314 0.075 0.123 1 3.379 0.072
BIO12 0.133 1 3.019 0.089 <0.001 1 0.009 0.925
BIO14 0.002 1 0.048 0.823 0.015 1 0.363 0.549
AGRIC 0.016 1 0.345 0.559 0.102 1 2.646 0.111
FOREST 0.201 1 4.702 0.035 0.178 1 4.793 0.033
MOOR 0.002 1 0.051 0.823 0.004 1 0.107 0.746
PASTURE <0.001 1 <0.001 0.981 <0.001 1 0.011 0.912
Table 2.
Results of the multiple regression models. For each predictor included in the model, sum
of squares (Sum Sq), degrees of freedom (Df), F value (F), and associated p-values (p) are reported.
Residual sum of squares (Res. Sum Sq) and residual degrees of freedom (Res. Df) of the full models
are also presented. Significant predictors and corresponding p-values are in bold.
Response Predictor Sum Sq Df F p Res. Sum Sq Res. Df
NMDS1 SVL 0.084 1 2.142 0.151 1.641 42
COLOUR 0.208 3 1.773 0.167
FOREST 0.021 1 0.536 0.468
NMDS2 SEX 0.136 1 4.255 0.045 1.347 42
COLOUR 0.223 3 2.317 0.089
FOREST 0.084 1 2.611 0.114
Following the same approach, we used single-predictor generalised linear models
(GLMs) to test the effect of each of the 14 selected predictors on the presence/absence of each
of the 10 polymorphic bands retrieved from the individual whole venom profiles. While no
significant effects on the presence/absence of bands 2, 5, and 6 were detected, at least one
predictor was significantly correlated with the presence/absence of the remaining seven
polymorphic bands. The effect of only one predictor was significant for bands 1 (
i.e., SVL
,
p= 0.022
) and 3 (i.e., POPULATION, p= 0.021); thus, we could not perform multiple-
predictor GLMs for these two bands. Similarly, we could not build a multiple-predictor
GLM for band 7 because single-predictor GLMs used for this band found significance
only for the predictors COLOUR (p= 0.002) and POPULATION (p= 0.026), which could
not be included in the same model because of the high level of association between them
(Cramér’s
V = 0.899
). Of all single predictor-GLMs used, some provided unreliable fit and
low prediction power and, thus, were not considered in subsequent analyses. Additional
details of the single-predictor binomial GLMs used are reported in Table S4 and Figure S2.
Multiple-predictor GLMs were built for bands 4, 8, 9, and 10 by including all the
predictors that had a significant effect in the single-predictor GLMs for these four bands.
Significant effects were found for the predictors SEX (p= 0.018) and AGRIC (p= 0.044)
on the presence/absence of band 4, and for the predictor COLOUR (p< 0.001) on the
presence/absence of band 10 (Table 3). Specifically, band 4 was significantly less detected
Toxins 2023,15, 371 8 of 22
in females than in males (Figure 4A), and the probability of detecting it appeared to
be inversely related to the percentage of agricultural areas (Figure 4B). Concerning the
significant effect of the predictor COLOUR on the presence/absence of band 10, the four
V. seoanei colour phenotypes considered (i.e., bilineata,cantabrica,classic,melanistic) were
related differently to the response variable (Figure 4C).
Table 3.
Results of the multiple predictor binomial GLMs. The models investigate the probability
of occurrence of bands 4, 8, 9, and 10 in an SDS-PAGE venom profile. Likelihood ratio chi-square
(LR
χ2
), degrees of freedom (Df), and associated p-values (p) are reported. Significant predictors and
corresponding p-values are in bold. The asterisks (*) indicate predictors that, although significant,
were not considered because the corresponding models did not converge.
Response Predictors LR χ2Df p
Band 4 SEX 5.544 1 0.018
AGRIC 4.068 1 0.044
Band 8 SVL 0.529 1 0.467
COLOUR 3.638 3 0.303
BIO5 0.059 1 0.807
BIO12 0.000 1 0.996
Band 9 SEX * 20.679 1 <0.001
GEN2 * 16.784 1 <0.001
Band 10 COLOUR 16.612 3 <0.001
BIO5 2.387 1 0.122
FOREST 3.765 1 0.052
Toxins 2023, 15, x FOR PEER REVIEW 8 of 22
than in males (Figure 4A), and the probability of detecting it appeared to be inversely
related to the percentage of agricultural areas (Figure 4B). Concerning the significant effect
of the predictor COLOUR on the presence/absence of band 10, the four V. seoanei colour
phenotypes considered (i.e., bilineata, cantabrica, classic, melanistic) were related differently
to the response variable (Figure 4C).
Figure 4. Multiple-predictor GLM predictions of occurrence of bands 4 and 10 in individual SDS-
PAGE venom profiles. The panels display the predicted probability of occurrence of band 4 in rela-
tion to the sex of the viper (A) and amount of cultivated fields (B), and of band 10 in relation to four
different colour phenotypes displayed by V. seoanei (C). The high standard deviations are attributa-
ble to model limitations owing to small sample size.
The effect of predictors SEX and GEN2 on the presence/absence of band 9 was also
significant, but the multiple-predictor GLM for this band showed unreliable fit and low
prediction power and was, therefore, not considered. None of the predictors included in
the model with the presence/absence of band 8 as the response variable were significant.
The p-value relative to Moran’s I calculated for all variables included in the four mul-
tiple-predictor GLMs was >0.1, meaning that there was no statistically significant spatial
autocorrelation in the data.
Table 3. Results of the multiple predictor binomial GLMs. The models investigate the probability of
occurrence of bands 4, 8, 9, and 10 in an SDS-PAGE venom profile. Likelihood ratio chi-square (LR
χ2), degrees of freedom (Df), and associated p-values (p) are reported. Significant predictors and
corresponding p-values are in bold. The asterisks (*) indicate predictors that, although significant,
were not considered because the corresponding models did not converge.
Response Predictors LR χ2 Df p
Band 4 SEX 5.544 1 0.018
AGRIC 4.068 1 0.044
Band 8 SVL 0.529 1 0.467
COLOUR 3.638 3 0.303
BIO5 0.059 1 0.807
BIO12 0.000 1 0.996
Band 9 SEX * 20.679 1 <0.001
GEN2 * 16.784 1 <0.001
Band 10 COLOUR 16.612 3 <0.001
BIO5 2.387 1 0.122
FOREST 3.765 1 0.052
3. Discussion
3.1. Protein Composition of V. seoanei Venom
Through the application of bottom-up venomics, we were able to provide the first
characterisation of the protein components present in the venom of V. seoanei. As shown
Figure 4.
Multiple-predictor GLM predictions of occurrence of bands 4 and 10 in individual SDS-
PAGE venom profiles. The panels display the predicted probability of occurrence of band 4 in relation
to the sex of the viper (
A
) and amount of cultivated fields (
B
), and of band 10 in relation to four
different colour phenotypes displayed by V. seoanei (
C
). The high standard deviations are attributable
to model limitations owing to small sample size.
The effect of predictors SEX and GEN2 on the presence/absence of band 9 was also
significant, but the multiple-predictor GLM for this band showed unreliable fit and low
prediction power and was, therefore, not considered. None of the predictors included in
the model with the presence/absence of band 8 as the response variable were significant.
The p-value relative to Moran’s I calculated for all variables included in the four
multiple-predictor GLMs was >0.1, meaning that there was no statistically significant
spatial autocorrelation in the data.
3. Discussion
3.1. Protein Composition of V. seoanei Venom
Through the application of bottom-up venomics, we were able to provide the first
characterisation of the protein components present in the venom of V. seoanei. As shown by
Toxins 2023,15, 371 9 of 22
the high number of peaks and bands retrieved from the RP-HPLC and SDS-PAGE profiles
produced for the venom pool (see Figure 1), the Iberian adder’s venom comprises at least
12 different toxin families. Five of these families, namely PLA
2
, svSP, DI, snaclec, and svMP,
constitute about 75% of the full venom proteome (see Figure 2and
Tables S1 and S2
). Simi-
lar compositional patterns have been described for the venoms of V. seoanei’s sister species
V. berus [
61
] and the Iberian endemic V. latastei [
16
], in line with the general composition of
Viperinae venoms, of which the four toxin families PLA
2
, svSP, snaclec, and svMP typically
compose 60–90% [13].
The effects generally caused by the five major toxin families composing the venom of
V. seoanei are concordant with the mainly haemorrhagic and cytotoxic symptoms typically
reported for viper envenomation [
62
,
63
]. Specifically, DI, snaclec, and svSP can affect blood
coagulation, fibrinolysis, angiogenesis, and platelet aggregation [
64
–
66
], and svMP (espe-
cially class PIII) are known to affect the coagulation cascade and platelet aggregation and
to cause severe haemorrhage [
67
,
68
]. Phospholipases of the PLA
2
type, the most abundant
component of the analysed venom pool, constitute a very diverse toxin family [
69
] that can
produce a plethora of different effects, including myotoxicity, cardiotoxicity, cytotoxicity,
and coagulotoxicity [
70
]. Of the several different PLA
2
proteoforms we detected in the
V. seoanei venom, some of them appear to correspond to ammodytin L and ammodytin I,
suggesting myotoxic and haemolytic effects (see [69]).
Among the less abundant toxins we found in V. seoanei venom, KUN and LAAO have
been reported to interfere with platelet aggregation, fibrinolysis, and angiogenesis [
71
–
73
],
with the latter also potentially being disrupted by CRISP and VEGF [
74
,
75
]. Finally, NGF
can increase vascular permeability in the envenomated organism, thus aiding the spread of
other toxins [
76
], while NP can decrease myocardial contractility and cause hypotension,
rapidly leading to the loss of consciousness [77].
3.2. Low Levels of Geographic Variation
In order to assess the presence and extent of geographic venom variation within
V. seoanei, we performed comparative analyses of the SDS-PAGE whole venom profiles
obtained from the 49 individuals collected across the species’ distribution in northern Iberia.
By visually comparing the profiles in question (Figure S1), we noticed the lack of marked
differences among them. Additionally, in the NMDS ordination plot generated using the
individual NMDS scores on the first two axes (i.e., NMDS1 and NMDS2), most of the points
were clustered in an almost central position, several of them were overlapping, and we did
not detect a defined pattern of geographic venom variation (see Figure 3). Although the
Shapiro–Wilk test performed for NMDS1 and NMDS2 values showed that only the latter
did not depart significantly from normal distribution (W= 0.951, p= 0.041 and W= 0.966,
p= 0.168, respectively), the statistical models supported the classification of NMDS1 and
NMDS2 values as normally distributed (41% of probability, against 20% of probability of
following a beta distribution). Taken together, these results suggest the presence of low
levels of variation within the adult V. seoanei venoms considered.
Being an ecologically critical functional trait, snake venom is under strong selective
pressures, shaping its composition and activity to facilitate the snake’s survival [
1
,
8
]. The
numerous reported cases of snake venom varying between different areas, likely in re-
sponse to differences in ecological and environmental conditions (e.g., prey availability
and susceptibility to venom), appear to provide consistent support for the adaptive value
of snake venom variation (see [
27
,
32
]). Nonetheless, cases of geographic snake venom
variation being almost undetectable at the intraspecific level are known. For example, Hof-
mann et al. (2018) and Rautsaw et al. (2019) did not detect a defined pattern of geographic
variation in the venoms of Crotalus cerastes from Arizona and California, possibly due to
stabilising selection favouring generalist venom phenotypes [
78
–
80
]. Similarly, Margres
et al. (2015) found no significant variation in the expression of toxins and toxin genes
across individuals of Micrurus fulvius from Florida, perhaps as result of relaxed selective
constraints or a recent range expansion [
81
]. Therefore, it would be interesting to investigate
Toxins 2023,15, 371 10 of 22
if the low variation we observed among the SDS-PAGE V. seoanei venom profiles analysed
could be a consequence of factors potentially preventing local venom adaptation, such as
the recent population expansion from the north-west towards the eastern Iberian Peninsula
suggested for this species [45], or considerable levels of gene flow (see [79]).
Additionally, in light of the arms race between snakes and their prey, typically con-
sisting of prey evolving resistance to venom and snake venom evolving to bypass this
resistance [
82
–
84
], it has been suggested that balancing selection might favour a diverse set
of venom alleles over a single optimal venom genotype in order to prevent fixed venom
alleles from becoming ineffective due to evolved prey resistance [
36
]. In this scenario,
considering that adult Iberian adders are small mammal specialists [
50
], the low venom
variation detected might align with balancing selection acting to allow the effective sub-
jugation of their potentially coevolving preferred prey while also avoiding evolutionary
dead ends.
3.3. Correlates of Venom Variation
In our regression models, we implemented the same predictors Espasandín et al. (2022)
used to study the effect of eco-geographic variables on the trophic ecology of
V. seoanei [50]
.
Considering the strong link between snake venom variation and diet, we aimed to test
whether they could be at play in determining the occurrence and extent of venom variation
within this species. Although we detected low overall levels of venom variation, the
univariate linear regression analyses performed showed significant effects of body size
(
i.e., SVL
), colour phenotype, sex, locality of origin of the vipers, and percentage of forested
area on the NMDS1 and NMDS2 scores (see Table 1for details). Nonetheless, in the multiple
regression models, only the predictor SEX had a significant effect on the variation detected
across the analysed venoms (see Table 2).
In gape-limited animals such as snakes [
85
], sexual dimorphism in body and/or head
size can define the spectrum of prey items each sex can feed on [
86
]. While no sexual
dimorphism concerning head size has been detected in V. seoanei, significant intersexual
variation in body size (i.e., SVL) has been reported for this species [
46
]. Interestingly,
significant differences in feeding frequencies have been found between male and female
Iberian adders, with females feeding more often than males, and males reducing their
feeding rates as they grow [
50
]. Although these considerations might suggest a potential
influence of intersexual differences in feeding ecology on the low variation detected across
the V. seoanei venoms analysed, we did not identify sex-specific venom protein bands.
Additionally, despite the NMDS ordination plot indicating less divergence within female
venoms than within male venoms (see Figure 3), the male and female venom profiles did
not differ significantly in their complexity (i.e., number of bands; Mann–Whitney U= 284,
p= 0.991
). The significance of the sex and body size of the vipers could, therefore, be related
to factors not included in our analyses (e.g., seasonality, reproductive stage), which might
be unravelled by studies developed on a finer scale.
Considering the strong correlation between the two predictors POPULATION and
COLOUR (see Section 5.5), and that the five colour phenotypes currently recognised for
V. seoanei are geographically structured, the significance of these two predictors in our
analyses might refer to different local selective regimes acting on the analysed venoms.
For instance, individuals presenting the four colour phenotypes included in our analyses,
i.e., bilineata
,cantabrica,classic, and melanistic, are known to differ in body proportions,
likely due to them being subjected to different ecological pressures (e.g., climatic conditions,
predatory pressure) in the habitats where they occur [
46
]. Additionally, we suspect the
significance of the predictor FOREST to be related to changes in the feeding ecology of
V. seoanei
associated with this habitat type (e.g., the consumption of amphibians by Iberian
adders appears to increase in rainy and forested areas [
50
]). Considering the lack of marked
genetic distinctness within
V. seoanei
, these results could hint that the little venom variation
we detected might be related to ecological factors determining differences in local selective
pressures acting on distinct colour phenotypes, and/or local changes in prey availability.
Toxins 2023,15, 371 11 of 22
The single-predictor binomial generalised linear models (GLMs) used for the 10 poly-
morphic SDS-PAGE bands showed that the presence/absence of eight bands (i.e., 1, 2, 3, 4,
7, 8, 9, 10) was significantly correlated to at least one of the 14 predictors tested. Specifically,
the models showed variation in band presence/absence in relation to body size (SVL),
sex (SEX), colour phenotype (COLOUR), locality of origin (POPULATION), the PC2 of
the SPCA of the interpolated genetic distances (GEN2), maximum temperature of the
warmest month (BIO5), annual precipitation (BIO12), and the amount of cultivated fields
(AGRIC) and forested area (FOREST) (see Figure S2 and Table S4). Multiple-predictor
GLMs for bands 4, 8, 9, and 10 supported the significance of the predictors SEX and AGRIC
for
band 4 and
of the predictor COLOUR for band 10, but did not provide significant or
reliable results for bands 8 or 9 (see Table 3).
Based on the molecular weight of the corresponding bands in the SDS-PAGE profile
of the venom pool (Figure 1), we suspect that bands 1, 2, 3, 4, and 10 of the whole venom
SDS-PAGE profiles are likely to represent svMPs. Based on the same criterion, the content
of the remaining five bands is more difficult to identify, since they likely include different
toxin families (e.g., CRISP, LAAO, svSP). It is interesting to note that while the probability of
occurrence of bands 2 and 10 increases with the amount of forested area (FOREST) and with
the increase in the maximum temperature of the warmest month (BIO5), respectively, the
probability of occurrence of band 4 is negatively correlated with the amount of cultivated
fields (AGRIC; see Figures 4and S2). Snake venom metalloproteinases are thought to allow
fast prey subjugation [
27
,
67
,
68
] and possibly aid prey digestion (controversial; see [
87
]). In
a number of snake species undergoing an ontogenetic dietary shift from an ectotherm-based
to an endotherm-based diet, smaller, younger individuals have been shown to produce
more svMPs than larger, older specimens (e.g., V. latastei [
16
]; Bothrops asper [
17
]; Crotalus
viridis [
18
]). While the number of ectotherms consumed by V. seoanei decreases as the
viper grows, adult Iberian adders were found to consume more reptiles and amphibians in
warm habitats and forests, respectively, and to feed on small mammals more frequently
in agricultural areas [50]. In this scenario, the results of the single-predictor GLMs appear
to at least partially support the importance of svMPs in subduing amphibians and rep-
tiles. Indeed, the probability of occurrence of svMP-related bands 2 and 10 appears to be
positively correlated with conditions that might favour an increase in the consumption of
ectotherm prey (i.e., BIO5, FOREST). Conversely, the probability of detecting svMP-related
band 4 decreases with the increase in conditions that have been linked to an increase in
the consumption of endotherm prey by V. seoanei (i.e., AGRIC). Therefore, it could be
interesting to investigate the toxins actually comprising these SDS-PAGE whole venom
profile bands and test whether their presence provides any functional advantages in the
subjugation of one prey type rather than the other (e.g., [88,89]).
4. Conclusions
This study provides the first proteomic characterisation of the venom of Vipera seoanei.
Our results show that it presents the typical compositional pattern of Viperinae venoms,
with the five toxin families PLA
2
, svSP, DI, snaclec, and svMP being predominant. Surpris-
ingly, we found little variation among the SDS-PAGE profiles of the 49 considered venoms,
which appeared instead to be remarkably similar. By performing regression analyses, we
discovered significant effects of body size, colour phenotype, sex, locality of origin of the
vipers, and percentage of forested area on the low levels of venom variation detected. In
light of the highly dynamic scenario of snake venom variation, the information gathered
here introduces new questions as to the selective drivers that underlie venom evolution
in the Iberian adder. The development of genomic and functional studies could help us
understand if the similarities we detected across the analysed venoms are the result of
factors such as balancing selection and stabilising selection, or the by-product of a recent
population expansion. Finally, proteomic analyses developed on a finer scale could likely
unravel patterns of individual venom variation (e.g., sex-related, age-related, etc.) that
remain undiscovered in V. seoanei.
Toxins 2023,15, 371 12 of 22
5. Materials and Methods
5.1. Sampling
Between 2018 and 2021, a total of 49 individuals of V. seoanei were collected from
20 localities
distributed across the species’ range in the Iberian Peninsula, with a maximum
of five per locality (Figure 5; Table S5). Venom samples were collected from each individual
following the protocol reported by Avella et al. (2022) [16]. After venom extraction, tissue
samples (i.e., buccal swabs) to be used for genetic analyses (see Section 5.4) were collected,
and the sex and snout–vent length (SVL) of each snake were recorded. All vipers had an
SVL > 325 mm, and were thus considered adults (see [46,90]).
Figure 5.
Map of the 20 localities sampled across the Iberian distribution of V. seoanei. Black stars
indicate the sampled localities. The area in green represents the species’ full distributional range.
Additional details about the sampling are reported in Table S5.
Information about the colour phenotypes displayed by the collected vipers was also
recorded, and each individual was assigned to one of the five categories reported by
Martínez-Freiría et al. (2017) [
48
] (see Table S5 and Figure S3). Each viper was then
released exactly where it had been captured. Venoms were lyophilised in a Scanvac
(Coolsafe, Lynge, Denmark) freeze dryer and stored at
−
20
◦
C until being transported
to the Süssmuth Laboratory of the Institut für Chemie, Technische Universität Berlin
(Germany) for proteomic analysis.
Vipers, venoms, and tissue samples were collected with the permission of the Instituto
da Conservação da Natureza e das Florestas, Portugal (ref. 537/2018, 362/2019, 295/2020
and 146/2021), Xunta de Galicia, Spain (ref. EB-017/2019, 018/2020 and 015/2021), Gobierno
del Principado de Asturias, Spain (ref. 2019/003003, 2020/682020 and CO/09/017/2020),
and Junta de Castilla y León, Spain (ref. EP/CyL/56/2018, 27/2019, 192/2020 and 54/2021).
5.2. Bottom-Up Venomics
5.2.1. Venom Fractionation through RP-HPLC
Due to the low amount of venom we collected from each viper, we performed the
proteomic analyses on a pool composed of equal amounts (i.e., 10
µ
L) of each of the
49 individual
V. seoanei venoms. For reverse-phase chromatography (RP-HPLC), 1 mg of
lyophilised venom from the pool was dissolved to a final concentration of 10 mg/mL in
aqueous 5% (v/v) ACN with 1% (v/v) formic acid (HFo) and centrifuged for 5 min at
10,000
×
g. The supernatant was then fractionated using an HPLC Agilent 1200 (Agilent
Toxins 2023,15, 371 13 of 22
Technologies, Waldbronn, Germany) chromatography system equipped with a reversed-
phase Supelco Discovery BIO wide Pore C18-3 (4.6
×
150 mm, 3
µ
m particle size) column.
The following gradient with ultrapure water with 0.1% (v/v) HFo (solvent A) and ACN
with 0.1% (v/v) HFo (solvent B) was used at 1 mL/min, with a linear solvent change given
at min (B%): 0 (5%), 5 (5%), 100 (40%), 120 (70%), 130 (70%), and 5 min re-equilibration to
5% B. For monitoring the chromatography run, a diode array detector (DAD) was used at a
λ
= 214 nm detection wavelength. Samples were collected through time-based fractionation
(1 fraction/min) and peak fractions were dried in a centrifugal vacuum evaporator.
5.2.2. SDS-PAGE Profiling
The dried fractions were redissolved in 10
µ
L reducing 2
×
SDS sample buffer, heated
for 10 min at 95
◦
C, and separated using 12% SDS-PAGE (SurePage Bis-Tris, Genscript,
Piscataway, NJ, USA) run with MES buffer at 200 V for 21 min. A PageRuler Unstained
Protein Ladder (Thermo Scientific, Waltham, MA, USA) was used as the protein standard.
Gels were short-washed with water three times. Proteins were fixed three times for 10 min
each with hot fixation buffer (aqueous, 40% (v/v) methanol, 10% (v/v) acetic acid), stained
for 45 min in hot fast staining buffer (aqueous, 0.3% (v/v) HCl 37%, 100 mg/L Coomassie
250G) under constant mild shaking, and kept overnight at 4
◦
C in storage buffer (aqueous,
20% (v/v) methanol, 10% (v/v) acetic acid) for destaining. The produced gels were then
scanned for documentation and quantification.
To produce profiles that could allow the assessment of similarities and differences
among the 49 individual venoms, 20
µ
g of each lyophilised venom sample was loaded in
10 µL reducing 2×SDS sample buffer and subjected to SDS-PAGE profiling following the
same protocol applied for the venom pool. The resulting gels were scanned for documenta-
tion, and the obtained digital images were used for statistical analysis. SDS-PAGE profiling
was performed once for each venom sample.
5.2.3. Tryptic Digestion
The bands of interest of the SDS-PAGE venom pool profile were cut, dried with
500 µL
ACN, and stored at
−
20
◦
C until tryptic digestion. Disulphide bridges were reduced
with 30
µ
L freshly prepared dithiothreitol DTT (100 mM in 100 mM ammonium hydrogen
carbonate (ABC) per gel band) for 30 min at 56
◦
C and dried with 500
µ
L ACN for 10 min.
Cysteines were alkylated with freshly prepared iodacetamid IAC (55 mM in 100 mM ABC)
for 20 min at RT in the dark to protect the reduced thiols from oxidation and washed with
500
µ
L ACN for 2 min. Samples were dried with 500
µ
L ACN for 15 min, followed by
30 min incubation on ice with 20–30
µ
L freshly activated trypsin (13.3 ng/
µ
L, 10% (v/v)
ACN in 10 mM ABC; Thermo, Rockfeld, IL, USA). When necessary, additional volumes
of trypsin were added. Samples were incubated for 90 min on ice, then 20
µ
L ABC buffer
(10 mM) was added to all of them, and they were incubated overnight at 37
◦
C. Peptides
were extracted with 100
µ
L elution buffer (aqueous, 30% (v/v) ACN MS grade, 5% (v/v)
HFo) pre-warmed at 37
◦
C for 30 min. The supernatant was transferred into a separate
microtube and vacuum-dried. Following a second HPLC purification of 1 mg crude venom,
smaller dried fractions were submitted to LC-MS for direct peptide detection, without any
SDS-PAGE separation or tryptic digestion.
5.2.4. Mass Spectrometry
For the mass spectrometry (MS) analysis, the excised SDS-PAGE bands of interest
were re-dissolved in 30
µ
L aqueous 3% (v/v) ACN with 1% (v/v) HFo, and 20
µ
L of
each was injected into an Orbitrap XL mass spectrometer (Thermo, Bremen, Germany)
via an Agilent 1260 HPLC system (Agilent Technologies, Waldbronn, Germany) using a
reversed-phase Grace Vydac 218MS C18 (2.1
×
150 mm; particle size, 5
µ
m) column. The
following gradient with ultrapure water with 0.1% (v/v) HFo (solvent A) and ACN with
0.1% (v/v) HFo (solvent B) was used at 0.3 mL/min, with a linear buffer change given at
min (B%): 0 (5%), 1 (5%), 11 (40%), 12 (99%), 13 (99%), and 2 min re-equilibration to 5% B.
Toxins 2023,15, 371 14 of 22
The parameters in the ESI positive modus were as follows: 270
◦
C capillary temperature,
45 L/min sheath gas, 10 L/min auxiliary gas, 4.0 kV source voltage, 100.0
µ
A source
current, 20 V capillary voltage, 130 V tube lens. FTMS measurements were performed
with 1
µ
scans and 1000 ms maximal fill time. AGC targets were set to 10
6
for full scans
and to 3
×
10
5
for MS2 scans. MS2 scans were performed with a mass resolution (R) of
60,000 (at m/z400) for m/z250–2000. MS2 spectra were obtained in data-dependent
acquisition (DDA) mode as top2 with 35 V normalized CID energy, and with 500 as the
minimal signal required with an isolation width of 3.0. The default charge state was set
to z = 2, and the activation time to 30 ms. Unassigned charge states and charge state
1 were rejected.
LC-MS/MS data RAW files were converted into the MASCOT generic file (MGF)
format using MSConvert (Version 3.0.22187) with peak picking (vendor msLevel = 1
−
) [
91
].
For an automated database comparison, files were analysed using pFind Studio [
92
], with
pFind (Version 3.1.5) and the integrated pBuild. The parameters used were as follows:
MS Data (format: MGF; MS instrument: CID-FTMS); identification with Database search
(enzyme: Trypsin KR_C, full specific up to 3 missed cleavages; precursor tolerance
+20 ppm
;
fragment tolerance +20 ppm); open search setup with fixed carbamidomethyl [C] and Result
Filter (show spectra with FDR
≤
1%, peptide mass 500–10,000 Da, peptide length 5–100,
and show proteins with number of peptides > 1 and FDR ≤1%).
The used databases included UniProt ‘Serpentes’ (ID 8750, reviewed, canonical and
isoform, 2674 entries, last accessed on 10 February 2022; available at: https://www.uniprot.
org/) and the common Repository of Adventitious Proteins (215 entries, last accessed on
10 February 2022
; available at: https://www.thegpm.org/crap/index.html). The results
were batch-exported as PSM score of all peptides identified with pBuild and manually
cleared from decoy entries, contaminations, and artifacts. Finally, a list of unique peptide
sequences per sample with the best final score was generated. For a second confirmation of
identified sequences, all unique ones were analysed using BLAST search [
93
], with blastp
against the non-redundant protein sequences (nr) of the Uniprot database “Serpentes”
(taxid: 8570). In case of non-automatically annotated bands, files were checked manually
using Thermo Xcalibur Qual Browser (version 2.2 SP1.4), de novo annotated, and/or com-
pared on MS1 and MS2 levels with other bands to confirm band and peptide identities.
Deconvolution of isotopically resolved spectra was carried out by using the XTRACT
algorithm of Thermo Xcalibur.
5.2.5. Relative Quantification of the Venom Pool Proteome
The quantification protocol is adapted from the three-step hierarchical venom pro-
teome quantification protocol developed at the Evolutionary and Translational Venomics
Laboratory of the Institute of Biomedicine of Valencia [
94
–
96
], based on a combination of
the separation and quantification of HPLC peaks, SDS-PAGE bands, and ion intensities.
Briefly, it defines the normalised toxin abundances within a single SDS-PAGE band, with
the normalised values of the RP-HPLC peak integral measured at 214 nm, the gel band
intensity, and, if necessary, the MS ion intensity of the most abundant peptides identified.
The gel band abundance measurement was performed through densitometry on the
scanned gel images. Non-highly compressed PNG images of the SDS-PAGE gels were
processed using the software Fiji [
97
]. Colour depth was changed to 8 bit grayscale, and
the area and integrated density of each SDS-PAGE band were measured. The integrated
relative densities were calculated using a representative background region as reference
and normalised considering the area of the corresponding chromatographic peaks. In case
of multiple toxin families identified within a single band, normalised toxin abundances
were estimated based on the difference between the sum of the relative ion intensities of the
three most abundant peptide ions of a toxin family and the sum of the three most abundant
peptide ions of any other comigrated toxin family. Mass spectrometry proteomics data
have been deposited with the ProteomeXchange Consortium10 via the MassIVE partner
Toxins 2023,15, 371 15 of 22
repository under the project name “Snake venomics of Vipera seoanei” with the dataset
identifier “MSV000091535”.
5.3. Non-Metric Multidimensional Scaling
To assess the presence and extent of variation among the 49 V. seoanei venoms collected,
we applied the individual-based approach reported by Zancolli et al. (2017) [
98
]. Thus, we
generated a presence–absence matrix of the bands present in the individual SDS-PAGE
venom profiles, excluding bands with frequency = 1 (i.e., not informative) or difficult to
identify. We preferred analysing SDS-PAGE profiles rather than RP-HPLC profiles because
in the latter, reliable peak identification can often be difficult due to venom complexity
and differences in protein elution times. Indeed, during different RP-HPLC runs, the same
proteins can present different elution times for a number of reasons (e.g., fluctuations in
room temperature), further complicating the analysis of the chromatograms [98].
The final binary matrix (Table S3) was used to analyse patterns of venom variation us-
ing non-metric multidimensional scaling (NMDS) based on pairwise Bray–Curtis similarity
distances among profiles (see [
99
]). In order to keep the stress value below the 0.2 threshold
(i.e., poor fit [
100
]), we opted for a three-dimensional NMDS analysis. We then used the
individual NMDS scores on the first two axes (i.e., NMDS1 and NMDS2) to produce an
ordination plot representing the dissimilarity between the individual SDS-PAGE profiles of
the V. seoanei venoms analysed.
5.4. Predictors
To investigate which factors could potentially influence the occurrence of variation
in the venoms of V. seoanei, we considered a total of 14 predictors: (i) three referring to
biological traits of the vipers; (ii) one corresponding to the locality of origin of each viper
(POPULATION); (iii) eight describing the climatic and habitat conditions of the geographic
position where each viper was collected; and (iv) two corresponding to the first two
principal components (GEN1 and GEN2) of a spatial principal component analysis (SPCA)
performed for the interpolated mtDNA genetic distances estimated within V. seoanei.
The predictors related to the biological traits of the vipers were snout–vent length (SVL,
in mm), sex (SEX), and colour phenotype (COLOUR) of each individual. The uniform colour
phenotype, displayed by only one of the sampled individuals (i.e., 20VS165), was excluded
from analyses involving the predictor COLOUR in order to avoid model performance
hindering due to small sample size. The predictors related to the bioclimatic and habitat
conditions of the geographic position where each viper was collected comprised four
bioclimatic variables (i.e., annual mean temperature, BIO1; maximum temperature of the
warmest month, BIO5; annual precipitation, BIO12; precipitation of driest month, BIO14)
and the percentage of ground cover for four habitat types (i.e., cultivated fields, AGRIC;
forest, FOREST; moors, MOOR; pastures and grasslands, PASTURE). These eight predictors,
used in previous works on the ecology of V. seoanei and other Vipera species (e.g., [
101
,
102
]),
have been found to influence the distribution and abundance of V. seoanei prey such as
amphibians, reptiles, and small mammals (e.g., [
103
,
104
]), as well as their frequency in
the species’ diet [
50
]. Values for bioclimatic predictors were extracted from raster layers
at a resolution of 30 arc seconds (~1 km) from WorldClim version 2.1 ([
105
]; available
at www.worldclim.org, accessed on 25 November 2022). Values for habitat types were
extracted from Corine Land Cover version 2020_20u1 (available at https://land.copernicus.
eu/pan-european/corine-land-cover/clc2018, accessed on 26 October 2022) after grouping
land cover categories describing similar structural habitat types and upscaling the resulting
rasters at a 2 km2resolution (see [50]).
In order to obtain information on the genetic structure within V. seoanei, to be used
as a predictor of venom variation in regression analyses, we estimated mtDNA genetic
distances within this species. While V. seoanei sequences are available from previous works
(see [
45
]), we aimed to obtain a better geographic coverage of the genetic diversity existing
across the species’ range and, thus, decided to also consider genetic data obtained from
Toxins 2023,15, 371 16 of 22
the individuals we collected. To this end, total genomic DNA was extracted from buccal
swabs from 20 vipers (1 specimen from each sampled locality; see Table S5), using a stan-
dard saline method. Two mitochondrial (mtDNA) gene fragments, NADH dehydrogenase
subunit 4 (ND4; 630 bp)
and cytochrome b (cyt b; 554 bp), were amplified through poly-
merase chain reaction (PCR) using primers VSnd4-F/VSnd4-R [
45
] and CB1/CB2 [
106
],
respectively. Laboratory procedures for DNA amplification and sequencing followed the
protocols used in [
45
]. In order to produce a geographically comprehensive genetic distance
matrix, we considered the sequences of 65 specimens, namely 20 retrieved from this study,
44 from GenBank, and one from a specimen included exclusively in the genetic analyses to
enhance geographic coverage (specimen code 20VS008). Sequences were manually aligned
and edited using Geneious v 4.8.5 [
107
]. Uncorrected p-distances between populations
were estimated for the resulting alignment using MEGA X [
108
]. For information about the
65 specimens used to produce the genetic distance matrix and GenBank sequence accession
numbers, see Table S6.
The genetic distances calculated in this way were then spatially interpolated using the
kriging method and summarised through spatial principal component analysis (SPCA).
The first two of the resulting principal components, hereafter named GEN1 and GEN2,
explained 54% and 40% of the variance, respectively (see Figure S4). GEN1 and GEN2
were used as predictors in the subsequent regression analyses, as a proxy of the genetic
structure of V. seoanei. Values of GEN1 and GEN2 were obtained in ArcGIS version 10.5 by
extracting information on the geographic position of each specimen from the original SPCA
raster. Interpolations and SPCA were performed in ArcGIS version 10.5 [
109
], following
the same procedure applied in [
45
]. The genetic distance matrix was calculated using the
ape package [110] in R, version 4.2.2 [111].
5.5. Regression Analysis
We applied linear regression models to investigate the relationship between the pre-
dictors considered and the individual scores on the first two axes of the NMDS analysis
performed on the polymorphic bands of the SDS-PAGE venom profiles. To this end, using
NMDS1 and NMDS2 as response variables, we performed univariate linear regression mod-
els for each of the 14 predictors. In order to separate the effects of the different predictors
tested, we then performed multiple regression models built using NMDS1 and NMDS2 as
response variables and including all the predictors that were shown to be significant in the
univariate linear regression models.
We then tested whether the presence or absence of the polymorphic bands identified
in the venom SDS-PAGE profiles was significantly correlated with the predictors consid-
ered. Thus, we used single-predictor binomial generalised linear models (GLMs) for each
predictor, considering the presence/absence of each band as a binomial dependent variable.
Following the same approach described above, we then built multiple-predictor GLMs
including all the predictors that were shown to be significant in the single-predictor GLMs.
In all models, all continuous predictors were scaled (i.e., mean = 0; SD = 1). The pres-
ence of spatial autocorrelation for all variables included in the models with more than one
predictor was tested by calculating Moran’s I [
112
]. Correlation between the predictors in-
cluded in the multiple regression models was generally low (<50%), but the predictors POP-
ULATION and COLOUR showed a strong association (Cramér’s
V = 0.899
; see [
113
]). Since
the number of levels of POPULATION (five times greater than COLOUR’s;
i.e., 20 vs. 4
)
was more likely to hinder model performance (see [
114
]), it was excluded from multiple
regression models and multiple-predictor GLMs including the predictor COLOUR.
We used the packages vegan [
115
] to perform the NMDS analysis, car to perform the
regression models [
116
], and ggeffects [
117
] to plot model predictions. All analyses were
performed in R, version 4.2.2 [108].
Supplementary Materials:
The following supporting information can be downloaded at: https:
//www.mdpi.com/article/10.3390/toxins15060371/s1, Table S1: Quantification of the V. seoanei
Toxins 2023,15, 371 17 of 22
venom pool proteome; Table S2: MS bottom-up annotation of the peptides detected in the V. seoanei
venom pool analysed; Table S3: Final binary matrix of the 10 polymorphic bands considered for
statistical testing; Table S4: Results of the single-predictor binomial GLMs; Table S5: List of the
V. seoanei specimens considered in the present work; Table S6: List of the 65 V. seoanei specimens
considered to produce the matrix of uncorrected p-distances estimated for cyt b and ND4 sequences;
Figure S1: Whole venom profiles of the 49 V. seoanei specimens pooled to produce the reference
proteome, under reducing conditions; Figure S2: Model predictions of occurrence of bands 1, 2, 4, 8,
9, and 10 in individual SDS-PAGE venom profiles in relation to the continuous predictors tested in
single-predictor GLMs; Figure S3: The five geographically structured colour phenotypes currently
recognised within V. seoanei; Figure S4: Geographic genetic variation in V. seoanei.
Author Contributions:
I.A.: Conceptualization, Methodology, Software, Formal analysis, Inves-
tigation, Resources, Data curation, Writing—Original draft, Writing—Review and editing, Visu-
alization, Project administration; M.D.: Methodology, Software, Formal analysis, Investigation,
Resources, Data curation, Writing—Review and editing, Visualization; I.F.: Methodology, Resources,
Software, Formal analysis, Investigation, Data curation, Writing—Review and editing, Visual-
ization;
W.W.: Conceptualization
, Methodology, Writing—Review and editing, Supervision; N.L.:
Methodology, Resources, Software, Formal analysis, Writing—Review and editing, Visualization;
Ó.Z.: Resources;
R.D.S.: Resources, Supervision, Writing—Review; F.M.-F.: Conceptualization,
Methodology, Software, Formal analysis, Investigation, Resources, Data curation, Writing—Review
and editing, Supervision, Project administration, Funding acquisition. All authors have read and
agreed to the published version of the manuscript.
Funding:
I.A., I.F., N.L. and F.M.-F. are supported by the FCT—Fundação para a Ciência e a Tecnolo-
gia, Portugal (ref. SFRH/BD/137797/2018 and COVID/BD/152929/2023, SFRH/BD/148514/2019,
2020.06609.BD, and DL57/2016/CP1440/CT0010, respectively). Fieldwork was partially supported
by the European Regional Development Fund through the COMPETE programme and by Portuguese
National Funds, within the scope of the project PTDC/BIA-EVL/28090/2017-POCI-01-0145-FEDER-
028090. Work co-funded by the project NORTE-01-0246-FEDER-000063, supported by the Norte
Portugal Regional Operational Programme (NORTE2020), under the PORTUGAL 2020 Partnership
Agreement, through the European Regional Development Fund (ERDF).
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement:
Mass spectrometry proteomics data have been deposited with the
ProteomeXchange Con-sortium10 via the MassIVE partner repository under the project name “Snake
venomics of Vipera seoanei” with the dataset identifier “MSV000091535”. Additional data presented in
this study are available in the Supplementary Materials section.
Acknowledgments:
The authors thank Luis Albero, Ismael Espasandín, David do Pereiro, Miguel
Puras, Adrián Suarez, and Rafael Vázquez for their help with the fieldwork. I.A. is grateful to José
Carlos Brito for providing old literature on V. seoanei venom, and to Prem Daswani Aguilar, Fulvio
Licata, Yuri Simone, and Pedro Tarroso for the helpful comments on the statistical workflow. Special
thanks are extended Juan J. Calvete and Libia Sanz for producing preliminary RP-HPLC profiles
of some V. seoanei venoms, and to Susana Casal Vicente (Faculty of Pharmacy, University of Porto,
Portugal) for providing the freeze dryer used to lyophilise the samples. The authors acknowledge
the financial support received from National funds through the FCT—Fundação para a Ciência e a
Tecnologia in the scope of the project UIDB/50027/2020.
Conflicts of Interest: The authors declare no conflict of interest.
References
1.
Casewell, N.R.; Wüster, W.; Vonk, F.J.; Harrison, R.A.; Fry, B.G. Complex cocktails: The evolutionary novelty of venoms. Trends
Ecol. Evol. 2013,28, 219–229. [CrossRef] [PubMed]
2.
Chan, Y.S.; Cheung, R.C.F.; Xia, L.; Wong, J.H.; Ng, T.B.; Chan, W.Y. Snake venom toxins: Toxicity and medicinal applications.
Appl. Microbiol. Biotechnol. 2016,100, 6165–6181. [CrossRef] [PubMed]
3.
Arbuckle, K. Evolutionary context of venom in animals. In Evolution of Venomous Animals and Their Toxins; Gopalakrishnokone, P.,
Malhotra, A., Eds.; Springer: Dordrecht, Germany, 2017; pp. 3–31.
Toxins 2023,15, 371 18 of 22
4.
Fry, B.G.; Vidal, N.; van der Weerd, L.; Kochva, E.; Renjifo, C. Evolution and diversification of the Toxicofera reptile venom
system. J. Proteom. 2009,72, 127–136. [CrossRef] [PubMed]
5.
Kerkkamp, H.M.I.; Casewell, N.R.; Vonk, F.J. Evolution of the snake venom delivery system. In Evolution of Venomous Animals and
Their Toxins; Gopalakrishnokone, P., Malhotra, A., Eds.; Springer: Dordrecht, Germany, 2017; pp. 303–316.
6. Tasoulis, T.; Isbister, G.K. A review and database of snake venom proteomes. Toxins 2017,9, 290. [CrossRef] [PubMed]
7.
Chippaux, J.P.; Williams, V.; White, J. Snake venom variability: Methods of study, results and interpretation. Toxicon
1991
,
29, 1279–1303. [CrossRef]
8.
Casewell, N.R.; Jackson, T.N.; Laustsen, A.H.; Sunagar, K. Causes and consequences of snake venom variation. Trends Pharmacol.
Sci. 2020,41, 570–581. [CrossRef]
9.
von Reumont, B.M.; Anderluh, G.; Antunes, A.; Ayvazyan, N.; Beis, D.; Caliskan, F.; Crnkovi´c, A.; Damm, M.; Dutertre, S.;
Ellgaard, L.; et al. Modern venomics—Current insights, novel methods, and future perspectives in biological and applied animal
venom research. GigaScience 2022,11, giac048. [CrossRef]
10.
Fry, B.G.; Casewell, N.R.; Wüster, W.; Vidal, N.; Young, B.; Jackson, T.N. The structural and functional diversification of the
Toxicofera reptile venom system. Toxicon 2012,60, 434–448. [CrossRef]
11.
Jackson, T.N.; Sunagar, K.; Undheim, E.A.; Koludarov, I.; Chan, A.H.; Sanders, K.; Ali, S.A.; Hendrikx, I.; Dunstan, N.; Fry, B.G.
Venom down under: Dynamic evolution of Australian elapid snake toxins. Toxins 2013,5, 2621–2655. [CrossRef]
12.
Jackson, T.N.; Koludarov, I.; Ali, S.A.; Dobson, J.; Zdenek, C.N.; Dashevsky, D.; Op den Brouw, B.; Masci, P.P.; Nouwens, A.;
Josh, P.; et al.
Rapid radiations and the race to redundancy: An investigation of the evolution of Australian elapid snake venoms.
Toxins 2016,8, 309. [CrossRef]
13.
Damm, M.; Hempel, B.F.; Süssmuth, R.D. Old World Vipers—A review about snake venom proteomics of Viperinae and their
variations. Toxins 2021,13, 427. [CrossRef] [PubMed]
14.
Pla, D.; Sanz, L.; Sasa, M.; Acevedo, M.E.; Dwyer, Q.; Durban, J.; Pérez, A.; Rodríguez, Y.; Lomonte, B.; Calvete, J.J. Proteomic
analysis of venom variability and ontogeny across the arboreal palm-pitvipers (genus Bothriechis). J. Proteom.
2017
,152, 1–12.
[CrossRef] [PubMed]
15.
Queiroz, G.P.; Pessoa, L.A.; Portaro, F.C.; Maria de Fátima, D.F.; Tambourgi, D.V. Interspecific variation in venom composition
and toxicity of Brazilian snakes from Bothrops genus. Toxicon 2008,52, 842–851. [CrossRef] [PubMed]
16.
Avella, I.; Calvete, J.J.; Sanz, L.; Wüster, W.; Licata, F.; Quesada-Bernat, S.; Rodríguez, Y.; Martínez-Freiría, F. Interpopulational
variation and ontogenetic shift in the venom composition of Lataste’s viper (Vipera latastei, Boscá1878) from northern Portugal. J.
Proteom. 2022,263, 104613. [CrossRef]
17.
Saldarriaga, M.M.; Otero, R.; Núñez, V.; Toro, M.F.; Dıaz, A.; Gutiérrez, J.M. Ontogenetic variability of Bothrops atrox and Bothrops
asper snake venoms from Colombia. Toxicon 2003,42, 405–411. [CrossRef] [PubMed]
18.
Saviola, A.J.; Pla, D.; Sanz, L.; Castoe, T.A.; Calvete, J.J.; Mackessy, S.P. Comparative venomics of the Prairie Rattlesnake (Crotalus
viridis viridis) from Colorado: Identification of a novel pattern of ontogenetic changes in venom composition and assessment of
the immunoreactivity of the commercial antivenom CroFab®.J. Proteom. 2015,121, 28–43. [CrossRef] [PubMed]
19.
Malina, T.; Krecsák, L.; Westerström, A.; Szemán-Nagy, G.; Gyémánt, G.; M-Hamvas, M.; Rowan, E.G.; Harvey, A.L.;
Warrell, D.A.
;
Pál, B.; et al. Individual variability of venom from the European adder (Vipera berus berus) from one locality in Eastern Hungary.
Toxicon 2017,135, 59–70. [CrossRef] [PubMed]
20.
Menezes, M.C.; Furtado, M.F.; Travaglia-Cardoso, S.R.; Camargo, A.C.; Serrano, S.M. Sex-based individual variation of snake
venom proteome among eighteen Bothrops jararaca siblings. Toxicon 2006,47, 304–312. [CrossRef]
21.
Simizo, A.; Kitano, E.S.; Sant’Anna, S.S.; Grego, K.F.; Tanaka-Azevedo, A.M.; Tashima, A.K. Comparative gender peptidomics of
Bothrops atrox venoms: Are there differences between them? J. Venom. Anim. Toxins Incl. Trop. Dis.
2020
,26, e20200055. [CrossRef]
22.
Boldrini-França, J.; Corrêa-Netto, C.; Silva, M.M.; Rodrigues, R.S.; De La Torre, P.; Pérez, A.; Soares, A.M.; Zingali, R.B.;
Nogueira, R.A.
; Rodrigues, V.M.; et al. Snake venomics and antivenomics of Crotalus durissus subspecies from Brazil: Assessment
of geographic variation and its implication on snakebite management. J. Proteom. 2010,73, 1758–1776. [CrossRef] [PubMed]
23.
Neri-Castro, E.; Lomonte, B.; del Carmen Gutiérrez, M.; Alagón, A.; Gutiérrez, J.M. Intraspecies variation in the venom of the
rattlesnake Crotalus simus from Mexico: Different expression of crotoxin results in highly variable toxicity in the venoms of three
subspecies. J. Proteom. 2013,87, 103–121. [CrossRef] [PubMed]
24.
Sunagar, K.; Undheim, E.A.; Scheib, H.; Gren, E.C.; Cochran, C.; Person, C.E.; Koludarov, I.; Kelln, W.; Hayes, W.K.;
King, G.F.; et al.
Intraspecific venom variation in the medically significant Southern Pacific Rattlesnake (Crotalus oreganus helleri): Biodiscovery,
clinical and evolutionary implications. J. Proteom. 2014,99, 68–83. [CrossRef]
25.
Zancolli, G.; Calvete, J.J.; Cardwell, M.D.; Greene, H.W.; Hayes, W.K.; Hegarty, M.J.; Herrmann, H.W.; Holycross, A.T.;
Lannutti, D.I.
; Mulley, J.F.; et al. When one phenotype is not enough: Divergent evolutionary trajectories govern venom
variation in a widespread rattlesnake species. Proc. R. Soc. B 2019,286, 20182735. [CrossRef] [PubMed]
26.
Barlow, A.; Pook, C.E.; Harrison, R.A.; Wüster, W. Coevolution of diet and prey-specific venom activity supports the role of
selection in snake venom evolution. Proc. R. Soc. B 2009,276, 2443–2449. [CrossRef] [PubMed]
27. Daltry, J.C.; Wüster, W.; Thorpe, R.S. Diet and snake venom evolution. Nature 1996,379, 537–540. [CrossRef]
28.
Holding, M.L.; Strickland, J.L.; Rautsaw, R.M.; Hofmann, E.P.; Mason, A.J.; Hogan, M.P.; Nystrom, G.S.; Ellsworth, S.A.;
Colston, T.J.
; Borja, M.; et al. Phylogenetically diverse diets favor more complex venoms in North American pitvipers. Proc. Natl.
Acad. Sci. USA 2021,118, e2015579118. [CrossRef] [PubMed]
Toxins 2023,15, 371 19 of 22
29.
Creer, S.; Malhotra, A.; Thorpe, R.S.; Stöcklin, R.S.; Favreau, P.S.; Hao Chou, W.S. Genetic and ecological correlates of intraspecific
variation in pitviper venom composition detected using matrix-assisted laser desorption time-of-flight mass spectrometry
(MALDI-TOF-MS) and isoelectric focusing. J. Mol. Evol. 2003,56, 317–329. [CrossRef]
30.
Barua, A.; Mikheyev, A.S. Toxin expression in snake venom evolves rapidly with constant shifts in evolutionary rates. Proc. R.
Soc. B 2020,287, 20200613. [CrossRef]
31.
Casewell, N.R.; Wagstaff, S.C.; Harrison, R.A.; Renjifo, C.; Wüster, W. Domain loss facilitates accelerated evolution and neofunc-
tionalization of duplicate snake venom metalloproteinase toxin genes. Mol. Biol. Evol. 2011,28, 2637–2649. [CrossRef]
32.
Holding, M.L.; Margres, M.J.; Rokyta, D.R.; Gibbs, H.L. Local prey community composition and genetic distance predict venom
divergence among populations of the northern Pacific rattlesnake (Crotalus oreganus). J. Evol. Biol.
2018
,31, 1513–1528. [CrossRef]
33.
Saviola, A.J.; Gandara, A.J.; Bryson, R.W., Jr.; Mackessy, S.P. Venom phenotypes of the Rock Rattlesnake (Crotalus lepidus) and the
Ridge-nosed Rattlesnake (Crotalus willardi) from México and the United States. Toxicon
2017
,138, 119–129. [CrossRef] [PubMed]
34.
Sousa, L.F.; Portes-Junior, J.A.; Nicolau, C.A.; Bernardoni, J.L.; Nishiyama, M.Y., Jr.; Amazonas, D.R.; Freitas-de-Sousa, L.A.;
Mourão, R.H.V.; Chalkidis, H.M.; Valente, R.H.; et al. Functional proteomic analyses of Bothrops atrox venom reveals phenotypes
associated with habitat variation in the Amazon. J. Proteom. 2017,159, 32–46. [CrossRef]
35.
Massey, D.J.; Calvete, J.J.; Sánchez, E.E.; Sanz, L.; Richards, K.; Curtis, R.; Boesen, K. Venom variability and envenoming severity
outcomes of the Crotalus scutulatus scutulatus (Mojave rattlesnake) from Southern Arizona. J. Proteom.
2012
,75, 2576–2587.
[CrossRef] [PubMed]
36.
Schield, D.R.; Adams, R.H.; Card, D.C.; Corbin, A.B.; Jezkova, T.; Hales, N.R.; Meik, J.M.; Perry, B.W.; Spencer, C.L.;
Smith, L.L.; et al.
Cryptic genetic diversity, population structure, and gene flow in the Mojave rattlesnake (Crotalus scutulatus).
Mol. Phylogenet. Evol. 2018,127, 669–681. [CrossRef] [PubMed]
37.
Strickland, J.L.; Smith, C.F.; Mason, A.J.; Schield, D.R.; Borja, M.; Castañeda-Gaytán, G.; Spencer, C.L.; Smith, L.L.; Trápaga, A.;
Bouzid, N.M.; et al. Evidence for divergent patterns of local selection driving venom variation in Mojave Rattlesnakes (Crotalus
scutulatus). Sci. Rep. 2018,8, 17622. [CrossRef]
38.
Ferquel, E.; De Haro, L.; Jan, V.; Guillemin, I.; Jourdain, S.; Teynié, A.; d’Alayer, J.; Choumet, V. Reappraisal of Vipera aspis venom
neurotoxicity. PLoS ONE 2007,2, e1194. [CrossRef]
39.
Malina, T.; Babocsay, G.; Krecsák, L.; Erdész, C. Further clinical evidence for the existence of neurotoxicity in a population of the
European adder (Vipera berus berus) in eastern Hungary: Second authenticated case. Wilderness Environ. Med.
2013
,24, 378–383.
[CrossRef]
40.
Di Nicola, M.R.; Pontara, A.; Kass, G.E.N.; Kramer, N.I.; Avella, I.; Pampena, R.; Mercuri, S.R.; Dorne, J.L.M.; Paolino, G. Vipers
of major clinical relevance in Europe: Taxonomy, venom composition, toxicology and clinical management of human bites.
Toxicology 2021,453, 152724. [CrossRef]
41.
Paolino, G.; Di Nicola, M.R.; Pontara, A.; Didona, D.; Moliterni, E.; Mercuri, S.R.; Grano, M.; Borgianni, N.; Kumar, R.;
Pampena, R.
Vipera snakebite in Europe: A systematic review of a neglected disease. J. Eur. Acad. Dermatol. Venereol.
2020
,34, 2247–2260.
[CrossRef]
42.
World Health Organization (WHO). Snakebite Information and Data Platform. 2020. Available online: https://www.who.int/
teams/control-of-neglected-tropical-diseases/snakebite-envenoming/snakebite-information-and-data-platform/overview#
tab=tab_1,2020 (accessed on 19 January 2022).
43.
Freitas, I.; Ursenbacher, S.; Mebert, K.; Zinenko, O.; Schweiger, S.; Wüster, W.; Brito, J.C.; Crnobrnja-Isailovi´c, J.;
Halpern, B.;
Fahd, S.; et al.
Evaluating taxonomic inflation: Towards evidence-based species delimitation in Eurasian vipers (Serpentes:
Viperinae). Amphib-Reptil. 2020,41, 285–311. [CrossRef]
44.
Martínez-Freiría, F.; Brito, J.C. Vipera seoanei (Lataste, 1879). In Fauna Ibérica, Vol. 10, Reptiles, 2nd ed.; Ramos, M.A., Alba, J.,
Bellés, X., Gosálbez, J., Guerra, A., Macpherson, E., Serrano, J., Templado, J., Eds.; Museo Nacional de Ciencias Naturales, CSIC:
Madrid, Spain, 2014; pp. 942–957.
45.
Martínez-Freiría, F.; Velo-Antón, G.; Brito, J.C. Trapped by climate: Interglacial refuge and recent population expansion in the
endemic Iberian adder Vipera seoanei.Divers. Distrib. 2015,21, 331–344. [CrossRef]
46.
Lucchini, N.; Kaliontzopoulou, A.; Aguado-Val, G.; Martínez-Freiría, F. Sources of intraspecific morphological variation in Vipera
seoanei: Allometry, sex, and colour phenotype. Amphib.-Reptil. 2020,42, 1–16. [CrossRef]
47.
Martínez-Freiría, F.; Brito, J.C. Integrating classical and spatial multivariate analyses for assessing morphological variability in the
endemic Iberian viper Vipera seoanei.J. Zoolog. Syst. Evol. Res. 2013,51, 122–131. [CrossRef]
48.
Martínez-Freiría, F.; i de Lanuza, G.P.; Pimenta, A.A.; Pinto, T.; Santos, X. Aposematism and crypsis are not enough to explain
dorsal polymorphism in the Iberian adder. Acta Oecol. 2017,85, 165–173. [CrossRef]
49.
Braña, F.; Bea, A.; Saint Girons, H. Composición de la dieta y ciclos de alimentación en Vipera seoanei Lataste, 1879. Variaciones en
relación con la edad y el ciclo reproductor. Munibe Cienc. Nat. 1988,40, 19–27.
50.
Espasandín, I.; Galán, P.; Martínez-Freiría, F. Sex, size and eco-geographic factors affect the feeding ecology of the Iberian adder,
Vipera seoanei.Amphib.-Reptil. 2022,1, 1–16. [CrossRef]
51. Galán, P. Segregación ecológica en una comunidad de ofidios. Doñana. Acta Vertebr. 1988,15, 59–78.
52.
Luiselli, L.M.; Agrimi, U. Composition and variation of the diet of Vipera aspis francisciredi in relation to age and reproductive
stage. Amphib.-Reptil. 1991,12, 137–144. [CrossRef]
Toxins 2023,15, 371 20 of 22
53.
Brito, J.C. Feeding ecology of Vipera latastei in northern Portugal ontogenetic shifts, prey size and seasonal variations. Herpetol. J.
2004,14, 13–19.
54.
Tomovi´c, L.; An ¯
delkovi´c, M.; Golubovi´c, A.; Arsovski, D.; Ajti´c, R.; Sterijovski, B.; Nikoli´c, S.; Crnobrnja-Isailovi´c, J.; Lakuši´c, M.;
Bonnet, X. Dwarf vipers on a small island: Body size, diet and fecundity correlates. Biol. J. Linn. 2022,137, 267–279. [CrossRef]
55.
Santos, X.; Pleguezuelos, J.M.; Brito, J.C.; Llorente, G.A.; Parellada, X.; Fahd, S. Prey availability drives geographic dietary
differences of a Mediterranean predator, the Lataste’s viper (Vipera latastei). Herpetol. J. 2008,18, 16–22.
56.
Detrait, J.; Bea, A.; Saint Girons, H.; Choumet, V. Les variations geographiques du venin de Vipera seoanei Lataste (1879). Bull. Soc.
Zool. Fr. 1990,115, 277–285.
57.
Archundia, I.G.; de Roodt, A.R.; Ramos-Cerrillo, B.; Chippaux, J.P.; Olguín-Pérez, L.; Alagón, A.; Stock, R.P. Neutralization of
Vipera and Macrovipera venoms by two experimental polyvalent antisera: A study of paraspecificity. Toxicon
2011
,57, 1049–1056.
[CrossRef] [PubMed]
58.
Olaoba, O.T.; Dos Santos, P.K.; Selistre-de-Araujo, H.S.; de Souza, D.H.F. Snake venom metalloproteinases (SVMPs): A structure-
function update. Toxicon X 2020,7, 100052. [CrossRef]
59.
Požek, K.; Leonardi, A.; Pungerˇcar, J.; Rao, W.; Gao, Z.; Liu, S.; Laustsen, A.H.; Bakija, A.T.; Reberšek, K.; Podgornik, H.; et al.
Genomic Confirmation of the P-IIIe Subclass of Snake Venom Metalloproteinases and Characterisation of Its First Member, a
Disintegrin-Like/Cysteine-Rich Protein. Toxins 2022,14, 232. [CrossRef] [PubMed]
60.
Eble, J.A. Structurally robust and functionally highly versatile—C-type lectin (-related) proteins in snake venoms. Toxins
2019
,
11, 136. [CrossRef]
61.
Latinovi´c, Z.; Leonardi, A.; Šribar, J.; Sajevic, T.; Žužek, M.C.; Frangež, R.; Halassy, B.; Trampuš-Bakija, A.; Pungerˇcar, J.;
Križaj, I.
Venomics of Vipera berus berus to explain differences in pathology elicited by Vipera ammodytes ammodytes envenomation:
Therapeutic implications. J. Proteom. 2016,146, 34–47. [CrossRef]
62.
Gutiérrez, J.M.; Calvete, J.J.; Habib, A.G.; Harrison, R.A.; Williams, D.J.; Warrell, D.A. Snakebite envenoming. Nat. Rev. Dis.
Primers 2017,3, 17063. [CrossRef]
63.
Warrell, D.A. Epidemiology, clinical features and management of snake bites in Central and South America. In Venomous Reptiles
of the Western Hemisphere, 2nd ed.; Campbell, J.R., Lamar, W.W., Eds.; Cornell University Press: Ithaca, NY, USA, 2004; pp. 709–761.
64.
Kini, R.M. Serine proteases affecting blood coagulation and fibrinolysis from snake venoms. Pathophysiol. Haemos. Thromb.
2005
,
34, 200–204. [CrossRef]
65.
Lazarovici, P.; Marcinkiewicz, C.; Lelkes, P.I. From snake venom’s disintegrins and C-type lectins to anti-platelet drugs. Toxins
2019,11, 303. [CrossRef]
66.
Marcinkiewicz, C. Functional characteristic of snake venom disintegrins: Potential therapeutic implication. Curr. Pharm. Des.
2005,11, 815–827. [CrossRef]
67.
Gutiérrez, J.M.; Rucavado, A.; Escalante, T. Snake venom metalloproteinases: Biological roles and participation in the pathophysi-
ology of envenomation. In Evolution of Venomous Animals and their Toxins, 1st ed.; Mackessy, S.P., Ed.; CRC Press: Boca Ratón, FL,
USA, 2017; pp. 115–130.
68.
Ramos, O.H.P.; Selistre-de-Araujo, H.S. Snake venom metalloproteases—Structure and function of catalytic and disintegrin
domains. Comp. Biochem. Physiol. C Toxicol. Pharmacol. 2006,142, 328–346. [CrossRef] [PubMed]
69.
Jan, V.M.; Guillemin, I.; Robbe-Vincent, A.; Choumet, V. Phospholipase A2 diversity and polymorphism in European viper
venoms: Paradoxical molecular evolution in Viperinae. Toxicon 2007,50, 1140–1161. [CrossRef]
70.
Rouault, M.; Rash, L.D.; Escoubas, P.; Boilard, E.; Bollinger, J.; Lomonte, B.; Maurin, T.; Guillaume, C.; Canaan, S.;
Deregnaucourt, C.; et al.
Neurotoxicity and other pharmacological activities of the snake venom phospholipase A2 OS2: The
N-terminal region is more important than enzymatic activity. Biochemistry 2006,45, 5800–5816. [CrossRef] [PubMed]
71. Du, X.Y.; Clemetson, K.J. Snake venom L-amino acid oxidases. Toxicon 2002,40, 659–665. [CrossRef]
72.
Morjen, M.; Honoré, S.; Bazaa, A.; Abdelkafi-Koubaa, Z.; Ellafi, A.; Mabrouk, K.; Kovacic, H.; El Ayeb, M.; Marrakchi, N.;
Luis, J.
PIVL, a snake venom Kunitz-type serine protease inhibitor, inhibits
in vitro
and
in vivo
angiogenesis. Microvasc. Res.
2014
,
95, 149–156. [CrossRef] [PubMed]
73.
Mukherjee, A.K.; Mackessy, S.P.; Dutta, S. Characterization of a Kunitz-type protease inhibitor peptide (Rusvikunin) purified
from Daboia russelii russelii venom. Int. J. Biol. Macromol. 2014,67, 154–162. [CrossRef]
74.
Ferreira, I.G.; Pucca, M.B.; de Oliveira, I.S.; Cerni, F.A.; da Silva Jacob, B.D.C.; Arantes, E.C. Snake venom vascular endothelial
growth factors (svVEGFs): Unravelling their molecular structure, functions, and research potential. Cytokine Growth Factor Rev.
2021,60, 133–143. [CrossRef] [PubMed]
75.
Lecht, S.; Chiaverelli, R.A.; Gerstenhaber, J.; Calvete, J.J.; Lazarovici, P.; Casewell, N.R.; Harrison, R.A.; Lelkes, P.I.; Marcinkiewicz,
C. Anti-angiogenic activities of snake venom CRISP isolated from Echis carinatus sochureki.Biochim. Biophys. Acta Gen. Subj.
2015
,
1850, 1169–1179. [CrossRef] [PubMed]
76.
Sunagar, K.; Fry, B.G.; Jackson, T.N.; Casewell, N.R.; Undheim, E.A.; Vidal, N.; Ali, S.A.; King, G.F.; Vasudevan, K.;
Vasconcelos, V.; et al
. Molecular evolution of vertebrate neurotrophins: Co-option of the highly conserved nerve growth factor
gene into the advanced snake venom arsenal. PLoS ONE 2013,8, e81827. [CrossRef]
77.
Oliveira, A.L.; Viegas, M.F.; da Silva, S.L.; Soares, A.M.; Ramos, M.J.; Fernandes, P.A. The chemistry of snake venom and its
medicinal potential. Nat. Rev. Chem. 2022,6, 451–469. [CrossRef] [PubMed]
Toxins 2023,15, 371 21 of 22
78.
Hofmann, E.P.; Rautsaw, R.M.; Strickland, J.L.; Holding, M.L.; Hogan, M.P.; Mason, A.J.; Rokyta, D.R.; Parkinson, C.L. Compara-
tive venom-gland transcriptomics and venom proteomics of four Sidewinder Rattlesnake (Crotalus cerastes) lineages reveal little
differential expression despite individual variation. Sci. Rep. 2018,8, 15534. [CrossRef] [PubMed]
79.
Rautsaw, R.M.; Hofmann, E.P.; Margres, M.J.; Holding, M.L.; Strickland, J.L.; Mason, A.J.; Rokyta, D.R.; Parkinson, C.L.
Intraspecific sequence and gene expression variation contribute little to venom diversity in sidewinder rattlesnakes (Crotalus
cerastes). Proc. R. Soc. B 2019,286, 20190810. [CrossRef] [PubMed]
80.
Margres, M.J.; Patton, A.; Wray, K.P.; Hassinger, A.T.; Ward, M.J.; Lemmon, E.M.; Lemmon, A.R.; Rokyta, D.R. Tipping the scales:
The migration–selection balance leans toward selection in snake venoms. Mol. Biol. Evol. 2019,36, 271–282. [CrossRef]
81.
Margres, M.J.; McGivern, J.J.; Seavy, M.; Wray, K.P.; Facente, J.; Rokyta, D.R. Contrasting modes and tempos of venom expression
evolution in two snake species. Genetics 2015,199, 165–176. [CrossRef]
82.
Gibbs, H.L.; Sanz, L.; Pérez, A.; Ochoa, A.; Hassinger, A.T.; Holding, M.L.; Calvete, J.J. The molecular basis of venom resistance in
a rattlesnake-squirrel predator-prey system. Mol. Ecol. 2020,29, 2871–2888. [CrossRef]
83.
Poran, N.S.; Coss, R.G.; Benjamini, E.L.I. Resistance of California ground squirrels (Spermophilus beecheyi) to the venom of the
northern Pacific rattlesnake (Crotalus viridis oreganus): A study of adaptive variation. Toxicon 1987,25, 767–777. [CrossRef]
84.
van Thiel, J.; Khan, M.A.; Wouters, R.M.; Harris, R.J.; Casewell, N.R.; Fry, B.G.; Kini, R.M.; Mackessy, S.P.; Vonk, F.J.;
Wüster, W.; et al. Convergent evolution of toxin resistance in animals. Biol. Rev. 2022,97, 1823–1843. [CrossRef]
85.
Forsman, A. Variation in sexual size dimorphism and maximum body size among adder populations: Effects of prey size. J. Anim.
Ecol. 1991,60, 253–267. [CrossRef]
86.
Shine, R. Intersexual dietary divergence and the evolution of sexual dimorphism in snakes. Am. Nat.
1991
,138, 103–122.
[CrossRef]
87.
McCue, M.D. Prey envenomation does not improve digestive performance in western diamondback rattlesnakes (Crotalus atrox).
J. Exp. Zool. A Ecol. Gen. Physiol. 2007,307, 568–577. [CrossRef]
88.
Calvete, J.J.; Bonilla, F.; Granados-Martínez, S.; Sanz, L.; Lomonte, B.; Sasa, M. Venomics of the Duvernoy’s gland secretion of
the false coral snake Rhinobothryum bovallii (Andersson, 1916) and assessment of venom lethality towards synapsid and diapsid
animal models. J. Proteom. 2020,225, 103882. [CrossRef]
89.
Modahl, C.M.; Mrinalini; Frietze, S.; Mackessy, S.P. Adaptive evolution of distinct prey-specific toxin genes in rear-fanged snake
venom. Proc. R. Soc. B 2018,285, 20181003. [CrossRef] [PubMed]
90.
Braña, F. Vipera seoanei Lataste, 1879. In Fauna Ibérica Vol. 10, Reptiles; Salvador, A., Ed.; Museo Nacional de Ciencias Naturales,
CSIC: Madrid, Spain, 1998; pp. 489–497.
91.
Chambers, M.C.; Maclean, B.; Burke, R.; Amodei, D.; Ruderman, D.L.; Neumann, S.; Gatto, L.; Fischer, B.; Pratt, B.;
Egertson, J.; et al.
A cross-platform toolkit for mass spectrometry and proteomics. Nat. Biotechnol.
2012
,30, 918–920. [CrossRef]
[PubMed]
92.
Guangcan, S.; Yong, C.; Zhenlin, C.; Chao, L.; Shangtong, L.; Hao, C.; Meng-Qiu, D. How to use open-pFind in deep proteomics
data analysis?—A protocol for rigorous identification and quantitation of peptides and proteins from mass spectrometry data.
Biophys. Rep. 2021,7, 207–226. [CrossRef]
93.
Altschul, S.F.; Wootton, J.C.; Gertz, E.M.; Agarwala, R.; Morgulis, A.; Schäffer, A.A.; Yu, Y.K. Protein database searches using
compositionally adjusted substitution matrices. FEBS J. 2005,272, 5101–5109. [CrossRef]
94.
Calderón-Celis, F.; Cid-Barrio, L.; Encinar, J.R.; Sanz-Medel, A.; Calvete, J.J. Absolute venomics: Absolute quantification of intact
venom proteins through elemental mass spectrometry. J. Proteom. 2017,164, 33–42. [CrossRef]
95.
Calvete, J.J. Next-generation snake venomics: Protein-locus resolution through venom proteome decomplexation. Expert Rev.
Proteom. 2014,11, 315–329. [CrossRef] [PubMed]
96.
Eichberg, S.; Sanz, L.; Calvete, J.J.; Pla, D. Constructing comprehensive venom proteome reference maps for integrative venomics.
Expert Rev. Proteom. 2015,12, 557–573. [CrossRef]
97.
Schindelin, J.; Arganda-Carreras, I.; Frise, E.; Kaynig, V.; Longair, M.; Pietzsch, T.; Preibisch, S.; Rueden, C.; Saalfeld, S.;
Schmid, B.; et al. Fiji: An open-source platform for biological-image analysis. Nat. Methods 2012,9, 676–682. [CrossRef]
98.
Zancolli, G.; Sanz, L.; Calvete, J.J.; Wüster, W. Venom on-a-chip: A fast and efficient method for comparative venomics. Toxins
2017,9, 179. [CrossRef] [PubMed]
99.
Minchin, P.R. An evaluation of relative robustness of techniques for ecological orderings. Vegetatio
1987
,71, 145–156. [CrossRef]
100.
Kruskal, J.B. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika
1964
,29, 1–27.
[CrossRef]
101.
Martínez-Freiría, F.; Sillero, N.; Lizana, M.; Brito, J.C. GIS-based niche models identify environmental correlates sustaining a
contact zone between three species of European vipers. Divers. Distrib. 2008,14, 452–461. [CrossRef]
102.
Chamorro, D.; Martínez-Freiría, F.; Real, R.; Muñoz, A.R. Understanding parapatry: How do environment and competitive
interactions shape Iberian vipers’ distributions? J. Biogeogr. 2021,48, 1322–1335. [CrossRef]
103.
Mira, A.; Marques, C.C.; Santos, S.M.; Rosário, I.T.; Mathias, M.L. Environmental determinants of the distribution of the Cabrera
vole (Microtus cabrerae) in Portugal: Implications for conservation. Mamm. Biol. 2008,73, 102–110. [CrossRef]
104.
Sillero, N.; Brito, J.C.; Skidmore, A.K.; Toxopeus, A.G. Biogeographical patterns derived from remote sensing variables: The
amphibians and reptiles of the Iberian Peninsula. Amphib.-Reptil. 2009,30, 185–206. [CrossRef]
Toxins 2023,15, 371 22 of 22
105.
Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol.
2017
,
37, 4302–4315. [CrossRef]
106.
Palumbi, S.R. The polymerase chain reaction. In Molecular Systematics; Hillis, D., Moritz, C., Mable, B.K., Eds.; Sinauer Associates:
Sunderland, MA, USA, 1996; pp. 205–247.
107.
Kearse, M.; Moir, R.; Wilson, A.; Stones-Havas, S.; Cheung, M.; Sturrock, S.; Buxton, S.; Cooper, A.; Markowitz, S.; Duran, C.; et al.
Geneious basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data.
Bioinformatics 2012,28, 1647–1649. [CrossRef]
108.
Kumar, S.; Stecher, G.; Li, M.; Knyaz, C.; Tamura, K. MEGA X: Molecular Evolutionary Genetics Analysis across computing
platforms. Mol. Biol. Evol. 2018,35, 1547–1549. [CrossRef]
109. ESRI. ArcGIS Desktop: Release 10.5; Environmental Systems Research Institute: Redlands, CA, USA, 2016.
110.
Paradis, E.; Claude, J.; Strimmer, K. APE: Analyses of phylogenetics and evolution in R language. Bioinformatics
2004
,20, 289–290.
[CrossRef] [PubMed]
111.
R Core Team. R: A Language and Environment for Statistical Computing; R Core Team: Vienna, Austria, 2022; Available online:
https://www.R-project.org (accessed on 9 November 2022).
112. Moran, P.A. Notes on continuous stochastic phenomena. Biometrika 1950,37, 17–23. [CrossRef] [PubMed]
113. Cramér, H. Mathematical Methods of Statistics; Princeton University Press: Princeton, NJ, USA, 1999.
114. Hastie, T.; Tibshirani, R.; Friedman, J.H.; Friedman, J.H. The Elements of Statistical Learning: Data Mining, Inference, and Prediction;
Springer: New York, NY, USA, 2009.
115.
Oksanen, J. Design Decisions and Implementation Details in Vegan. Vignette of the Package Vegan. Available online: https:
//cran.r-project.org/web/packages/vegan/vignettes/decision-vegan.pdf (accessed on 9 November 2022).
116. Fox, J.; Weisberg, S. An R Companion to Applied Regression, 3rd ed; Sage Publications: Washington, DC, USA, 2018.
117.
Lüdecke, D. ggeffects: Tidy data frames of marginal effects from regression models. J. Open Source Softw.
2018
,3, 772. [CrossRef]
Disclaimer/Publisher’s Note:
The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.
