DOI: 10.5061/dryad.nzs7h44qs
title: Ecological and biogeographic processes drive the proteome evolution of snake venom

description: The emergence of venom is an evolutionary innovation that favored the
 diversification and survival of snakes. The composition of snake venoms is
 known in detail from venom gland proteomic data. However, there is still a
 gap of knowledge about the forces that lead to the expression of different
 toxins in different proportions in the venom cocktail across space and
 time. LOCATION World. TIME PERIOD Modern. MAJOR TAXA STUDIED Elapidae and
 Viperidae. We integrated proteomic data with phylogenetic comparative
 methods to understand how ecological and biogeographical processes drive
 the evolution of snake venom. We observed that more productive
 environments favor a more complex venom, which presents a higher diversity
 and similarity on the toxin proportions in its composition. We found that
 the taxa that live in islands, where there is lower variability of
 resources, tended to present  less complex venom dominated by few
 toxins. In this case, the island's isolation seems to be a relevant
 element for a faster fixation of specific venom compositions. We show that
 ecological and biogeographic processes, which can act differentially over
 time and space, affect the gene expression of toxins in snake venoms.
description: We
 obtained proteomic data of snake venom from Tasoulis and Isbister (2017).
 Proteomic data was available for 112 taxa, 74 of which belonging to the
 Viperidae family, and 38 to the Elapidae family. The data includes
 information on the proteomes of some subspecies and isolated populations
 on islands (Table S1). In addition, this dataset corresponds to ten types
 of toxins found in these species: phopolispases A2 (PLA2), serine
 proteases (SVSP), metalloproteases (SVMP), Kunitz peptides (KUN), 3-finger
 toxins (3FT), L-amino acid oxidases (LAAO), cysteine-rich secretory
 proteases (CRiSP), C-type lectins (CTL), disintegrins (DIS), and
 natriuretic peptides (NP).  The defensins (DEF), metalloprotease
 inhibitors (Mpi), vascular endothelial growth factors (VEGF), and
 cystatins (CYS) toxins were not considered for analysis due to their low
 frequency in proteomes. To understand how snake venom compositions vary,
 we performed a Pearson correlation analysis between the toxins from the
 venoms in R platform v.4.0.2 (R Core Team, 2020). Additionally, we
 estimated the complexity of the species’ venom using the standardized
 Levin index (Winemiller & Pianka, 1990). The Levin index is a
 measure that calculates the niche trophic breadth of an organism, where
 the number of items and the proportion of each item are used in the
 calculation. In this case, species that feed on a greater variety of items
 and with more equitable proportions have a broader trophic niche. However,
 this measure can be used to quantify the complexity of the venom, since we
 expect that venoms with a greater diversity of toxins and in more
 equitable proportions will have more diverse pathways of action against
 their prey. The Levin index is calculated as
 follows: Where<i> B</i>i is the Levin’s standardized index of the snake species i, <i>P</i>ij is the proportion of the toxin j in the species i; and n is the total number of toxins in the venom. Thus, Levin's index values oscillate between 0 and 1, where a value close to 0 represents a less complex venom, while values close to 1 are a highly complex venom. We collected georeferenced occurrences of each species from the GBIF platform (https://www.gbif.org/, associated references are in the Supporting Information), with the help of the “spocc” package (Chamberlain, 2020) in R v.4.0.2 (R Core Team, 2020). We also collected occurrence data of holotypes of all species that are not present in the GBIF database, using the Reptile Database (https://reptile-database.reptarium.cz/). All spatial coordinates were validated with the proposed distributions on the IUCN Red List to eliminate possible erroneous points. From the geographic coordinates, we obtained data of Average Annual Temperature (Bio1) and Precipitation Seasonality (Bio15) from WorldClim (Booth et al., 2014; https://www.worldclim.org/), Net Primary Productivity (NPP) from the FAO database (www.fao.org), and the Speed of Temperature Change in the last 21 thousand years (Velocity) by Sandel et al. (2011). For each of the variables (Bio1, Bio15, NPP, and Velocity), we estimate the average for each species. We used the variables Bio1 and NPP as proxies for the diversity of prey available in the environment (<i>Food availability hypothesis</i>) (Waide et al., 1999; Allen et al., 2002). We didn't include dietary data of Viperidae and Elapidae due to the lack of data for more than 50% of the analyzed taxa (for the majority of species there are few records of the dietary resources consumed by them) (Grundler, 2020). We used the variables Bio15 and Velocity as proxies for the resource variation over the year and historical time (Brito, 2004; Brown & Shine, 2007; Morales-Barbero et al., 2021), respectively (<i>Food seasonality hypothesis</i>). To assess the effect of body size on snake venoms (<i>Body </i><i>size hypothesis</i>), we obtained the body mass estimate in grams available from Feldman et al. (2016). Additionally, to assess the biogeographic effect on snake venoms (<i>Biogeography </i><i>hypothesis</i>), we collected for each proteome data point the locality of occurrence of the taxon to classify whether it occurs on islands or the mainland (Table S2). <i>Phylogenetic Reconstruction</i> We used the genetic alignments of Figueroa et al. (2016), which is based on five mitochondrial genes and five nuclear genes. Given that some species for which we had proteome data were not included in the study by Figueroa et al. (2016), we supplemented our data with CytB, ND2, ND4, 12S, and 16S sequences available in the GenBank database (Table S2). We used the species <i>Rena humilis</i>, <i>Boa constrictor</i>, and <i>Iguana iguana</i> as external groups. The final alignment of the sequences was performed with the MAFFT software (https://mafft.cbrc.jp/alignment/server/). We conducted a Bayesian phylogenetic reconstruction using the BEAST v 2.5 software (Bouckaert et al., 2019), using the evolutionary model of GTR + G + I, a relaxed lognormal molecular clock, and the Yule process prior.  For the calibration of the tree, we used five fossils associated with the first register of Elapidae, Viperinae, Crotalinae, Scolecophidia, and Boidae (see Table S3 for prior distribution and fossil sources).  Furthermore, we used 50 million generations of MCMC chains, and the chain convergence was evaluated using Tracer software v 1.7.1 (Rambaut et al., 2018), using a 10% burn-in and an ESS&gt; 200. Finally, we used TreeAnnotator v.2.6.3 (http://beast.community/treeannotator) to obtain the Major Credibility Tree. Taxa that lacked genetic data (Table S2) were added randomly within the genus, with the help of the “phytools” package (Revell, 2012) in the R v.4.0.2 environment (R Core Team, 2020). <i>Evolutionary Models</i> To evaluate the mode of evolution of venom complexity in the snake families Elapidae and Viperidae, we applied an evolutionary model-fitting approach to both families, together, as well as to each of the lineages. We evaluated different evolutionary models derived from the stochastic processes of Brownian motion (BM) (Felsenstein, 1985) and Ornstein-Uhlenbeck (OU) (Hansen & Martins, 1996) with the help of the “mvMORPH” package (Clavel et al., 2015). The BM model establishes that phenotypic differences are the result of a stochastic accumulation of changes, which occur independently and at a specific rate (σ2) along the lineages, increasing the trait variance through time (Felsenstein, 1985). The OU model is a modified version of the BM model that incorporates a selection process. Thus, the OU model presents adaptive zones that function as attractors of biological characteristics, representing evolution through adaptive processes (Butler & King, 2004). Thus, two parameters are incorporated in the model, one that describes the strength of the selection (α) while another describes an adaptive optimal (θ) (Hansen & Martins, 1996). When α=0, the OU model is reduced to a BM. The four evolutionary models that we evaluated were: (i) BM1 (one rate), where there is a single drift rate (σ2); (ii) BM2 (two rates), where there are different σ2 for species that inhabit islands and continents; (iii) OU1 (one optimal), where there is a single adaptive optimal (θ); and (iv) OU2 (two optima), where there are different adaptive peaks for species that inhabit islands and continents. Therefore, BM2 and OU2 models allowed us to validate the existence of different rates and optima, respectively, for the complexity of the venoms in islands and continents. To perform the BM (two rates) and OU (two optima) models, we estimate the ancestral states of the binary characteristic island/continent (1000 simulations) using the <i>make.simmap</i> in the “phytools” package (Revell, 2012). To take phylogenetic uncertainty into account, we generated 100 trees where we included the missing species (n = 7) randomly within the genus with the “phytools” package (Revell, 2012) in the R v.4.0.2 environment (R Core Team, 2020), as recommended by Martinez et al. (2015). Subsequently, we ran the model with these 100 trees and estimated the average of AIC values. We selected the model using the lowest AIC value (average AIC) and the highest AIC weight (wi) (average wi) (Wagenmakers & Farrell, 2004; Johnson & Omland, 2004). <i>Phylogenetic Path Analysis</i> To assess the hypotheses <i>Food availability</i>, <i>Food seasonality</i>, <i>Body size</i>, and <i>Biogeography</i>, we performed a Phylogenetic Path Analysis (PPA) (von Hardenberg & Gonzalez-Voyer, 2013). We built nine possible causal models, represented by acyclic graphs, which illustrate the relationship between our variables (Figure S1). The causal models were developed to consider the separated or combined effects of the predictors on toxins and on the venom complexity. These causal models illustrate different paths to describe the <i>Food availability hypothesis</i> (Figure S1a, b, f, h), the <i>Food seasonality hypothesis</i> (Fig. S1a, b, d, e, f, g, h, i), the <i>Body size hypothesis</i> (Fig. S1a, b, c) and the <i>Biogeography hypothesis</i> (Fig. S1a, b, i) (more details of the causal models in the legend of Fig. S1).  For each response variable, we select the best causal model based on the C statistic's Information Criterion (CICc), which is a modified version of the Akaike Information Criterion, AIC (von Hardenberg & Gonzalez-Voyer, 2013). Models with ΔCICc lower than two were considered to be the best fitting ones. The PPAs were built and tested with the support of the “phylopath” package (van der Bijl, 2018) in R v.4.0.2 (R Core Team, 2020). For the analysis, our venom complexity index was used on a logarithmic scale to comply with the assumption of normality. To take into account phylogenetic uncertainty, PPAs were performed using 100 phylogenetic trees. Subsequently, we ran the PPAs with these 100 trees and estimated the average, minimum, and maximum of the standardized coefficients and the 95% confidence intervals.