Snakes on a plain: biotic and abiotic factors determine venom compositional variation in a wide-ranging generalist rattlesnake
BMC Biology volume 21, Article number: 136 (2023)
Snake venoms are trophic adaptations that represent an ideal model to examine the evolutionary factors that shape polymorphic traits under strong natural selection. Venom compositional variation is substantial within and among venomous snake species. However, the forces shaping this phenotypic complexity, as well as the potential integrated roles of biotic and abiotic factors, have received little attention. Here, we investigate geographic variation in venom composition in a wide-ranging rattlesnake (Crotalus viridis viridis) and contextualize this variation by investigating dietary, phylogenetic, and environmental variables that covary with venom.
Using shotgun proteomics, venom biochemical profiling, and lethality assays, we identify 2 distinct divergent phenotypes that characterize major axes of venom variation in this species: a myotoxin-rich phenotype and a snake venom metalloprotease (SVMP)-rich phenotype. We find that dietary availability and temperature-related abiotic factors are correlated with geographic trends in venom composition.
Our findings highlight the potential for snake venoms to vary extensively within species, for this variation to be driven by biotic and abiotic factors, and for the importance of integrating biotic and abiotic variation for understanding complex trait evolution. Links between venom variation and variation in biotic and abiotic factors indicate that venom variation likely results from substantial geographic variation in selection regimes that determine the efficacy of venom phenotypes across populations and snake species. Our results highlight the cascading influence of abiotic factors on biotic factors that ultimately shape venom phenotype, providing evidence for a central role of local selection as a key driver of venom variation.
Understanding the processes that generate and maintain phenotypic diversity is central to understanding how natural selection shapes the evolution of traits and the variation in traits observed across populations. Predator–prey dynamics are a core feature of ecological systems, and the traits involved in these antagonistic interactions are under strong selective pressures that may give rise to novel, extreme, and complex phenotypes, which may also vary substantially across populations [1, 2] The resulting trophic adaptations that mediate predator–prey relationships represent a coalescence of morphological, behavioral, and physiological characteristics whose adaptive qualities are context-dependent and respond to a variety of influences . Elucidating the biotic and abiotic determinants that shape these traits requires an integrative approach that examines possible evolutionary influences, their relationship to one another, and their combined effects on the individual components that comprise complex phenotypes.
A variety of extreme trophic traits have evolved among snakes, including a chemical means of prey capture via the production, storage, and delivery of heterogeneous mixtures of toxins in venom-producing lineages. As such, snake venoms provided a rich system for the investigation of the evolutionary drivers and mechanisms shaping complex trophic adaptations due to their clear genotype–phenotype relationships and their potent measurable bioactive properties [4, 5]. Though evolutionary history tends to correspond with broad compositional trends in venom among venomous snake lineages [2, 6, 7], venoms may also vary substantially within species and populations [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. What factors drive venom variation, particularly within species, remains poorly understood [23, 24].
The association between diet and venom variation in some species points to prey availability as a major biotic determinant of venom composition and variation [22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37]. Indeed, the multifaceted influence of diet has been demonstrated by the expression of prey-specific toxins [25,26,27, 35, 37], correlations between toxinological diversity and prey diversity [2, 28], and the co-evolutionary dynamics of some populations to venom resistance in prey [23, 29, 38]. However, other studies have shown evidence that diet did NOT correlate with patterns of venom variation, indicating that other forces may also play major roles in shaping venom composition in different contexts [20, 39]. These varying conclusions indicate that the interplay of various determinants shaping venom composition is likely multifactorial, but only recently have venom studies explored the associations between biotic and abiotic factors and their relationship to venom phenotypes [19, 20]. Ultimately, the influences of population structure, gene flow, diet, and environment are complex and likely interacting, and it has remained unclear how these forces may synergistically contribute to shaping venom variation.
Here, we examine venom compositional diversity in a common, wide-ranging species and explore multiple potential drivers of variation in venom as observed across its range. The Prairie Rattlesnake (Crotalus viridis viridis) is a habitat generalist species found from northern Mexico, through the Great Plains of the western United States, to southern Canada (Fig. 1; . Its venom is known to induce hemorrhage and tissue degradation due to the abundance of larger enzymes such as snake venom metalloproteases (SVMPs) and snake venom serine proteases (SVSPs), and smaller nonenzymatic myotoxin peptides result in tetanic paralysis of prey as well as myonecrosis [41,42,43]. Venom of the Prairie Rattlesnake has also been shown to contain many other classes of enzymes that are characteristic of rattlesnake venoms in general, including phospholipases A2 (PLA2), L-amino acid oxidases (L-AAO), and phosphodiesterases (PDE) , and ontogenetic shifts in venom composition, transitioning from a higher SVMP-level phenotype as juveniles to a myotoxin-rich venom as adults, occur in some populations . Previous research on C. v. viridis venom has shown population-level differences in myotoxin presence  and PLA2 enzyme variation , but no studies have examined proteome-scale venom variation across the range of this species.
The current study uses a multipronged approach that examines phylogenetic and population genetic structure, diet, and environmental variables to provide a broad foundational context to understand factors that shape geographic patterns of venom variation. First, we characterize venom phenotypes range-wide with chromatographic profiling, shotgun proteomics, and biochemical and prey lethality assays. We then compare these patterns to range-wide phylogeographic structure, geographic trends in diet, and analyses of environmental niche variables and environmental niche modeling (ENM) to decipher the influences of biotic and abiotic factors that likely shape variation in venom composition.
In total, 200 venom samples were collected from nine states (Fig. 1; Additional file 1: Table S1). Thirteen snakes collected were neonates (6.5%), 18 snakes were juveniles (9%), 21 were subadults (10.5%), and 147 were adults (73.5%). Because C. v. viridis venom composition varies ontogenetically , only adults were analyzed for trends in venom variation.
High-performance liquid chromatography
HPLC analysis revealed complex toxin profiles from all snakes analyzed, and myotoxins eluted at approximately 22 min, SVSPs, cysteine-rich secretory proteins (CRiSPs), and PLA2s between 50 and 70 min, L-AAO at 78–79 min, and SVMPs between 85 and 90 min (Fig. 2a, b; based on Saviola et al., 2015  and SDS-PAGE/MS data on peaks). Considerable complexity of toxin composition was observed in peaks that eluted between 50 and 70 min, so these percentages were not used in analysis because quantification of specific toxin families was unclear. The proportion of myotoxin ranged from 2.2 to 60.7% (average = 32.2, median = 35.3), L-AAO percentages ranged from 0 to 2.9% (average = 0.96, median = 0.9), and SVMP proportion ranged from 3.5 to 50.1% (average = 19.8, median = 13.4); for all of these peaks, the stated toxin family was exceedingly the most abundant representative. Snakes from northern regions (ex. Montana, Wyoming, Nebraska) displayed higher percentages of myotoxin and lower proportions of SVMPs (Fig. 2a) than snakes from southern areas (e.g., Texas, New Mexico; Fig. 2b). Snakes from southern regions in general had lower percentages of myotoxin and higher percentages of SVMPs (Fig. 2b). All percentages determined by area under the curve of HPLC chromatograms were used for geographic venom analysis using Mantel tests and linear regression analysis.
The 147 adult venoms tested displayed high levels of variation in all enzymes. Adult SVMP specific activity ranged from 0.09 to 1.36 ΔA342nm/min/mg venom protein (average = 0.51, median = 0.39). Thrombin-like specific activity ranged from 346.1 to 2607.4 nmol product produced/min/mg venom (average = 1364.4, median = 1285.3) and kallikrein-like specific activity ranged from 64.6 to 2851.9 nmol product produced/min/mg venom (average = 1012.2, median = 962.1). Venoms ranged from 6.7 to 125.6 nmol product formed/min/mg venom (average = 42.7, median = 39.4) for L-amino acid oxidase activity. Phosphodiesterase activity ranged from 0.11 to 1.6 ΔA400nm/min/mg venom (average = 0.5, median = 0.46). Phospholipase A2 activity ranged from 4.7 to 72.3 nmol product/min/mg (average = 25.4, median = 22.3). These values were investigated for geographic patterns in activity with Mantel tests and linear regression analysis against latitude.
Proteins identified from shotgun proteomic analyses were binned into toxin protein families, and total intensity for each family was compared across samples and used for principal component analysis. We detected 94 proteins across 20 protein families’ total (Additional file 2: Tables S2-S3). The toxin families with the largest number of proteins assigned to them were SVSPs , SVMPs , PLA2s , C-type lectins (CTL; 9), disintegrins , and myotoxins . Venoms of snakes from New Mexico and southern Colorado had a higher abundance of SVMPs than those from northern Colorado and Montana, and the venom from the Montana snake had the highest amount of SVSPs. The individual from Socorro Co., NM, had the highest levels of PLA2, while the snakes from Luna Co., NM, had the lowest levels. CTLs were most abundant in NM snakes and least abundant in northern Colorado snakes. Disintegrin levels were highest in snakes from NM and southern Colorado. With the exception of one individual from Luna Co., snakes in New Mexico had the lowest levels of myotoxins. Myotoxin abundance was also lower in southern Colorado than northern Colorado or Montana. Principal component analysis revealed the presence of two distinct clusters that corresponded to a northern group containing MT and northern CO venoms and a southern group with southern CO and NM venoms (Fig. 2c). Between these groupings, we detected significantly higher abundances of cystatin (p < 0.001), L-AAO (p = 0.03), myotoxin (p = 0.03), PDE (p = 0.04), thrombin-like serine proteases (p = 0.04), VF (p = 0.01), and VNGF (p = 0.03) in northern snakes and higher levels of disintegrins (p < 0.001) and SVMP (p < 0.001) in southern snakes (Fig. 2d).
Based on a Mantel test, the percentage of myotoxin and percentage of SVMPs significantly varied based on geographic location (r = 0.323, p = 0.001; r = 0.316, p = 0.001, respectively), while the percentage of L-AAO did not vary with geographic location (r = 0.032, p = 0.232). Of other enzymes tested, metalloproteases (r = 0.28, p = 0.001), thrombin-like serine proteases (r = 0.17, p = 0.001), and kallikrein-like serine proteases (r = − 0.03, p = 0.016) showed significant geographic variation.
The percentage of SVMPs (Fig. 3a; R2 = 0.4412, p < 0.001) and the percentage of myotoxin (Fig. 3b; R2 = 0.408, p < 0.001) varied significantly with latitude in linear regression analysis. Strikingly, rather than an increasing or decreasing gradient of abundance, both toxins have an apparent high or low “plateau,” with a drastic shift in abundance occurring in central Colorado (approximately the 39th parallel north). Azocasein metalloprotease activity (Additional file 3: Figure S1a; R2 = 0.323, p < 0.001), thrombin-like (Additional file 3: Figure S1b; R2 = 0.181, p < 0.001), and kallikrein-like (Additional file 3: Figure S1c; R2 = 0.13, p < 0.001) serine protease activity, as well as phosphodiesterase activity (Additional file 3: Figure S1d; R2 = 0.052, p < 0.001), were significantly correlated with latitude. Neither PLA2 (Additional file 3: Figure S1e; R2 = 0.001, p = 0.65) nor L-AAO (Additional file 3: Figure S1f; R2 = 0.001, p = 0.68) activities significantly correlated with latitude. Using only toxins that were significant in the both the Mantel test and in the linear regressions, cluster analysis in PC-ORD 7.07 broadly defined two major venom phenotypes that explained 100% of the variation (Table 1). These groups represented a high-myotoxin venom profile and a high-SVMP (and also high SVSP) profile and were used in MAXENT analysis, which requires inputs to be categorical (by species typically, but in this case venom phenotype).
Phylogeographic and population genetic analyses
The phylogeny generated from our nuclear SNP dataset shows strong posterior support (> 0.95 PP) for a split between C. v. nuntius and C. v. viridis, with the C. v. viridis clade containing three subclades of varying support. One of these subclades corresponds to samples from more northern localities, including individuals from Montana, Nebraska, and northern Colorado; this northern subclade is supported with low posterior support (> 0.50 PP; Fig. 4a). Also, within the C. v. viridis clade, populations from southern Colorado, New Mexico, Texas, and Oklahoma form clades that diverge more basally within the larger C. v. viridis clade (Fig. 4a). Our population cluster results from STRUCTURE indicated that a K value of 2 was the best supported model based on the Evanno method (ΔK) across all K values tested (Fig. 4b), although we show results for additional values of K further to investigate evidence for potential structure across these populations. Results from the K = 2 model shows substantial differentiation of ancestry coefficients between C. v. viridis and C. v. nuntius, with moderate levels of admixture/differentiation between C. v. nuntius and southern C. v. viridis localities. Our results also show sub-structuring between northern and southern C. v. viridis samples for various models of K (K = 2–4), although no values of K tested clearly demarcate southern and northern populations of C. v. viridis. However, these analyses do highlight gradual shifts in population differentiation along a North–South axis of C. v. viridis populations. Thus, while STRUCTURE analyses provide some evidence for the cohesiveness of northern populations, they do not support a clearly defined break in population differentiation between northern and southern C. v. viridis.
In total, 638 C. v. viridis specimens were examined from seven natural history collections. Stomach contents were recovered from 130 snake stomachs (20.4% of snakes) and comprised four broad taxonomic classes: Aves, Mammalia, Amphibia, and Reptilia. Of these prey classes, Mammalia was the best represented, comprising 78.2% of adult rattlesnake stomach contents. Lizards represented only 10.3% of all prey items recovered from adults; avian prey items were also recovered from five food boli. A single amphibian was found in an adult snake from Hidalgo Co., New Mexico. Other contents found (3.4% of total) included nematode roundworms as well as unidentifiable material likely not prey related.
Only three prey items (all mammals) were recovered from Montana snakes (Fig. 5a). One mammal each was collected from Oklahoma, South Dakota, and Wyoming. Of the 25 prey items collected from Colorado, one was avian (4%) and 24 (96%) were mammalian. Three birds, one lizard, and 18 rodents were recovered from 22 snakes in Kansas (14%, 4.5%, 82%, respectively). Of the identifiable material from New Mexico snakes, one amphibian (4%), one bird (4%), 8 lizards (32%), and 15 (60%) mammals were recovered from 25 snakes (Fig. 5a).
Based on the northern and southern venom classifications that are separated by the 39th parallel north, we investigated the dietary composition of snakes above and below this latitude. The majority of prey items (95%) from snakes in the northern region were mammals, with only 1 bird recovered (Fig. 5a). The diet of snakes south of the 39th latitude was also dominated by mammals (72.6%); however, we also recovered reptilian (14.5%), amphibian (1.6%), and avian prey (6.5%), along with 4.8% unidentifiable material. Discriminant analysis based on latitude confirmed the dietary clusters identified as north and south of the 39th latitude that corresponded to north and south venom phenotypes (p < 0.001; F = 67.94).
We also found that diet has a significant relationship with venom toxin activities and composition overall (F = 207.9, p < 0.001), and diet showed a statistically significant relationship with the toxins that varied geographically (SVMP activity (F = 333), thrombin-like SVSP activity (F = 17.6) and kallikrein-like SVSP activity (F = 37), and percentage of myotoxin (F = 419) and percentage of SVMP (F = 872); p < 0.001 for all).
Lethality (LD50) assays
SVMP-rich venoms had LD50 values of 2.5 μg/g and 2.0 μg/g towards NSA mice and Anolis lizards, respectively (n = 12; Fig. 5b). Myotoxin-rich venoms were more toxic towards lab mice (1.25 μg/g; n = 12) than towards lizards (LD50 = 5.1 μg/g; n = 12). Purified myotoxin a was quite toxic towards mice (2.0 μg/g; n = 9) and was essentially nontoxic to lizards (> > 50 μg/g; n = 4).
Environmental niche modeling
The ENM area under the curve (AUC) was high for both the high myotoxin phenotype (AUC = 0.881) and the high SVMP phenotype (AUC = 0.942), indicating high predictive power of occurrence likelihood (Fig. 6a, b). The high myotoxin phenotype had a higher probability occurrence in the northern regions of the known range of C. v. viridis but extended towards northern New Mexico, northern Arizona, Utah, and eastern Nevada (Fig. 6a). The high SVMP phenotype had a higher probability of occurrence in the known southern range of this species (Fig. 6b) from Colorado to southern New Mexico and southwestern Texas. The predicted range extended partially into Utah and Nevada and had a high probability up to central Colorado.
The environmental variables responsible for these predicted distributions differed between venom phenotypes (Fig. 6a, b; Additional file 4: Table S4). Mean temperature of driest quarter (BIO 9) contributed the most to the range of the northern phenotype (percent contribution = 19.9%), followed by mean diurnal range (BIO 2; 17.5%) and annual mean temperature (BIO 1; 15.2%). Mean diurnal range contributed the most to the predicted southern phenotype distribution (BIO 2; 42.5%), followed by mean temperature of coolest quarter (BIO 11; 22.8%). Interestingly, the third most significant variable predicting the northern phenotype (BIO 1: annual mean temperature) contributed the least (0%) to predicting the southern phenotype.
Reptile and mammal richness was significantly predicted by environmental variables (p < 0.001; adjusted R2 = 0.947, R2 = 0.859, respectively). Specifically, reptile species richness was significantly correlated with BIO 1 (p < 0.001), BIO 2 (p = 0.001), BIO 4 (p < 0.001), BIO 10 (p < 0.001), BIO 11 (p < 0.001), BIO 15 (p < 0.001), BIO 17 (p < 0.001), BIO 18 (p = 0.044), and BIO 19 (p < 0.001), while mammal species richness was correlated only with BIO 12 (p < 0.001), BIO 17 (p = 0.008), and BIO 18 (p < 0.001). When determining relationships of dietary availability to venom phenotypes, there was a significant contribution to the final multinomial logistic regression model (p < 0.001) of both mammal species richness (p = 0.007) and reptile species richness (p < 0.001).
In order to investigate the influence of diet and environment on venom composition we performed hierarchical linear regressions with two predictor models: (1) only reptile and mammal species richness and (2) Species richness combined with environmental variables. Both models significantly predicted the percentage of myotoxin and SVMP (p < 0.001), but the second model had a higher R2 value in both cases (0.696 > 0.419, and 0.804 > 0.799, respectively). Azocasein SVMP, thrombin-like SVSP, kallikrein-like SVSP activities were also significantly predicted by both models (p < 0.001) with a higher R2 value for model 2 compared to model 1 (0.747 > 0.433; 0.273 > 0.1, 0.236 > 0.181, respectively). PDE activity and percentage and activity of L-AAO was only predicted by the second model (p < 0.001, p = 0.044, and p < 0.001, respectively) and PLA2 did not have a significant relationship with either model.
Unraveling the multifaceted and nuanced abiotic and biotic factors that shape complex phenotypes is a major hurdle to understanding complex trait evolution, and our results demonstrate the influence of two of these factors, diet and environment, in shaping wide-ranging geographic venom variation in the Prairie Rattlesnake. Previous systems have identified a relationship between temperature-related variables and venom composition; however, these patterns did not extend to diet . Therefore, our results provide a critical link between abiotic and biotic influences and their combined effects on venom composition and demonstrate that selection on venom variation is context-dependent even in closely related populations.
We found two distinct venom phenotypes in C. v. viridis: a northern myotoxin-rich phenotype and a southern SVMP-rich phenotype, which are distributed along a latitudinal gradient that corresponds with climatic and prey variation. The diversity of components in pit viper venoms has been shown to correlate with dietary complexity ; however, this does not appear to be a trend in our study. In this case, both venom phenotypes are dominated by a single toxin family (~ 50% myotoxins or SVMPs), distantly followed by other toxins of much lower abundance. Our analyses also identified the toxin responsible (myotoxin a) for the mammal-specific toxicity observed in the northern phenotype, which to our knowledge is the first prey-specific toxin identified from a rattlesnake. These geographic patterns in venom variation in C. viridis add an intriguing new level of complexity to a broader trend observed in rattlesnake venom variation [18, 21, 46,47,48]. The dichotomy observed effectively represents a trade-off between rapid-acting neurotoxic venom (type 2) and enzyme-rich degradative venom (type 1; sensu ). Though the myotoxin-rich northern phenotype lacks the classical PLA2-based neurotoxin characteristic of type 2 venoms, its dominant component is a rapidly-acting toxin that causes near-immediate muscular tetanic paralysis .
The maintenance of type 1 and 2 venom phenotypes has been studied in other rattlesnake species and are thought to be driven by either disruptive selection against type 1 + 2 venoms  or localized directional selection for a particular venom type , resulting in only rare cases of intermediate phenotypes. Regardless of the underlying selection regime, we observe similar evidence in C. v. viridis for the presence of either type 1 or 2 venom across populations, and a rapid transition between these phenotypic extremes within a narrow geographic range of central Colorado, indicating a steep cline in adaptive advantage across this latitudinal threshold. Interestingly, we found some evidence for North–South variation in population genetic structure, although the present-day population genetic differentiation across this narrow region (central Colorado) is not strong or clearly defined. It is possible that past population structure across this North–South axis, which is now potentially obscured by more recent gene flow, may have contributed to the differentiation and divergence of venom phenotypes across this axis. It therefore remains an open question whether past population structure contributed to venom variation divergence, but evidence that population structure is not strong currently across this region suggests that this phenotype transition is likely maintained by relatively strong selection for divergent venom phenotypes on either side of this transition.
In spite of its large geographic range and substantial variation in venom composition, population genetic analyses suggest relatively little present-day population genomic structure across geographically distant populations of C. v. viridis. This absence of well-defined population structure yet high venom variation is similar to patterns observed in the closely related Mojave Rattlesnake (C. s. scutulatus; [19, 20]), which also possesses populations showing type 1 or type 2 venoms. Taken together, these similar examples suggest that selective forces that operate on venom variation are likely very strong, to the extent that distributions of these traits are not readily explained by population genetic structure, and instead that alleles that encode variants underlying venom composition may readily defy gene flow, due primarily to the impacts of local and regional selection regimes on related venom phenotypes.
Our findings highlight that a latitudinal shift in diet corresponds to a switch from SVMP-rich venom to myotoxin-rich venom. All gut contents from the northern range (Montana, South Dakota, northern Colorado) were rodents, with the exception of a single bird. A higher diversity of prey type and a considerably higher abundance of lizards were found in snakes from southern areas (New Mexico, southern Colorado, Kansas, Oklahoma). Prey availability and diversity play a critical role in the selective processes that shape venom composition in some systems [2, 31, 32]; therefore, it is important to note that these dietary patterns likely do not reflect a difference in prey preference but a difference in prey availability, as the proportion of available reptile prey species (particularly lizards) decreases precipitously moving north as a general rule [50,51,52]. In fact, only one species of lizard is found in the northern extreme of C. v. viridis range in southern Canada , making rodents the dominant food source available at higher latitudes (cf. 53). The relationship between diet and venom composition in C. v. viridis represents a departure from the trends displayed in the closely related C. s. scutulatus, where diet did not correlate with venom variation ; however, C. s. scutulatus does not reach the same northern latitudinal extremes of C. v. viridis.
Crotalus s. scutulatus, a sympatric rattlesnake in parts of the extreme southern range of C. v. viridis, still predominantly preys on mammals ; however, lizards have been shown to represent a much larger proportion of diet (up to 44% ). Another potentially sympatric snake in southern regions, C. lepidus klauberi, is found in more isolated and higher elevation regions and predominantly feed on lizards (55.4%), with a much lower proportion of diet consisting of mammals (13.8%; ). In a larger rattlesnake, C. atrox, sympatric in the southern range, mammals comprise the vast majority of diet (94.8%); however, this species is large enough regularly to consume adult lagomorph prey , a prey niche likely not available to the smaller C. v. viridis. Therefore, what prey species are available for consumption by sympatric rattlesnakes and the level of resource niche partitioning is determined by a number of factors, including environmentally constrained species distributions, habitat specialization, and snake size and gape limitations.
A closely related species, C. o. oreganus, which does not overlap C. v. viridis except potentially in narrow zones of Idaho, still shows a comparable wide latitudinal range, reaching from southern California to southern Canada . Though this species does not extend as far south as C. v. viridis, a similar latitudinal shift in diet composition is observed, where mammals constitute virtually the entire diet at northern extremes and lizards represent a more significant (but still not dominant) proportion of prey consumed in the south. A study at the northern extreme (British Columbia) showed that diet consisted only of mammals (95.6%) and birds (4.4%; , and in British Columbia, not a single lizard was recovered from gut contents of roadkill (). In Idaho, diet consisted of 98.1% mammalian, 0.9% avian, and 0.9% reptilian (represented by a single lizard; ) prey. Farther south in California, the diet still consisted primarily of mammals (76.1%); however, lizards comprised a larger proportion of the diet (14.8%), though juveniles consumed lizards more often than adults .
Venomous snakes rely on a chemical rather than a physical means of prey capture; thus, venoms that are effective at rapidly incapacitating available prey might be assumed to be highly advantageous. Indeed, broad relationships between higher venom toxicity and prey specialization in snakes have been observed in fish-eaters , bird-eaters [25, 37], lizard-eaters , and arthropod specialists [32, 62]. The evolution of these taxon-specific venoms can even lead to cases where venom is extremely toxic towards one or a few prey items and is less toxic towards others [25,26,27, 35, 37] due to the production of prey-specific toxins. Crotalus v. viridis follows this trend but in a unique way, by producing a venom dominated by myotoxin a, perhaps the only demonstrated mammal-specific toxin in a rattlesnake, in northern latitudes where mammals represent the vast majority of available prey. Purified myotoxin a is potently toxic towards mice and is known to induce rapid tetanic paralysis ; however, it is virtually nontoxic towards lizards, and the myotoxin-rich venoms are considerably less toxic toward lizards. Conversely, in southern areas, where nonmammals make up a larger proportion of the diet, C. v. viridis produces a venom phenotype effective against both lizards and mammals. This suggests (1) dietary specialization on mammals is a driving factor determining the venom phenotype of C. v. viridis at higher latitudes and (2) a more diverse diet with a higher proportion of non-mammalian prey in lower latitudes is related to expression of SVMP-rich venoms.
In this study, we utilize lethal toxicity as a proxy for diet-specific venom specialization. However, time-to-incapacitation as opposed to 24-h survivability may be a better overall indicator of venom functionality. While we did not quantify time-to-incapacitation during our LD50 experiments, we noticed rapid incapacitation due to apparent tetany of mice injected with myotoxin-rich venom or purified myotoxin a (as previously observed ), but we did not notice a similar effect in lizards even at high doses of purified myotoxin a. This suggests that the myotoxin-rich venom is both more immobilizing and more toxic to mice than SVMP-rich venoms. However, other factors may compensate for a minimally effective venom towards less susceptible (or more dangerous) prey types. Snakes can modify their predatory repertoire by altering post-contact release behavior, strike frequency, and distance depending on the prey type [35, 63] or by injecting a larger volume of venom . Therefore, snakes in northern regions may behaviorally compensate for the lower effectiveness of their venom towards lizards, retaining their potential as a feasible prey type when encountered.
A complete view of community ecology incorporates both biotic and abiotic factors, but the abiotic processes controlling species interactions have historically been underappreciated . High myotoxin and high SVMP venoms are predicted in distinct ecological niches that mirror the latitudinal shifts in phenotype we observed, indicating a relationship between environment and venom composition. We find a variety of temperature variables contributed to our ENM models (though the most significant variables differed between venom phenotype models) and that temperature fluctuations appear particularly important . The overall trend suggesting that temperature is relevant to the distribution of venom phenotypes across a broad range may be explained by two non-exclusive factors: (1) the temperature-dependence of venom toxin functionality and/or (2) the importance of abiotic factors determining aspects of both predator and prey ecology and distribution .
SVMP-rich type 1 venoms have been suggested to play a part in prey predigestion in cooler environmental conditions [18, 66, 67]; however, their role as a digestive agent remains unclear [68, 69]. Our findings do not support a predigestive role of SVMPs in cooler environments, as SVMP-dominated degradative venoms in C. v. viridis are more common in warmer lower latitudes. Enzyme-dependent predigestion is optimal at higher temperatures  and SVMP-rich venoms are likely more efficient at southern latitudes in general. However, the repeated recovery of variables associated with temperature fluctuations (e.g., mean diurnal range , mean temperature of coldest quarter) in enzyme-rich Type 1 venom niche models could indicate the adaptive utility of an enzyme-rich venom in areas of greater temperature variation.
Rattlesnake predatory behavior includes strike-induced chemosensory searching [71,72,73,74], a potentially lengthy process where prey must be relocated after the classic viperid strike-and-release. This may involve following a chemical trail for hours without success, incurring the metabolic costs of activity and later venom production, which may be high regardless of external conditions [75, 76]. Furthermore, it is possible that in cooler northern climates, snakes have a slower response time during predatory events , which may also be problematic for prey relocation. In C. cerastes, less than half of bites result in successful subjugation and feeding in part due to the antipredator maneuvers displayed by some rodent prey that cause premature fang withdrawal and incomplete envenomation [78, 79]. Venom that rapidly incapacitates mammals (virtually the only food source at northern extremes) even at low doses via the mammal-specific myotoxin a would provide an adaptive edge over SVMP-rich venoms since they require a lower volume of venom and less chemosensory searching, resulting in more efficient feeding for low-metabolism snakes in cooler climates.
A more likely scenario is that temperature-related abiotic variables that shift with latitude have a cascading influence on prey distribution and therefore have an indirect impact on numerous toxin families that define venom phenotypes. Latitude is a major defining factor in the distribution of nearly all species ; however, global ectotherm species richness is more strongly related to latitudinal shifts in temperature [52, 53]. Due to the steep decline in ectotherm diversity as climactic temperatures decrease, snakes become increasingly reliant on mammalian prey at higher latitudes. Our data suggests that this pressure drove the evolution of mammal-specific toxicity in northern venoms and generalist toxicity in southern venoms. This points to dietary availability, as determined by latitudinal temperature gradients, as a key driving factor shaping the geographic variation in venom composition. The fact that model performance predicting venom composition increased with the addition of environmental variables to species richness suggests the importance of abiotic factors combined with dietary availability in shaping venom composition. The majority of environmental variables that predicted venom phenotype also predicted species richness, suggesting that venom phenotype is largely indirectly impacted by environment and is due to the more direct influence of abiotic variables on prey availability.
Decades of research has contributed to our contemporary understanding that snake venom composition is remarkably variable between and even within species [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22], yet the biotic and abiotic factors that drive this extreme variation have remained largely speculative [13, 18,19,20, 80]. Our findings provide another key example detailing the extensive venom variation within a single species. We find evidence that venom variation correlates both with abiotic climatic factors and prey type distributions and that variation in venom composition may be linked to differential toxicity for distinct prey taxa (i.e., reptiles versus mammals). These links between venom variation and variation in biotic and abiotic factors further suggest that venom variation likely results from substantial geographic variation in selection regimes that dictate the efficacy of venom phenotypes at the toxin family level across the range of snake species. This highlights the implicit influence of abiotic factors on venom phenotype and provides further evidence for a dominant role of local selection as a key driver of venom variation.
Snake collection and venom extraction
Snakes and venom samples were field collected over a five-year period across the central United States from nine states that span the wide latitudinal range of C. v. viridis (Fig. 1). Snakes were either field-extracted or processed in the UNC Animal Resource Facility and returned to their location of capture within 3 days, or they were maintained in the UNC Animal Resource Facility per UNC IACUC approval (protocol #1302D-SM-S-16). All snake collections and field-collected samples were conducted according to permits obtained from Colorado Parks and Wildlife (scientific collecting permits #14HP974, 17HP0974, 20HP0974, to SPM), Wyoming Game and Fish Department (scientific collecting permit #1544, SPM), New Mexico Game and Fish Department (permit #3418, SPM), Utah Division of Wildlife Resources (permit #1COLL10567, SPM), Texas Parks and Wildlife Department (permit #SPR-0604–391, SPM), and Montana Fish, Wildlife and Parks (permits 2019–075-W and 2020–076-W, TAC). Venom was manually extracted as described ; samples were either dried in the field over calcium carbonate and placed on ice until storage at − 20 °C or were placed in liquid nitrogen until lyophilization and storage at − 20 °C. Due to the potential for ontogenetic venom variation, only venoms from adult snakes were utilized in this study.
Protein concentration determination
Dried venom samples were reconstituted at an approximate concentration of 4.0 mg/mL in Millipore-filtered water, vortexed, and centrifuged for 5 min at 9500 × g. Samples were diluted to an approximate concentration of 0.5 mg/mL and centrifuged again for 5 min at 9500 × g. Protein concentration was then determined in triplicate with a Nanodrop using the Absorbance 280 protocol, and 260/280 ratio was recorded. The amount of material used in all subsequent assays was based on these determinations. Reconstituted stock samples at 4 mg/mL were frozen at − 20 °C until use and then thawed and centrifuged at 9500 × g for 5 min to pellet cellular debris.
Azocasein metalloprotease assay
The azocasein metalloprotease assay protocol was performed as outlined in Aird and da Silva . Briefly, 20 μg of venom was incubated with 1.0 mg of azocasein substrate in 0.5 mL buffer (50 mM HEPES, 100 mM NaCl, pH 8.0) for 30 min at 37 °C. The reaction was stopped with 250 μl of 0.5 M trichloroacetic acid, vortexed, brought to room temperature, and centrifuged at 2000 rpm (~ 700 × g) for 10 min. Absorbance of the supernatant was read at 342 nm, and all values were expressed as ΔA342nm/min/mg venom protein.
L-amino acid oxidase assay
L-amino acid oxidase assays were performed as outlined in Weissbach . L-kynurenine substrate was solubilized at 1.04 mg/mL buffer (50 mM HEPES, 100 mM NaCl, pH 8.0), and 75 μl was added to 20 μg of venom in 645 μl buffer. The reaction was incubated at 37 °C for 30 min and then terminated with 750 μl of 10% trichloroacetic acid. Absorbances were read at 331 nm and specific activity was calculated as nanomoles product formed/min/mg venom protein from a standard curve of the reaction product, kynurenic acid, using the same buffer and acid termination conditions as indicated above.
Phosphodiesterase assays were performed based on Laskowski’s  modification of Björk . Twenty micrograms of crude venom (5 μL) was added to 220 μl of buffer (100 mM tris–HCl, 10 mM MgCl2, pH 9.0). One hundred fifty μl of 1.0 mM bis-p-nitrophenylphosphate substrate was added, and the reaction was incubated at 37 °C for 30 min. The reaction was terminated with 375 μl of 100 mM NaOH containing 20 mM disodium-EDTA, vortexed, brought to room temperature and absorbance read at 400 nm. Specific activity was expressed as ΔA400nm/min/mg venom protein.
Thrombin-like and kallikrein-like serine protease assays
Serine protease assays were performed as described by Mackessy . Eight μg of crude venom (2 μL) was added to 373 μl of buffer (50 mM HEPES, 100 mM NaCl, pH 8.0). Tubes were incubated at 37 °C for approximately 3 min. Fifty μl of thrombin-like substrate (1.0 mM BzPheValArg-pNA; Sigma) or kallikrein-like substrate (1.0 mM Bz-ProPheArg-pNA; Bachem) was added and the sample was vortexed and placed back at 37 °C. Reactions were stopped after 5.0 min with 50% acetic acid. The absorbances for each of the samples were read at 405 nm, and specific activity was calculated from a standard curve of p-nitroaniline and expressed as nanomoles product produced/min/mg of venom.
Phospholipase A2 assay
Phospholipase A2 assays were performed as outlined in Holzer and Mackessy . Briefly, 12.5 μL of crude venom was combined with 583 μL of assay buffer (10 mM tris–HCl, pH 8.0 with 10 mM CaCl2, 100 mM NaCl) followed by 50 μL of 4-nitro-3-(octanoyloxy)benzoic acid substrate. Samples were incubated at 37 °C in duplicate for 20 min and the reaction was stopped with 50 μL of 2.5% Triton X-100. Samples were read at 425 nm and specific activity was recorded as nmol product/min/mg venom protein.
High-performance liquid chromatography
One milligram of venom was subjected to reverse phase HPLC using a Waters system, Empower software, and a Phenomenex Jupiter C18 (250 × 4.6 mm, 5 μm, 300 Å pore size) column as outlined in Smith and Mackessy . Protein/peptide was detected at 280 nm and 220 nm with a Waters 2487 Dual λ Absorbance Detector. Fractions corresponding to each peak were then frozen at − 80 °C overnight, lyophilized and then analyzed along with 20 μg crude venom via SDS-PAGE as described  in order to determine peak complexity, mass and toxin families . Percent peak area at 280 nm was recorded as a proxy for relative toxin abundance.
Myotoxin a purification
Based on protein size (SDS-PAGE), peak 5 was identified as myotoxin-containing and was further purified using reverse phase high-performance liquid chromatography using a Phenomenex Jupiter C18 (250 × 4.6 mm, 5 mm, 300 Å pore size) column, and 1-min fractions were collected at a flow rate of 1 mL/min for a total of 61 min. To elute proteins, a gradient of 95% solution A (0.1% trifluoroacetic acid in Millipore-filtered water) to 100% solution B (0.1% TFA in 100% acetonitrile) was used (5% solution B for 6 min; 5–22% B over 2 min; 22–35% B over 37 min; 35–100% over 2 min; and returning to 5% B over 2 min). Protein/peptide was detected at 220 nm and 280 nm with a Waters 2487 Dual l Absorbance Detector. Fractions corresponding to each peak were pooled, then frozen overnight at − 80 °C and lyophilized.
The dominant peak at 21 min and a secondary peak at 28 min were analyzed via Bruker Microflex LRF MALDI-TOF Proteomics (Proteomics and Metabolomics Facility at Colorado State University (Fort Collins, CO)) to confirm isoform mass identity previously determined by Griffin and Aird . Approximately 1.0 mg protein was dissolved in 1.0 mL 50% acetonitrile containing 0.1% TFA, mixed with 1.0 mL sinapinic acid matrix (10 mg/mL, dissolved in the same solvent), and spotted onto target plates.
A Mantel test was run using Euclidean distance measures in PC-ORD 7.07  for all toxins (using specific enzyme activity or toxin abundance based on HPLC) in order to assess the association between venom components and geographic location (p < 0.05). The PC-ORD enzyme specific activities were normalized to the maximum value for each toxic activity. Linear regression was performed in GraphPad Prism 9® between all enzyme activities and toxin abundances with latitude to determine if toxins varied significantly over latitude. Linear regression t-tests were run to determine if the regression line slope differed significantly from 0. Cluster analysis was run in PC-ORD 7.07  with a Euclidean Distance Measure and Ward’s Group Linkage Method to bin samples into defined venom phenotypes based on the venom toxin abundances and activities that were significant (p < 0.05) both the Mantel test and the linear regressions.
We investigated the diet composition in the geographic ranges that corresponded to the identified venom phenotype groups and performed discriminant analysis in SPSS® to test group membership using latitude as a predictor and bootstrapping with a 95% confidence interval. In order to test the relationship between diet groups and venom variables that significantly differed geographically, we performed a one-way ANOVA in SPSS® with a significance level of 0.05.
Eight representative venom samples from New Mexico, Colorado, and Montana were subjected to shotgun proteomics. Venoms (20 μg) were resuspended in 8 M urea/0.1 M Tris (pH 8.5) and reduced with 5 mM TCEP (tris (2-carboxyethyl) phosphine) for 20 min at room temperature. Samples were alkylated with 50 mM 2-chloroacetamide for 15 min in the dark at room temperature and then diluted fourfold with 100 mM Tris–HCl (pH 8.5) and trypsin digested at an enzyme/substrate mass ratio of 1:20 overnight at 37 °C. To stop the reaction, samples were acidified with formic acid (FA), and digested peptides were purified with Pierce™ C18 Spin Tips (Thermo Scientific) according to the manufacturer’s protocol. Samples were dried in a speed vacuum and resuspended in 0.1% FA. Liquid chromatography tandem mass spectrometry (LC–MS/MS) was performed using an Easy nLC 1000 instrument coupled with a Q Exactive™ HF Mass Spectrometer (both from ThermoFisher Scientific). For each venom, 3 μg of digested peptides were loaded on a C18 column (100 μM inner diameter × 20 cm) packed in-house with 2.7 μm Cortecs C18 resin and separated at a flow rate of 0.4 μl/min with solution A (0.1% FA) and solution B (0.1% FA in ACN) under the following conditions: isocratic at 4% B for 3 min, followed by 4–32% B for 102 min, 32–55% B for 5 min, 55–95% B for 1 min and isocratic at 95% B for 9 min. Mass spectrometry was performed in data-dependent acquisition (DDA) mode. Full MS scans were obtained from m/z 300 to 1800 at a resolution of 60,000, an automatic gain control (AGC) target of 1 × 106, and a maximum injection time (IT) of 50 ms. The top 15 most abundant precursors with an intensity threshold of 9.1 × 103 were selected for MS/MS acquisition at a 15,000 resolution, 1 × 105 AGC, and a maximal IT of 110 ms. The isolation window was set to 2.0 m/z and ions were fragmented at a normalized collision energy of 30. Dynamic exclusion was set to 20 s.
Fragmentation spectra were searched against an in-house Crotalus database containing translated sequences derived from multiple venom transcriptomes (derived from public databases) using the MSFragger-based FragPipe computational platform [91, 92]. Reverse decoys and contaminants were included in the search database. Cysteine carbamidomethylation was selected as a fixed modification, oxidation of methionine was selected as a variable modification, and precursor-ion mass tolerance and fragment-ion mass tolerance were set at 20 ppm and 0.4 Da, respectively. Up to 2 missed tryptic cleavages were allowed and the protein-level false discovery rate (FDR) was set to < 1%. The relative abundance of major toxin families was compared across samples using log-transformed total spectral intensity . Raw mass spectrometry proteomics data has been deposited to the Mass Spectrometry Interactive Virtual Environment (MassIVE) repository with the dataset identifier MSV000089601.
Phylogeographic and population genetic analyses
We incorporated previously published and analyzed ddRADseq data  to investigate how phylogeographic and population genetic structure correspond with patterns of venom variation (Additional file 5: Table S5; see 94 for details). We mapped individual samples to the C. viridis reference genome  using the program BWA v.0.7.10  specifying the “mem” option with default settings. We then used a combination of SAMtools and BCFtools to generate pileup alignments of all samples and called variant sites . We filtered out sites that had read depths lower than 5 and retained only biallelic SNPs that were present across at least 90% of samples using VCFtools. We used the program BEAST V.2.5.3  to infer phylogenetic relationships from our SNP dataset by using BEAUti2 to format input files using a JC69 substitution model and a Gamma Site model. We performed 10,000,000 MCMC iterations with a burn-in of 100,000 generations under a Yule model prior. We included data from several samples of C. v. nuntius, historically considered a subspecies of C. viridis, and we designated C. oreganus concolor (a closely related sister species) as the outgroup and visualized phylogenetic trees using TreeAnnotator to produce a maximum clade credibility tree with posterior probabilities (PP) to assess node support. To evaluate population genetic structure, we used the program STRUCTURE  to infer probabilities of individuals of being assigned to K genetic clusters. STRUCTURE uses a Bayesian clustering approach on multilocus genotype data both to identify populations and to assign individuals to a population . For multiple values of K (two through four), we ran 20 iterations using the program StrAuto  with 10,000 burn-in generations followed by 75,000 logged generations. We first used the Evanno method which examines the overall change in likelihood scores across K values to determine the most supported value (ΔK; ) and then combined iterations from the highest supported K value using the “greedy” method in the program CLUMPP  for each K value.
Natural history collections in multiple locations containing fluid-preserved Prairie Rattlesnake specimens were visited to examine specimens (Additional file 6: Table S6). Specimens were considered usable if accompanying data did not indicate that animals were captive or were held in captivity prior to preservation. Before sampling for prey remains, snakes were sexed via presence or absence of hemipenes, and snout-vent length (SVL) and tail length (TL) were recorded by measuring snakes with a soft metric tape measure.
Prey remains were sampled via inspection of stomach contents by making a small (2–6 cm) midventral incision through ventral scales of the snake. Once located, an additional incision was made through the stomach wall to determine if a food bolus was present. If present, direction of ingestion of the prey item was recorded (inferred from orientation in the stomach), and the bolus was removed for identification to greatest taxonomic resolution possible. For the purposes of this study, prey remains were identified based on the presence of fur, scales, or feathers into broad categories (amphibian, reptile, mammal, bird) and were stored in 70% ethanol.
Lethal toxicity (LD50) assays
Venoms from snakes found in Lincoln County, Colorado (to represent the southern venom phenotype), Weld County, Colorado (to represent the northern venom phenotype), and purified myotoxin a were tested against NSA mice and Anolis sp. lizards. Doses from 0 to 10 μg/g (3 mice/dose) were given with an intraperitoneal injection (100 μl bolus of 0.9% saline for mice, 50 μl for lizards) on the right side of the body. NSA mice were used because we have bred them in our facility for many years and they are traceable to Jackson Labs, they originate from an outcrossed line, and they have been a standard strain for much of our research over the last 20 + years; hence, toxicity data of venoms from different species are comparable (same model strain). Mouse LD50 determinations used the Trimmed Spearman–Karber (TSK) program version 1.5 (U.S. Environmental Protection Agency, 1990). Toxicity data was analyzed using “Quest Graph™ LD50 Calculator,” AAT Bioquest, Inc., 1 Feb. 2023. [https://www.aatbio.com/tools/ld50-calculator] For lizards, additional myotoxin a doses of 25 and 50 μg/g were administered in a 50 μL bolus.
Environmental niche modeling
To determine if differences in venom composition correlated with environmental variables, ecological niche models (ENMs) were constructed in MAXENT v3.4.1  for the occurrence of venom phenotypes (as determined by cluster analysis of toxins with significant geographic variation). MAXENT uses occurrence points and a set of input environmental variables to generate a probability of presence distribution map . Locality points were filtered in the R package spThin to a minimum nearest neighbor distance of 10 km in order to reduce sampling bias and improve model performance . Environmental variables were taken from the WorldClim dataset v 1.4  with a 1-km resolution. WorldClim data was processed in ArcGIS Pro to the region of interest, and all layers were modified to the same cell size (1 km), extent, and projected coordinate system. Highly correlated environmental variables (BIO 3, BIO 5, BIO 6, BIO 7, BIO 13, BIO 14, BIO 16), including those that were calculated based on other variables, were removed to avoid biasing model fitting [19, 106]. MAXENT was run with 15 replicates and 5000 iterations. Cross-validation with random seeding was used to determine the best model and area under the curve (AUC) was used to evaluate model performance.
In order to investigate further the relationship between environment and diet with venom composition, we extracted WorldClim dataset v1.4 environmental variables and species richness values from raster layers that corresponded to exact collection localities of our 147 adult samples used for venom analysis [105, 107] and linked these to each venom sampling locality in ArcGIS Pro 3.0. We tested the relationship between diet and environment on all venom components individually with hierarchical linear regressions between each venom component and predictor variables from two models: (1) mammal and reptile species richness alone  and (2) mammal and reptile species richness combined with the WorldClim dataset v 1.4 described above . We also performed linear regression to test the relationship between WorldClim variables and mammal and reptile species richness and further investigated the relationship between dietary availability and venom phenotype with a multinomial logistic regression with a chi-square goodness-of-fit test and a confidence interval of 95% between mammal and reptile species richness  with both venom phenotypes.
Availability of data and materials
All data generated or analyzed during this study are included in this published article, its supplementary information files and publicly available repositories. Proteomics results are available in Additional file 2: Tables S2-S3. Raw mass spectrometry proteomics data has been deposited to the Mass Spectrometry Interactive Virtual Environment (MassIVE) repository with the dataset identifier MSV000089601. The address is as follows: https://massive.ucsd.edu/ProteoSAFe/datasets.jsp#%7B%22query%22%3A%7B%7D%2C%22table_sort_history%22%3A%22createdMillis_dsc%22%7D and the ID number is MSV000089601. Information is also available through Proteome Xchange with the dataset identifier PXD034274 (http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD034274). Final alignment files for genetic analyses can be found on Dryad (https://doi.org/10.5061/dryad.3j9kd51ph).
Snake venom metalloprotease
Snake venom serine protease
- PLA2 :
L-amino acid oxidase
Environmental niche modeling
High-performance liquid chromatography
Area under the curve
Cysteine-rich secretory protein
Meyer JR, Kassen R. The effects of competition and predation on diversification in a model adaptive radiation. Nature. 2007;446(7134):432–5.
Holding ML, Strickland JL, Rautsaw RM, Hofmann EP, Mason AJ, Hogan MP, et al. Phylogenetically diverse diets favor more complex venoms in North American pitvipers. Proc Natl Acad Sci U S A. 2021;118(17):e2015579118.
Schmitz O. Predator and prey functional traits: understanding the adaptive machinery driving predator-prey interactions. F1000Research. 2017;6:1767.
Zancolli G, Casewell NR. Venom systems as models for studying the origin and regulation of evolutionary novelties. Mol Biol Evol. 2020;37(10):2777–90.
Vonk FJ, Casewell NR, Henkel CV, Heimberg AM, Jansen HJ, McCleary RJR, et al. The king cobra genome reveals dynamic gene evolution and adaptation in the snake venom system. Proc Natl Acad Sci U S A. 2013;110(51):20651–6.
Tasoulis T, Isbister GK. A review and database of snake venom proteomes. Toxins. 2017;9(9):290.
Chippaux JP, Williams V, White J. Snake venom variability: methods of study, results and interpretation. Toxicon. 1991;29(11):1279–303.
Tasoulis T, Silva A, Veerati PC, Baker M, Hodgson WC, Dunstan N, et al. Intra-specific venom variation in the australian coastal taipan Oxyuranus scutellatus. Toxins (Basel). 2020;12(8):485.
Alape-Girón A, Sanz L, Escolano J, Flores-Díaz M, Madrigal M, Sasa M, et al. Snake venomics of the lancehead pitviper Bothrops asper. Geographic, individual, and ontogenetic variations. J Proteome Res. 2008;7(8):3556–71.
Daltry JC, Ponnudurai G, Shin CK, Tan NH, Thorpe RS, Wüster W. Electrophoretic profiles and biological activities: Intraspecific variation in the venom of the Malayan pit viper (Calloselasma rhodostoma). Toxicon. 1996;34(1):67–79.
Kalita B, Mackessy SP, Mukherjee AK. Proteomic analysis reveals geographic variation in venom composition of Russell’s Viper in the Indian subcontinent: implications for clinical manifestations post-envenomation and antivenom treatment. Expert Rev Proteomics. 2018;15(10):837–49.
Margres MJ, Walls R, Suntravat M, Lucena S, Sánchez EE, Rokyta DR. Functional characterizations of venom phenotypes in the eastern diamondback rattlesnake (Crotalus adamanteus) and evidence for expression-driven divergence in toxic activities among populations. Toxicon. 2016;119:28.
Sousa LF, Holding ML, Del-Rei THM, Rocha MMT, Mourão RHV, Chalkidis HM, et al. Individual variability in Bothrops atrox snakes collected from different habitats in the Brazilian Amazon: new findings on venom composition and functionality. Toxins (Basel). 2021;13(11):814.
Amazonas DR, Portes-Junior JA, Nishiyama-Jr MY, Nicolau CA, Chalkidis HM, Mourão RHV, et al. Molecular mechanisms underlying intraspecific variation in snake venom. J Proteomics. 2018;181:60–72.
Rashmi U, Khochare S, Attarde S, Laxme RRS, Suranse V, Martin G, et al. Remarkable intrapopulation venom variability in the monocellate cobra (Naja kaouthia) unveils neglected aspects of India’s snakebite problem. J Proteomics. 2021;242(4):383–97.
Gibbs LH, Chiucchi JE. Deconstructing a complex molecular phenotype: population-level variation in individual venom proteins in eastern massasauga rattlesnakes (Sistrurus c. catenatus). J Mol Evol. 2011;72:383–97.
Smiley-Walters SA, Farrell TM, Lisle Gibbs H. High levels of functional divergence in toxicity towards prey among the venoms of individual pigmy rattlesnakes. Biol Lett. 2019;15(2):20180876.
Mackessy SP. Evolutionary trends in venom composition in the Western Rattlesnakes (Crotalus viridis sensu lato): Toxicity vs. tenderizers. Toxicon. 2010;55(8):1463–74.
Strickland JL, Smith CF, Mason AJ, Schield DR, Borja M, Castañeda-Gaytán G, et al. Evidence for divergent patterns of local selection driving venom variation in Mojave Rattlesnakes (Crotalus scutulatus). Sci Rep. 2018;8(1):17622.
Zancolli G, Calvete JJ, Cardwell MD, Greene HW, Hayes WK, Hegarty MJ, et al. When one phenotype is not enough: Divergent evolutionary trajectories govern venom variation in a widespread rattlesnake species. Proc R Soc B Biol Sci. 2019;286(1898):20182735.
Sunagar K, Undheim EAB, Scheib H, Gren ECK, Cochran C, Person CE, et al. Intraspecific venom variation in the medically significant Southern Pacific Rattlesnake (Crotalus oreganus helleri): biodiscovery, clinical and evolutionary implications. J Proteomics. 2014;99:68–83.
Smiley-Walters SA, Farrell TM, Gibbs HL. Evaluating local adaptation of a complex phenotype: reciprocal tests of pigmy rattlesnake venoms on treefrog prey. Oecologia. 2017;184(4):739–48.
Holding ML, Biardi JE, Gibbs HL. Coevolution of venom function and venom resistance in a rattlesnake predator and its squirrel prey. Proc R Soc B Biol Sci. 2016;283(1829):20152841.
Gibbs HL, Mackessy SP. Functional basis of a molecular adaptation: prey-specific toxic effects of venom from Sistrurus rattlesnakes. Toxicon. 2009;53(6):672–9.
Pawlak J, Mackessy SP, Sixberry NM, Stura EA, Le Du MH, Ménez R, et al. Irditoxin, a novel covalently linked heterodimeric three-finger toxin with high taxon-specific neurotoxicity. FASEB J. 2009;23(2):534–45.
Heyborne WH, Mackessy SP. Identification and characterization of a taxon-specific three-finger toxin from the venom of the Green Vinesnake (Oxybelis fulgidus; Family Colubridae). Biochimie. 2013;95(10):1923–32.
Modahl CM, Mrinalini, Frietze S, Mackessy SP. Adaptive evolution of distinct prey-specific toxin genes in rear-fanged snake venom. Proc R Soc B Biol Sci. 2018;285(1884):20181003.
Davies EL, Arbuckle K. Coevolution of snake venom toxic activities and diet: evidence that ecological generalism favours toxicological diversity. Toxins (Basel). 2019;11(12):711.
Robinson KE, Holding ML, Whitford MD, Saviola AJ, Yates JR, Clark RW. Phenotypic and functional variation in venom and venom resistance of two sympatric rattlesnakes and their prey. J Evol Biol. 2021;34(9):1447–65.
Jorge Da Silva N, Aird SD. Prey specificity, comparative lethality and compositional differences of coral snake venoms. Comp Biochem Physiol - C Toxicol Pharmacol. 2001;128(3):425–56.
Daltry JC, Wüster W, Thorpe RS. Diet and snake venom evolution. Nature. 1996;379(6565):537–42.
Barlow A, Pook CE, Harrison RA, Wüster W. Coevolution of diet and prey-specific venom activity supports the role of selection in snake venom evolution. Proc R Soc B Biol Sci. 2009;276(1666):2443–9.
Andrade DV, Abe AS. Relationship of venom ontogeny and diet in Bothrops. Herpetologica. 1999;55(2):200–4.
Sanz L, Lisle Gibbs H, Mackessy SP, Calvete JJ. Venom proteomes of closely related Sistrurus rattlesnakes with divergent diets. J Proteome Res. 2006;5(9):2098–112.
Mackessy SP, Sixberry NM, Heyborne WH, Fritts T. Venom of the Brown Treesnake, Boiga irregularis: Ontogenetic shifts and taxa-specific toxicity. Toxicon. 2006;47(5):537–48.
Cipriani V, Debono J, Goldenberg J, Jackson TNW, Arbuckle K, Dobson J, et al. Correlation between ontogenetic dietary shifts and venom variation in Australian brown snakes (Pseudonaja). Comp Biochem Physiol Part - C Toxicol Pharmacol. 2017;1(197):53–60.
Pawlak J, Mackessy SP, Fry BG, Bhatia M, Mourier G, Fruchart-Gaillard C, et al. Denmotoxin, a three-finger toxin from the colubrid snake Boiga dendrophila (mangrove catsnake) with bird-specific activity. J Biol Chem. 2006;281(39):29030–41.
Holding ML, Margres MJ, Rokyta DR, Gibbs HL. Local prey community composition and genetic distance predict venom divergence among populations of the northern Pacific rattlesnake (Crotalus oreganus). J Evol Biol. 2018;31(10):1513–28.
Modahl CM, Mukherjee AK, Mackessy SP. An analysis of venom ontogeny and prey-specific toxicity in the Monocled Cobra (Naja kaouthia). Toxicon. 2016;119:8–20.
Campbell JA, Lamar WW, Brodie ED. The venomous reptiles of the Western Hemisphere. Ithaca: Comstock Pub. Associates; 2004.
Ownby CL, Cameron D, Tu AT. Isolation of myotoxic component from rattlesnake (Crotalus viridis viridis) venom. Electron microscopic analysis of muscle damage. Am J Pathol. 1976;85(1):149–66.
Ownby CL, Colberg TR. Classification of myonecrosis induced by snake venoms: venoms from the Prairie rattlesnake (Crotalus viridis viridis), Western diamondback rattlesnake (Crotalus atrox) and the Indian cobra (Naja naja naja). Toxicon. 1988;26(5):459–74.
Saviola AJ, Pla D, Sanz L, Castoe TA, Calvete JJ, Mackessy SP. 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 Proteomics. 2015;121:28–43.
Bober MA, Glenn JL, Straight RC, Ownby CL. Detection of myotoxin a-like proteins in various snake venoms. Toxicon. 1988;26(7):665–73.
Tsai IH, Wang YM, Chen YH. Variations of phospholipases A2 in the geographic venom samples of pitvipers. J Toxicology - Toxin Reviews. 2003;22(4):651–62.
Massey DJ, Calvete JJ, Sánchez EE, Sanz L, Richards K, Curtis R, et al. Venom variability and envenoming severity outcomes of the Crotalus scutulatus scutulatus (Mojave rattlesnake) from Southern Arizona. J Proteomics. 2012;75(9):2576–87.
Glenn JL, Straight RC, Wolt TB. Regional variation in the presence of canebrake toxin in Crotalus horridus venom. Comp Biochem Physiol Part C Pharmacol. 1994;107(3):337–46.
Borja M, Neri-Castro E, Castañeda-Gaytán G, Strickland JL, Parkinson CL, Castañeda-Gaytán J, et al. Biological and proteolytic variation in the venom of Crotalus scutulatus scutulatus from Mexico. Toxins (Basel). 2018;10(1):35.
Bieber AL, Nedelkov D. Structural, biological and biochemical studies of myotoxin a and homologous myotoxins. J Toxicol - Toxin Rev. 1997;16(42006):33–52.
Schall JJ, Pianka ER. Geographical trends in numbers of species. Science. 1978;201(4357):679–86.
Jetz W, Fine PVA. Global gradients in vertebrate diversity predicted by historical area-productivity dynamics and contemporary environment. PLoS Biol. 2012;10(3):e1001292.
Darlington PJ. The geographical distribution of cold-blooded vertebrates. Q Rev Biol. 1948;23(1):1–26.
Bleakney S. A zoogeographical study of the amphibians and reptiles of Eastern Canada. Musée Natl du Canada. 1958;Bull no 155:119 p.
Salazar J, Lieb C. Geographic diet variation of Mojave rattlesnake (Crotalus scutulatus). El Paso: University of Texas at El Paso; 2003.
Holycross AT, Painter CW, Prival DB, Swann DE, Schroff MJ, Edwards T, et al. Diet of Crotalus lepidus klauberi (Banded Rock Rattlesnake). J Herpetol. 2002;36(4):589–97.
Loughran CL, Nowak EM, Schofer J, Sullivan KO, Sullivan BK. Lagomorphs as prey of western diamond-backed rattlesnakes (Crotalus atrox) in Arizona. Southwest Nat. 2013;58(4):502–5.
Macartney JM. Diet of the Northern Pacific rattlesnake, Crotalus viridis oreganus, in British Columbia. Herpetol. 1989;45(3):299–304.
McAllister JM, Maida JR, Dyer O, Larsen KW. Diet of roadkilled Western Rattlesnakes (Crotalus oreganus) and Gophersnakes (Pituophis catenifer) in southern British Columbia. Northwest Nat. 2016;97(3):181–9.
Wallace RL, Diller LV. Feeding ecology of the rattlesnake, Crotalus viridis oreganus, in northern Idaho. J Herpetol. 1990;24(3):246.
Sparks AM, Lind C, Taylor EN. Diet of the Northern Pacific rattlesnake (Crotalus o. oreganus) in California. Herpetol Rev. 2015;46(2):161–5.
Glodek GS, Voris HK. Marine snake diets: prey composition, diversity and overlap. Copeia. 1982;1982(3):661.
Richards DP, Barlow A, Wüster W. Venom lethality and diet: differential responses of natural prey and model organisms to the venom of the saw-scaled vipers (Echis). Toxicon. 2012;59(1):110–6.
Farrell TM, Smiley-Walters SA, McColl DE. Prey species influences foraging behaviors: rattlesnake (Sistrurus miliarius) predation on little brown skinks (Scincella lateralis) and giant centipedes (Scolopendra viridis). J Herpetol. 2018;52(2):156–61.
Hayes W. The snake venom-metering controversy: levels of analysis, assumptions, and evidence. In: W.K. Hayes, K.R. Beaman, M.D. Cardwell, and S.P. Bush (editors), The Biology of Rattlesnakes. Loma Linda University Press, Loma Linda, CA. 2008;1–30.
Dunson WA, Travis J. The role of abiotic factors in community organization. Am Nat. 1991;138(5):1067–91.
Thomas RG, Pough FH. The effect of rattlesnake venom on digestion of prey. Toxicon. 1979;17(3):221–8.
Mackessy SP. Venom composition in rattlesnakes: trends and biological significance. In: W.K. Hayes, K.R. Beaman, M.D. Cardwell, and S.P. Bush (editors), The Biology of Rattlesnakes, Loma Linda University Press, Loma Linda, CA. 2008;495-510.
Chu CW, Tsai TS, Tsai IH, Lin YS, Tu MC. Prey envenomation does not improve digestive performance in Taiwanese pit vipers (Trimeresurus gracilis and T. stejnegeri stejnegeri). Comp Biochem Physiol - A Mol Integr Physiol. 2009;152(4):579–85.
McCue MD. Prey envenomation does not improve digestive performance in western diamondback rattlesnakes (Crotalus atrox). J Exp Zool Part A Ecol Genet Physiol. 2007;307(10):568–77.
Secor SM. Digestive physiology of the Burmese python: broad regulation of integrated performance. J Exp Biol. 2008;211(24):3767–74.
Reinert HK, Cundall D, Bushar LM. Foraging behavior of the Timber rattlesnake, Crotalus horridus. Copeia. 1984;1984(4):976.
Chiszar D, Radcliffe CW, Byers T, Stoops R. Prey capture behavior in nine species of venomous snakes. Psychol Rec. 1986;36(4):433–8.
Chiszar D, Radcliffe CW, Scudder KM. Analysis of the behavioral sequence emitted by rattlesnakes during feeding episodes. I. Striking and chemosensory searching. Behav Biol. 1977;21(3):418–25.
Schwenk K. Of tongues and noses: chemoreception in lizards and snakes. Trends Ecol Evol. 1995;10(1):7–12.
McCue MD. Cost of producing venom in three North American pitviper species. Copeia. 2006;2006(4):818–25.
Pintor AFV, Krockenberger AK, Seymour JE. Costs of venom production in the common death adder (Acanthophis antarcticus). Toxicon. 2010;56(6):1035–42.
Vincent SE, Mori A. Determinants of feeding performance in free-ranging pit-vipers (Viperidae: Ovophis okinavensis): key roles for head size and body temperature. Biol J Linn Soc. 2008;93(1):53–62.
Whitford MD, Freymiller GA, Clark RW. Managing predators: the influence of kangaroo rat antipredator displays on sidewinder rattlesnake hunting behavior. Ethology. 2019;125(7):450–6.
Whitford MD, Freymiller GA, Clark RW. Avoiding the serpent’s tooth: predator–prey interactions between free-ranging sidewinder rattlesnakes and desert kangaroo rats. Anim Behav. 2017;130:73–8.
Mackessy SP, Williams K, Ashton KG. Ontogenetic variation in venom composition and diet of Crotalus oreganus concolor. A case of venom paedomorphosis? Copeia. 2003;4:769–82.
Mackessy SP. Venom ontogeny in the Pacific rattlesnakes Crotalus viridis helleri and C. v. oreganus. Copeia. 1988;1988(1):92.
Aird SD, da Silva NJ. Comparative enzymatic composition of Brazilian coral snake (Micrurus) venoms. Comp Biochem Physiol - Part B Biochem. 1991;99(2):287–94.
Weissbach H, Robertson AV, Witkop B, Udenfriend S. Rapid spectrophotometric assays for snake venom l-amino acid oxidase based on the oxidation of l-kynurenine or 3,4-dehydro-l-proline. Anal Biochem. 1960;1(4–5):286–90.
Laskowski M. Purification and properties of venom phosphodiesterase. Methods Enzymol. 1980;65(C):276–84.
Bjork W. Purification of phosphodiesterase from Bothrops atrox venom, with special consideration of the elimination of monophosphatases. J Biol Chem. 1963;238:2487–90.
Mackessy SP. Kallikrein-like and thrombin-like proteases from the venom of the northern Pacific rattlesnake (Crotalus viridis oreganus). J Nat Toxins. 1993;2(2):223–39.
Holzer M, Mackessy SP. An aqueous endpoint assay of snake venom phospholipase A2. Toxicon. 1996;34(10):1149–55.
Smith CF, Mackessy SP. The effects of hybridization on divergent venom phenotypes: characterization of venom from Crotalus scutulatus scutulatus × Crotalus oreganus helleri hybrids. Toxicon. 2016;120:110–23.
Griffin PR, Aird SD. A new small myotoxin from the venom of the Prairie rattlesnake (Crotalus viridis viridis). FEBS Lett. 1990;274(1–2):43–7.
Grandin U. PC-ORD version 5: a user-friendly toolbox for ecologists. J Veg Sci. 2006;17(6):843–4.
Yu F, Teo GC, Kong AT, Haynes SE, Avtonomov DM, Geiszler DJ, et al. Identification of modified peptides using localization-aware open search. Nat Commun. 2020;11(1):1–9.
Kong AT, Leprevost FV, Avtonomov DM, Mellacheruvu D, Nesvizhskii AI. MSFragger: ultrafast and comprehensive peptide identification in mass spectrometry-based proteomics. Nat Methods. 2017;14(5):513–20.
Zhu W, Smith JW, Huang CM. Mass spectrometry-based label-free quantitative proteomics. J Biomed Biotechnol. 2010;2010: 840518.
Schield DR, Perry BW, Adams RH, Card DC, Jezkova T, Pasquesi GIM, et al. Allopatric divergence and secondary contact with gene flow: a recurring theme in rattlesnake speciation. Biol J Linn Soc. 2019;128(1):149–69.
Schield DR, Card DC, Hales NR, Perry BW, Pasquesi GM, Blackmon H, et al. The origins and evolution of chromosomes, dosage compensation, and mechanisms underlying venom regulation in snakes. Genome Res. 2019;29(4):590–601.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, et al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10(4):e1003537.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.
Chhatre VE, Emerson KJ. StrAuto: Automation and parallelization of STRUCTURE analysis. BMC Bioinformatics. 2017;18(1):1–5.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14(8):2611–20.
Jakobsson M, Rosenberg NA. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23(14):1801–6.
Hijmans RJ, Phillips S, Leathwich, Elith J. dismo: species distribution modeling. R version 1.1–4. Reference Module in Life Sciences. 2017.
Aiello-Lammens ME, Boria RA, Radosavljevic A, Vilela B, Anderson RP. spThin: an R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography (Cop). 2015;38(5):541–5.
Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37(12):4302–15.
Dormann CF, Elith J, Bacher S, Buchmann C, Carl G, Carré G, et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography (Cop). 2013;36(1):27–46.
Jenkins CN, Van Houtan KS, Pimm SL, Sexton JO. US protected lands mismatch biodiversity priorities. Proc Natl Acad Sci U S A. 2015;112(16):5081–6.
The authors thank the Loose Family, Palmer Family, and Gottsch Family for access to sampling locations. We also thank Joe Ehrenberger, Bryon Shipley, and Ryan Borgmann with Adaptation Environmental Services and Rocky Mountain Rattlesnake Services for their help collecting snakes in the Denver area. The assistance of Graham Dawson and Kaleb Hill in the field is gratefully acknowledged, as is the assistance in the lab (SPM) of undergraduates Keira Lopez, Parker Gottsch, Forrest Quinn Tappy, David Ceja Galindo and Quinn Cochran.
This work was supported by a grant from the NSF (DEB Phylogenetic Systematics; DEB-1655571) to TAC and SPM. CFS received further financial support from the UNC Graduate Student Association and the UNC College of Natural and Health Sciences.
Ethics approval and consent to participate
All snake collections and field-collected samples were conducted according to permits obtained from the Colorado Parks and Wildlife (scientific collecting permits #14HP974, 17HP0974, 20HP0974, to SPM), the Wyoming Game and Fish Department (scientific collecting permit #1544, SPM), New Mexico Game and Fish Department (permit #3418, SPM), Utah Division of Wildlife Resources (permit # 1COLL10567, SPM), Texas Parks and Wildlife Department (permit #SPR-0604–391, SPM), Montana Permit (2020–076-076-W and 2019–075-W, TAC). Venom extractions and toxicity tests were performed in accordance with protocols (1905D-SM-SbirdsLM-22, 2004D-SM-S-23) approved by the University of Northern Colorado Institutional Animal Care and Use Committee.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Venom Sample Information. Collection information for all individuals included in study.
Table S2. Proteomics Data. List of venom proteins identified by LC-MS/MS with total signal intensity values for each protein. Table S3. Proteomics Data. List of venom toxin families identified by LC-MS/MS with the summed total signal intensity values for all proteins in a family.
Figure S1. Enzyme Activities. Linear regressions of latitude with major enzyme toxins. a) Azocasein snake venom metalloprotease*, b) Thrombin-like serine protease*, c) Kallikrein-like serine protease*, d) Phosphodiesterase*, e) Phospholipase A2 and f) L-amino acid oxidase. * = p<0.05.
Table S4. Environmental Modeling. Percent contributions of WorldClim environmental layers to ENM predictions.
Table S5. Genomic Information. Species ID, project ID, locality, tree label, SRA accession, and reference study for genomic data.
Table S6. Diet Data Collection Information. Museum collections investigated, collection location, and number of specimens examined.
About this article
Cite this article
Smith, C.F., Nikolakis, Z.L., Ivey, K. et al. Snakes on a plain: biotic and abiotic factors determine venom compositional variation in a wide-ranging generalist rattlesnake. BMC Biol 21, 136 (2023). https://doi.org/10.1186/s12915-023-01626-x
- Adaptive trait
- Phenotypic variation