Next Article in Journal
Calcination Improves the In Vivo Efficacy of a Montmorillonite Clay to Bind Aflatoxin G1 in Broiler Chickens: A Toxicokinetic Approach
Previous Article in Journal
Immunotoxin Screening System: A Rapid and Direct Approach to Obtain Functional Antibodies with Internalization Capacities
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Gradual and Discrete Ontogenetic Shifts in Rattlesnake Venom Composition and Assessment of Hormonal and Ecological Correlates

1
Department of Biological Sciences, Florida State University, Tallahassee, FL 32304, USA
2
School of Natural Resources and Environment, University of Florida, Gainesville, FL 32611, USA
3
Department of Organismic and Evolutionary Biology, Harvard University, Cambridge, MA 02138, USA
4
Department of Integrative Biology, University of South Florida, Tampa, FL 33620, USA
*
Author to whom correspondence should be addressed.
Joint lead authors.
Submission received: 18 August 2020 / Revised: 7 October 2020 / Accepted: 13 October 2020 / Published: 16 October 2020
(This article belongs to the Section Animal Venoms)

Abstract

:
Ontogenetic shifts in venom occur in many snakes but establishing their nature as gradual or discrete processes required additional study. We profiled shifts in venom expression from the neonate to adult sizes of two rattlesnake species, the eastern diamondback and the timber rattlesnake. We used serial sampling and venom chromatographic profiling to test if ontogenetic change occurs gradually or discretely. We found evidence for gradual shifts in overall venom composition in six of eight snakes, which sometimes spanned more than two years. Most chromatographic peaks shift gradually, but one quarter shift in a discrete fashion. Analysis of published diet data showed gradual shifts in overall diet composition across the range of body sizes attained by our eight study animals, while the shifts in abundance of different prey classes varied in form from gradual to discrete. Testosterone concentrations were correlated with the change in venom protein composition, but the relationship is not strong enough to suggest causation. Venom research employing simple juvenile versus adult size thresholds may be failing to account for continuous variation in venom composition lifespan. Our results imply that venom shifts represent adaptive matches to dietary shifts and highlight venom for studies of alternative gene regulatory mechanisms.
Key Contribution: We found that differences in venom composition, including the phospholipase A2 and snake venom metalloproteinase enzymes that contribute significantly to venom function, were correlated with increases in body size. These changes tended to occur gradually rather than discretely, although several clear examples of discrete shifts in individual peaks occurred, suggesting regulatory mechanisms may differ among venom proteins. Both gradual and discrete shifts in venom protein expression may be adaptive for the similar modes of shift in diet composition.

1. Introduction

Ontogenetic shifts are common in nature [1], and can optimize resource use, defense, or other ecological function across life stages when they are adaptive. Optimal trait performance is expected to differ most across life stages in organisms that metamorphose, or in size-structured populations that are defined by altered patterns of resource use throughout development [2]. Ontogenetic shifts can be described by the nature of their trajectory: namely, whether they occur suddenly or gradually [3]. Discrete ontogenetic shifts involve a number of rapidly occurring alterations, such as in the metamorphosis of insects [4] or the shift in the diet of juvenile yellowfin tuna when they reach 1.5 kg [5], and correspond to sudden changes in life-history features [6]. Gradual shifts are known from many plant species, where defense traits such as palatability to herbivores and phytochemical profiles change slowly [7]. Gradual ontogenetic shifts in animal phenotypes that match gradually shifting resource use would provide models for understanding the proximate and ultimate causes of gradual ontogeny in animals’ systems.
Ontogenetic shifts are well-known in the venom composition of animals such as snakes [6,8,9,10,11,12,13,14,15,16,17,18,19] and scorpions [20]. Animal venoms are composed of an assortment of proteins and peptides that are produced in a gland and injected into prey or in defense against predators [21]. Ontogenetic shifts in composition can yield measurable variation in venom function as stark as shifts from primarily neurotoxic and myotoxic to primarily hemotoxic and hemorrhagic venoms in snakes [22,23,24]. Venom ontogeny has medical relevance, as snakebite outcomes and antivenom efficacy can depend upon venom composition [25,26,27,28]. It also has biological relevance, as ontogenetic shifts in venom composition can impact the ability of snakes to incapacitate certain species of prey [22]. Although the presence of shifts has been demonstrated across many species, less is known about their tempo, either as gradual or sudden processes, and the proximate physiological signals that might initiate shifts in venom composition remain largely unexplored [10,29].
Studies on both the evolution and biomedical application of venom often categorize animals as juvenile or adult based on size cutoffs, such as size at reproductive maturity [24,30,31,32,33,34,35]. This logic inherently implies a rapidly shifting mode for the ontogenetic shift in venom composition. If the trajectory of the shift were more continuous, this would confound studies of venom variation by introducing variation from mid-sized animals that blurs the lines between the juvenile and adult study groups. Previous work in Pacific rattlesnakes (Crotalus oreganus ssp.) suggested that protease activity in venoms begins to increase between 50 and 60 cm in snout–vent length (SVL) [22], while in Crotalus adamanteus there is evidence for a shift in venom composition at a length of 102 cm, the size that indicates sexual maturity in the species [30]. Prior works on ontogenetic shifts in venom across developmental stages tend to employ a “venom census”, surveying populations by collecting single samples from many individuals of different sizes. However, the census approach is limited by both significant individual variation in venoms [32,36], the potential presence of a terminal venom phenotype in otherwise indeterminately growing snakes [37], and the potential for size-biased sampling. Therefore, census or population-snapshot surveys cannot account for potential individual variation in tempo and mode of ontogenetic shifts, obscuring our ability to distinguish between gradual and discrete shifts in venom phenotype. The problem of individual variation in the ontogenetic trajectory and overall venom phenotype is compounded if there is indeed a terminal phenotype, beyond which further growth is no longer accompanied by shifts in venom. In a venom census, terminal phenotypes appear to represent sudden change followed by a cease in change, but whether or not apparent changes are gradual from neonate to subadult and adult sizes is difficult to determine. Time-series data on captive individuals that span small juvenile to large adult sizes would control for extraneous factors and individual variation to allow a detailed look at the trajectory of venom compositional change.
The timber rattlesnake (Crotalus horridus) and eastern diamondback rattlesnake (Crotalus adamanteus) are heavy-bodied pit vipers native to the southeastern United States that reach adult body sizes over 1 m in length. Both species undergo ontogenetic shifts in venom composition [16,18,30]. We profiled the shift in venom composition across the lifespan of three individual eastern diamondback rattlesnakes (C. adamanteus) and five timber rattlesnakes (Crotalus horridus) including representatives of Type A (neurotoxic) and Type B (hemorrhagic and hemotoxic) phenotypes [38]. By assessing changes in proportion of venom proteins and shifts in chromatographic peak abundances within individuals over time, we characterize the mode and tempo of ontogenetic changes in venom of two rattlesnake species. Additionally, we investigated correlations between changes in venom composition and serum testosterone to explore the hormone’s potential as a mediator of ontogenetic venom shifts, as testosterone expression also exhibits ontogenetic changes with body size in pit vipers [39,40]. Under the hypothesis that the shift in venomic composition happens discreetly at a distinct size [30], we would predict that the majority of the magnitude of compositional change should occur over a brief window of time. In contrast, the gradual change hypothesis predicts slower continuous change in composition as snakes grow from small juveniles to large adults, possibly with a plateau when a terminal adult phenotype is reached. Finally, when venom ontogeny is adaptive, we would expect shifts in snake diet composition that are similar in mode to those of venom. We examined this possibility by mining diet data for both species from previously published studies of gut contents of wild specimens and determine the mode of ontogenetic shifts in diet as snakes grow.

2. Results & Discussion

2.1. Venom Composition Changes Slowly over Time

We detected statistically significant ontogenetic shifts in overall venom composition in all individual snakes in this study using permutational mutivariate analysis of variance (PERMANOVA) (Figure 1, Table 1). Time explained between 34 and 80% of multivariate venom variation in these snakes. Time series plots of Aitchison’s compositional distance from the small juvenile venom phenotype (Figure 2) to that of latter samples revealed ontogenetic change that occurs incrementally over years as snakes grew in body size, where shifts sometimes reached a plateau as snakes became very large adults (e.g., Figure 2, snake KW1576). We classified the mode of ontogenetic shifts as linear, sigmoid-gradual (occurs over more than the time of a single active season based on slope coefficient), and sigmoid-discrete (occurs in less than the time of a single active season), we found that the linear model was the top model in four of the eight snakes examined, although the linear model was significantly better than the sigmoid model (DeltaAIC > 2) in only a single case (Type B C. horridus KW1594; Table 2). The sigmoid model was significantly better than the linear model in the other four individuals (Table 2). Of these four cases, only one was classified as showing a rapid, discrete shift (Type B C. horridus MM0027; 242 days); the other three cases showed evidence of a gradual shift that would take longer than a single active season, although KW1753 (Type A C. horridus) was more rapid than others (309 days; Table 2). Similarly, the linear model was the best fit (but not significantly better than sigmoid, DeltaAIC = 0.6) for MM0114 (C. adamanteus), but this individual also showed a rapid change (213 days; Table 2). The tempo for the other five individuals was substantially longer (≥677 days; Table 2). Overall, we found evidence for gradual shifts in six of the eight snakes tested (four of which were more than one year in length) with MM0114 and MM0027 (and, to a lesser extent, KW1753) showing evidence of discrete ontogenetic shifts. Some unique characteristics of these two snakes are worth noting. MM0114 is unique among the three sampled C. adamanteus individuals studied herein, in that venom sampling ended before the snake reached 100 cm SVL. Margres et al. [30] estimated that juvenile and adult venoms can be differentiated in wild snakes best at 100 cm SVL, and so a significant portion of this animal’s ontogenetic change may have been yet to come at the conclusion of sampling. Finally, MM0014 is from a different population than the other C. adamanteus, so shifts could occur differently in different populations. Snake MM0027 differs from the seven other snakes in this study in that it was the largest individual at time of first sampling; it was wild caught and first sampled at 64 cm in length. Based on the profiles of the two other Type B C. horridus (which were first sampled when only half the size of MM0027 at capture), we likely missed between ⅓ and ½ of the ontogenetic shift experienced by MM0027, suggesting that earlier sampling may have led to detecting a gradual shift in this animal.
We used non-metric multidimensional scaling of binned high performance liquid chromatography (HPLC) retention time blocks to allow between-individual comparison and visualize the position of all samples from all snakes in a single multivariate space (Figure 3). Such comparisons were not possible before binning as homologous peaks were called within, but not between, individuals. Although sample size limited our ability to statistically compare movement among individuals through venom space, several patterns are apparent. First, C. adamanteus and both venom types of C. horridus occupied unique areas of venom phenotype space, and the major directions of within-individual variation in each species (or venom type) tend to have a similar pattern of orientation. For example, within-individual variation in the three C. adamanteus is oriented largely along NMDS axis 2, whereas type B C. horridus venom varies along NMDS axis 1. Meanwhile, the two type A C. horridus vary over an apparently smaller overall space, and that variation is oriented diagonally along Figure 3, and therefore lies along both axes. In seven of eight individuals, we again see apparent effects of ontogeny driving within-individual variation, as the first and last samples taken in time tend to polarize each individual in the NMDS space (beginning and end of blue arrows in Figure 3). The exception is the KW1425, which was wild-caught at moderate size (61 cm SVL). The first sample taken from this snake occupies a unique area of venom space, whereas the second and third samples show the same trend as seen in the other seven individual snakes, in that they are the two samples with the highest NMDS axis 2 values, and so lie furthest away from the samples taken as a large adult. The first sample taken from KW1425 could be therefore be biased by uncontrolled impacts of initial capture and prior free-living.
The clear implication of these analyses is that ontogenetic shifts in the complete venom profile are not wholly discrete processes. Instead, the change in composition from a neonate to adult venom tends to occur across a longer period of time. It is important first to acknowledge our small overall sample size of eight individual snakes (two species) as a limitation in our ability to generalize our findings as representative of each species, and particularly to rattlesnakes as a group. However, we are confident in the power of our repeated samples to provide useful information about the individual with high resolution, and similar gradual shifts across those individuals over which we had the most experimental control encourages the view that gradual ontogenetic shifts are commonplace and worth further exploration. The passing of time, which is likely a proxy for increases in body size, led to increasing Aitchison distance from the neonate baseline phenotype with 90% of observed change being completed in 213–989 days (Table 2). Post-hoc model comparison of time vs. body size as predictors of venom compositional change revealed that body size was a better predictor of venom change than time in five of the eight snakes studied here; body size and time explained a comparable amount of overall variation in the models (Supplemental Table S1). The result that body size can be a better predictor than time suggests a direct role of growth in venom compositional shifts despite the fact that body size measurements were subject to clear observer error associated with measuring an unanesthetized restrained snake (visible in how some SVLs appeared to decrease between subsequent measurements in Figure 2). Most samples appear to reach a plateau in amount of differentiation from the neonate venom phenotype, indicating that a terminal adult phenotype does exist; the path toward that phenotype is generally (but not always) achieved through gradual shifts in protein composition that eventually reach a plateau where change drastically slows or ceases. The power of body size as a covariate for venom compositional change is encouraging for researchers studying free-ranging snake populations, as age cannot be determined directly but body size can be easily measured.
The slow tempo of compositional change over time and body size in at least half of our samples challenges the notion that ontogenetic shifts occur discretely in C. adamanteus, as proposed by Margres et al. [30]. There are also multiple plausible biological explanations for the distinction between this study and that of Margres et al. [30]. First, most studies focused on the evolution of venom categorized individuals as being juvenile or adult based on predetermined size cutoffs in terms of total SVL [31,34]. The issue with this practice is that it inherently implies that a sudden ontogenetic shift in composition occurs once the size cutoff between designations is reached, and thus omits the potential gains of treating body size as a continuous predictor of venom variation. In our study, we leverage continuous measurement of size and sampling of venom variability. Second, any individual or geographic variation in the timing or magnitude of the ontogenetic shift would further hinder attempts to discern the tempo of shifts in census-type studies. Our study collected time-series based data over the course of several hundreds of days from individuals in captivity and controlled for unknown factors such as geography, diet, and sex by conducting analyses within single individuals. Maintaining specimens in this manner permitted concrete observation of the tempo and mode at which changes were occurring. Margres et al. [30] used a fairly restrictive definition to delineate support for the hypotheses of gradual vs. discrete shifts in multivariate venom composition, namely whether or not a sigmoid fit outperforms a linear fit of body size to venom composition across the full range of body sizes sampled from free-ranging snakes. We show here that sigmoid best fits are expected to explain overall venom composition, particularly if a terminal phenotype is eventually reached but the snake continues to grow. These sigmoid patterns of venom variability across lifespan can still take quite long to complete, which would not be obvious from sampling only wild-caught individuals a single time.

2.2. Peak-by-Peak Analysis Reveals Both Modes of Expression Variation

Gradual shifts in overall venom composition could result from many gradually shifting peaks, or multiple discretely shifting individual peaks whose shifts occur at different time points. Our repeated sampling within the same individuals permitted peak-by-peak analysis of compositional shifts, allowing an estimate of the overall percentage of venom proteins that shift expression levels and determine if both gradual and discrete modes of expression variation operate among venom genes. These analyses showed that ontogenetic shifts involved large numbers of venom proteins, with between 55% and 89% of the peaks changing significantly in relative abundance over time (Table 3, peaks marked with asterisks in Figure 4). We found evidence for both gradual and discrete modes of shift among the peaks within individual snakes (Table 3; Figure 4, Figure 5 and Figure 6, Supplemental Figures S9–13). On average, 25% of peaks showed discrete shifts. A discrete shift is exemplified by Peak 18 in Figure 6, which shifted from a low abundance juvenile state to a maintained high abundance in adulthood over the course of three sampling periods spanning 221 days. Some of the gradual shifts could be artifacts of our use of chromatographic peaks as measures of abundance since a single peak may contain multiple proteins. It is therefore possible that multiple discrete shifts spaced over time could generate an overall pattern of gradual change in peak relative abundance. We consider this unlikely to explain the bulk of our observations, given the fact that the observed transitions in relative abundance of many peaks are markedly smooth in the nature of shifts (e.g., Figure 4, peak 18, Figure 5, peak 19) and a given peak usually contains no more than 3 protein isoforms [10]. Future work that combines transcriptomic and mass-spectrometric characterization of these changes could rule out varied peak contents as an artifact and reveal whether the mode and tempo of shifts varies among venom protein families.
The repeated-samples nature of our study allowed for the control necessary to call peaks in a highly standardized fashion for each individual, providing the resolution to visualize slow shifts in peak abundance over time and classify their mode of shift based on our statistical framework. Previously, proteomic analysis of C. adamanteus and C. horridus described by Wray et al. [10] identified HPLC peak contents based on shared peptide evidence. Their analyses included individuals or siblings included in the present study, and therefore permit us to infer venom protein contents of peaks in this study (Supplemental Tables S2 and S3). The protein contents of these peaks indicate that gradual shifts occur for peaks containing multiple important venom gene families, including snake venom metalloproteinases (SVMP), snake venom serine proteases (SVSP), phospholipases A2 (PLA2), bradykinin potentiating peptides (BPP), myotoxins (MYO), and C-type lectins (CTL). While identification of specific protein contents of peaks was not a goal of our study, we emphasize the usefulness of combining serial sampling as done here with mass-spectrometric approaches for a more detailed understanding of whether different modes of shift tend to occur in different gene families. Such a pattern would suggest alternative regulatory mechanisms among gene families.
Margres et al. [41] discussed the compositional (simplex) nature of the venom phenotype extensively, pointing out that a rise in relative abundance of one protein in the venom necessitates a drop in abundance for other proteins. Within the framework of compositional data, it is important to remember that the shifts in relative abundance we see do not necessitate variation in transcription or translation of a protein across an animal’s life span. Instead, major up- or down-regulation of other venom genes can change the relative abundance of other proteins in a compositional phenotype such as venom. In other words, a relationship between time and peak relative abundance like that for KW1594 peak 11 (Figure 4; containing VEGF and NGF) can be the result of two distinct processes. First, the action of cis- or trans-regulatory mechanisms may be slowly reducing transcription over time in one or both of these genes, as post-transcriptional regulation of expression is not a large contributor to venom variation [42], leading to decreased production of this peak relative to other peaks that are produced at a constant rate. Alternatively, the reverse can be true, where upregulation of other genes via gene regulatory mechanisms is leading to a slow depression in relative abundance of these proteins in the venom, despite constant transcription and translation. It may be that those peaks displaying the steepest slopes, such as KW1594 peak 22 (Figure 4) are the most likely candidates to be targets of regulatory mechanisms over the course of development. Studies combining transcriptomics and proteomics with novel tools to study gene regulation are forthcoming in snake venomics [43,44] and will shed light on the proportion of venom that is directly regulated vs. impacted in a relative manner by changes in other venom components. More broadly, the presence of ontogenetic shifts in peak abundance that vary from highly linear to highly discrete and sigmoidal in shape suggests that alternate regulatory mechanisms are at play among the genes contributing to the same predatory phenotype of venom. Our work here therefore indicates that venom provides ostensible targets for studies of different gene regulatory mechanisms operating within a common functional module.

2.3. Assessing Testosterone as a Predictor of Venom Change

Our results do not meet the prediction of a strong linear relationship between the change in venom composition and the testosterone concentration across the course of our study (Figure 7), and a high correlation between body size and testosterone violates the assumption of a lack of multicollinearity in predictors, preventing us from conducting informative multiple regression analysis with both body size and testosterone level in the models. The change in venom composition and ln-testosterone concentration were significantly correlated in four of the six snakes for which testosterone was assayed (KW1425: R2 = 43%, p = 0.01; KW1576: R2 = 55%, p = 0.0003; KW1753: R2 = 18%, p = 0.04; KW1726: R2 = 62%, p = 0.001; MM0071: R2 = 23%, p = 0.07; KW1594: R2 = 15%, p = 0.18). These significant correlations underscore the patterns of maturation in both these snake species, where both testosterone concentrations and venom distance from the juvenile venom phenotype increase as snakes mature. We expected both venom to change and testosterone to increase into adulthood, as ontogenetic changes in venom composition are known from both species studied [18], and testosterone is known to increase with body size in snakes [39,40,45] as it mediates sexual maturity in rattlesnakes in both sexes [46]. As a result, a correlative study like ours required a specific prediction of a very strong relationship. Anything but a very strong linear relationship can be explained by other factors. Clearly, Figure 7 shows that the samples taken during periods of high testosterone are also samples that have changed in venom composition from the earliest sampling periods. However, the range in Aitchison distances associated with moderate to high testosterone levels spans nearly the full spectrum of the ontogenetic shift in venom. While seasonal testosterone shifts are known from adult rattlesnakes [46,47,48,49,50], there is little information on seasonal testosterone in young, small individuals [39]. Therefore, we suspect that the patterns in Figure 7 are best explained by slow ontogenetic shifts in venom that clearly occur over the snakes’ lifespans (Figure 2) and that these are eventually accompanied by changes in testosterone as snakes reach a size of sexual maturity.
The proximate physiological causes of intra-individual venom variation have received very little attention, despite the potential of such information as a tool in both studies of venom gene regulation and in proper care of animals used in venom production for medical uses. Our investigation of testosterone as a possible physiological trigger of venom variation is among the first to assess the proximate factors inducing the extreme changes seen in animal venoms [12,51]. Generally, our lack of power to draw stronger conclusions about the potential for a relationship between venom and testosterone reflects how our study, and any correlative study that involves hormones with complex roles, will be limited compared to an experimental manipulation of hormone levels. We emphasize that our inconclusive testosterone results are thus a demonstration of how experimental manipulations, instead of purely correlative studies, are the key to uncovering the proximate causes of ontogenetic shifts in venom. Experimental manipulation of levels of several metabolic and maturation hormones is now commonplace in reptiles [29,39]. Future studies that manipulate hormone levels in neonate snakes and conduct before and after sampling of venom will be particularly valuable in uncovering the triggers of shifts in venom phenotypes as snakes grow.

2.4. Mode of Dietary Shifts

We found evidence that ontogenetic changes in diet reported in prior studies occurred across body size in C. horridus (n = 179 total prey counted), where there were significant differences in diet composition between both the 2nd and 3rd body size quartiles (p = 0.006), as well as between the 3rd and 4th (chi-square test; p = 0.05, Figure 8). The shift of diet in C. horridus manifests as consumption of shrews and mice at small sizes, with gradual increases in the abundance of larger mammals such as squirrels and rats at larger body sizes. Inspecting the shifts in abundance of the individual prey categories suggests that alternative modes of dietary shifts may occur within C. horridus. Rat and vole abundances show relatively steady increases over time, while mouse abundances shift in a discrete fashion between the 2nd and 3rd body size quartiles. The overall nature of the shift in C. adamanteus is more discrete in nature (Figure 8), with significant compositional variation between only the 3rd and 4th body size quartiles (p = 0.002), although failure to detect differences among smaller size classes could be due to low power (n = 58 total prey counted with particularly low numbers in the smallest size class). The 4th size quartile present in the diet study actually lies outside the body size range of the snakes whose venom we sampled in this study, meriting examination of the taxon-specific patterns of abundance variation in the 1st through 3rd size quartiles. There is a propensity toward the consumption of larger prey taxa over time in C. adamanteus, evidenced by the C. adamanteus diet exclusively (100%) consisting of mice and rats at the 512–782 mm range whereas rabbits and squirrels appear in the 2nd size quartile (783–1053 mm), and rabbits double in abundance between the 2nd and 3rd quartiles.
As dietary breadth has been previously shown to impact venom composition in taxa such as Echis vipers and Sistrurus rattlesnakes [8,12], the specific nature of dietary shifts such as those in C. adamanteus and C. horridus may be expected to predict the mode and tempo of venom ontogeny, but like the mode and tempo of venom shifts, the nature of dietary variation is complex, with the overall nature of these shifts underscored by apparent gradual and discrete shifts in the abundance of individual prey types. Our ability to link dietary shifts to the mode of venom shifts is likely limited by the differential scope of the data collected, as we compare range wide diet data with venom sampled from individuals from a single region in the southeastern U.S. Given the evolutionary lability of snake venoms such as those of C. adamanteus [30], it may be that local diet studies may provide an even stronger match to the mode and tempo of venom ontogeny. High resolution diet studies now exist for several species of vipers [52,53,54], and these studies often include information about the body size of snakes that had consumed particular prey items. Notably, neutrality or near-neutrality has not been rejected on a protein-by-protein basis for any venom, so functional studies of venom ontogeny’s impact on prey toxicity will be required to clearly demonstrate the hypothesized adaptive nature of venom ontogeny as a phenotypic match to a changing diet. Future studies that combine diet data with profiling of the mode and tempo of ontogenetic shifts in the venom of each snake species would allow a phylogenetic comparative test of the hypothesis that the ontogenetic shifts in venom evolve to match the nature of a snake species’ movement through a dietary ecological niche. Additionally, treating the presence or strength of ontogenetic shifts as a trait for population-level comparisons may help differentiate adaptive and neutral explanations for ontogenetic venom variation among venom proteins [32].

3. Conclusions

Ontogenetic shifts in snake venom are more the rule than exception, changing the functional profile of venom in both biologically and biomedically important ways [22,27,55,56]. Our study is the first to undertake venom compositional profiling at the level of an individual spanning neonate to large adult sizes with a high density of repeated samples, which allowed several novel insights into the nature of the ontogenetic shifts in the species studied. In particular, we show multiple shifts in venom composition that occurred gradually over the course of over 1 m, or multiple years, of growth. The vast majority of previous work on ontogenetic shifts in venom represent population sampling efforts, where venom traits are regressed against individual body size, and therefore leave individual and population variability as an uncontrolled variable [22]. Studies such as those are used in documenting the presence and intensity of a shift, but the presence of extensive individual variability in venoms [36] leads to poor inference into the specific nature of ontogeny. Still, many studies of venom in both natural and laboratory settings merely report venoms as juvenile or adult, and thus bin highly continuous variation in venom traits over lifespan into groupings of binary predictor variables.
We therefore suggest that one of the most important implications of our work is the need to record more detailed life history data, particularly body size, as venom is collected. The only comparably long-term study to ours involving another large pit viper (Bothrops atrox) provides a demonstration of the utility of considering body size information when making inferences about venom ontogeny [57]. These authors captured snakes in the wild classified as adults and used HPLC and enzymatics to profile change in venom during long-term captivity. While minor compositional shifts were reported in all snakes, only 4 of 13 snakes showed shifts in venom proteins classified as “core function toxins”, and the authors provided several potential explanations for the major changes in these select individuals, including dietary shifts, collection during the rainy season, and that two of the four snakes were the two smallest females collected. Repeated measurements of body size were not provided in the study by Amazonas et al. [57], but three of the four animals showing the large compositional shifts were indeed three of the four animals with the smallest body sizes at collection. Therefore, extrapolating our results in Crotalus to Bothrops leads us to hypothesize that these animals may have simply had the most opportunity for greater compositional change in venom, whereas the larger animals were already close to their terminal adult venom phenotype at first collection.
Animal venom provides a unique window into the path from genotype to phenotype to function, as the trait is broadly accessible through modern ‘omic’ and enzymological technologies [9,58]. Venomous snakes in particular are accompanied by extensive studies of their behavioral [59], physiological [39,60,61], and dietary ecology [22,52,54]. Venom systems therefore present an incredible opportunity for integrative biology that addresses both the proximate and ultimate causes of variation in a key functional trait, as long as venom research is done in conjunction with collecting quality life-history information. Shifts in venom are likely adaptive responses to shifts in diet that are natural consequences of gape-limited predation in snakes. Closely related species of venomous animals that show extreme variation in the range of body sizes experienced and the steepness and intensity of ontogenetic shifts in diet represent fruitful ground for studying the genomic basis and evolution of ontogenetic shifts in venom, and will inform us further on the interplay between the evolution of development and the evolution of diversity.
Data accessibility: raw data are provided as supplemental information.

4. Materials and Methods

4.1. Acquisition of Specimens and Samples

Venom samples (n = 148 total samples) were obtained from three C. adamanteus and five C. horridus. Individuals were either originally caught in Florida or Georgia (Supplemental Figure S1) or raised from previously established cohorts at Florida State University; detailed specimen information is provided in Table 4. These snakes were maintained in captivity for multiple years, growing from juvenile size (newborn or wild caught juveniles) to large adult sizes. Collection and handling of these animals was approved under Florida Fish and Wildlife Conservation Commission permits LSSC-13-00004, LSSC-09-0399, and LSSC-12-00071A and Institutional Animal Care and Use Committee (IACUC) protocols #0924 and #1333. We collected venom and blood samples every second month by inducing snakes to bite through a parafilm-covered beaker, and venom was then lyophilized and stored at −80 °C until use. Blood samples were centrifuged at 1500× g for 10 min, and serum was stored at −40 °C.

4.2. High-Performance Liquid Chromatography

Venom samples were redissolved in LCMS water and diluted to a concentration of 0.18 mg/mL. We determined venom protein concentration by the A280 reading from a NanoDrop Photospectrometer (ThermoFisher Scientific, Waltham, MA, USA). We then loaded 15 μg of venom protein onto a Prominence HPLC apparatus (Shimadzu Corp., Columbia, MD, USA). Samples were run through a Phenomenex Aeris widepore column (3.6 micron, C18, part # 00G-4482-AN) maintained at 25 °C. Samples were eluted with solution A as 0.1 mM TFA in water and solution B as 0.06 mM TFA in acetonitrile, and absorbance at 220 nm was measured over 150 min with the following run parameters: 5 min at 10% B, 110 min increasing from 10 to 55% B, 5 min increasing from 55 to 75% B, 5 min at 75% B, 5 min decreasing from 75–10% B, and 5 min at 10% B. Venom peak areas were quantified using Shimadzu Lab Solution v2 software, combining default peak calls by the iPeakFinder algorithm with manual curation of the peak baseline and filtering of peaks making up less than 0.1% total area as noise. Peaks were called uniquely per individual snake, by aligning the serial HPLC profiles from all samples from a single individual and calling peaks that were present. The relative abundance of HPLC peaks in each sample will be henceforth referred to as venom composition.

4.3. Testosterone ELISA

We used a commercially available EIA kit [62] to measure serum testosterone. Because we anticipated variation between adult and neonate testosterone would differ more than between species, we validated the kit for us in adult and neonate rattlesnakes by pooling serum of both species for adults and neonates separately. Both adult and neonate pools were serially diluted and assessed for parallelism against results from a serial dilution of testosterone at known concentrations to verify the kit’s accuracy across a range of concentrations for adult and neonate serum. We determined the optimal dilution for assessing both age classes simultaneously to be 9:100 unextracted serum to assay buffer. Samples were extracted according to the protocol provided by the kit manufacturer [62]. Briefly, 12 μL serum was extracted by vortexing with diethyl ether at a 5:1 ether:sample ratio, frozen in a dry ice and ethanol bath, and the remaining ether solution pipetted into a new tube. This was repeated for each sample for maximum efficiency. The pooled ether-extractions for each sample were dried in a speedvac, stored at −20 overnight, and redissolved in 133 μL of assay buffer. Samples were plated and treated according to the kit manufacturer protocol [62]. The optical density of each well was read at 450 nm by a SpectraMax iD5 multi-mode microplate reader (Molecular Devices); results are reported in ng/mL. Samples were haphazardly assigned to each plate. There was a subset of samples that did not have sufficient serum to run in duplicate (n = 9); for these samples, 6.5 μL of serum was extracted as above and redissolved in 71.5 μL of assay buffer and were run singly. All other samples were run in duplicate, and the values of each duplicate averaged for use in analysis. Samples falling below the limit of detection (n = 5) were set at the lower limit of detection for analysis (0.001 ng/mL). The mean intra-assay coefficient of variation was 1.96 ± 1.74% and the inter-assay coefficient of variation was 4.7%.

4.4. Dietary Shifts

We accessed previously published diet studies for both C. horridus [63] and C. adamanteus [64] to assess the extent to which these snakes exist as size-structured populations accessing different resources. Furthermore, each study included SVL of the snakes containing diet items, allowing us to assess the extent to which diet composition changes suddenly or more gradually over time. To do this, we binned the SVLs of both species into quartiles based on the shortest and longest snake in the diet datasets, and tabulated counts of prey in the diet within a given size quartile. Both snakes are known to specialize on mammals, and we were able to categorize the mammals presented in the diet studies into general categories that varied by size [63,64]. In order of smallest to largest average prey size, C. horridus fed on shrews, mice, voles, chipmunks, rats, squirrels, and rabbits, while C. adamanteus fed on only mice, rats, squirrels, and rabbits. We tested for differences in the overall proportion of these prey classes between adjacent size class quartiles using chi-square tests.

4.5. Statistical Analyses

All statistical analyses were performed in R v. 3.6.0 [65]. Venom phenotype as measured using HPLC represents compositional data, and thus exists in the simplex [30,66] and requires transformation prior to analysis. The presence of an ontogenetic shift in the whole venom HPLC profiles of each individual snake was tested with permutational MANOVA as implemented in the adonis function of the vegan 2.5 R package [67], with isometric-log-ratio (ilr) transformed venom composition as a multivariate response and day of study (i.e., time) as a continuous independent variable. To determine the mode of the ontogenetic shift, we visualized venom compositional change by calculating differences between each venom sample and the mean composition of the three earliest juvenile samples taken from a given individual. In this way, our compositional distances for each time point were relative to an average of early composition and not relative to a potentially noisy single first sample. Relative abundance values of zero within a sample were replaced by fractional values using the cmultRepl function in the zcompositions R package [68]. Next, the aDist function in zCompositions was used to calculate the Aitchison compositional distance between the focal sample and the juvenile sample average [68]. The Aitchison distance is a Euclidean distance between data that has been previously subjected to centered log–ratio (clr) transformation, and is an appropriate measure of distance between compositional datasets [66]. Gradual modes of ontogenetic shift are expected to be best fit by a linear model or a sigmoid model with a gradual slope. We used the lm function from base R and the nlsLM function from the minpack.lm package [69] to fit linear regressions and sigmoid curves, respectively, to the relationship between mean-scaled Aitchison distance and time for each snake. We then compared the fit of each model using AICc values. If a nonlinear sigmoid curve was a better fit than a linear model, we used the slope at the sigmoid inflection point (the scale parameter) to determine if 90 percent of the compositional change to the asymptote estimate occurred in less than one active season (number of days from April through the end of November: 244 days). If the change occurred in less than one active season, we counted the peak as having a discrete shift mode, whereas non-linear changes that took longer were counted as gradual. Margres et al. [30] defined support for a discrete shift as any case where a sigmoid curve fits the data better than a linear relationship, and thus our current work is differentiated from Margres et al. [30] in recognizing that a change manifesting as a sigmoid relationship may still take a long time from initiation to asymptote. Last, we placed all samples from all snakes in a comparable multivariate space for visualization in the same analysis of phenotype. To do this, we binned peak areas by retention time eleven retention time windows: R1 (10–20 min), R2 (20–30 min), R3 (30–40 min), R4 (40–50 min), R5 (50–60 min), R6 (60–70 min), R7 (70–80 min), R98 (80–90 min), R9 (90–100 min), R10 (100–110 min), and R11 (110–125 min). Next, we clr-transformed the percent area data for these eleven windows across all 148 venom samples used in this study. We then calculated pairwise Aitchison distance between all samples, and projected them into a 2 dimensional multivariate venom space using NMDS implemented in the metaMDS function in the vegan v.2.5 [67] package in R.
Upon detecting significant shifts in the overall venom profile, we followed with modeling the clr transformed abundance of each peak versus time. This peak-by-peak analysis provided a post-hoc test of which peaks varied in abundance significantly over the course of the study and their mode of change, where mode (no shift, linear shift, gradual sigmoid, or discrete sigmoid) was determined as above using linear and non-linear regression. Non-linear fits were again obtained using the nlsLM function from the minpack.lm package [69].
Last, we determined whether testosterone concentrations were closely associated with the change in venom composition as snakes grew via linear regression. Given the correlative nature of these tests, the best evidence for a potential role of testosterone would be a particularly strong linear relationship between the change in testosterone concentration between the focal and initial reads and the change in venom composition as measured by Aitchison distance. Other forms of a relationship would be inconclusive, as testosterone is expected to be higher in adult than juvenile snakes.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2072-6651/12/10/659/s1, Figure S1: All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus adamanteus—KW1425, Figure S2: All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus adamanteus— KW1726, Figure S3: All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus adamanteus—MM0114, Figure S4: All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus horridus—KW1576 (Type B venom), Figure S5, All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus horridus—KW1594 (Type B venom), Figure S6: All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus horridus—MM0027 (Type B venom), Figure S7: All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus horridus—KW1753 (Type A venom), Figure S8: All HPLC chromatograms used to profile the ontogenetic shifts in the venom of Crotalus horridus—MM0071 (Type A venom), Figure S9: Peak-by-peak analysis of relative abundance shifts over time for an individual Crotalus adamanteus, Figure S10: Peak-by-peak analysis of relative abundance shifts over time for an individual Crotalus adamanteus, Figure S11: Peak-by-peak analysis of relative abundance shifts over time for an individual type B Crotalus horridus, Figure S12: Peak-by-peak analysis of relative abundance shifts over time for an individual type B Crotalus horridus, Figure S13: Peak-by-peak analysis of relative abundance shifts over time for an individual type A Crotalus horridus, Table S1: Comparison of time (day) and body size (SVL) as predictors of Aitchison distance of venom samples from the juvenile average venom composition. Bold AIC scores indicate the top model for that individual, Table S2: Protein contents of peaks based on HPLC profiles and mass-spectrometry results presented in Wray et al. [10], Table S3: Protein contents of peaks based on HPLC profiles and mass-spectrometry results presented in Wray et al. [10].

Author Contributions

Conceptualization, M.L.H., K.W., J.M., M.J.M. and D.R.R.; methodology, M.L.H., N.M.C., K.W., J.M., M.J.M. and D.R.R.; formal analysis, E.M.H., R.B.S., M.L.H. and T.J.C.; investigation, E.M.H., R.B.S., M.L.H., N.M.C., S.A.E. and M.P.H.; data curation, E.M.H., R.B.S., N.M.C., S.A.E. and M.P.H.; writing—original draft, E.M.H., R.B.S., M.L.H. and N.M.C.; visualization, E.M.H., R.B.S. and M.L.H.; supervision, R.B.S. and D.R.R.; writing—review & editing, M.L.H., N.M.C., S.A.E., M.P.H., M.J.M., T.J.C. and D.R.R.; funding acquisition, M.L.H. and D.R.R.; resources, M.P.H., K.W., J.M., M.J.M., T.J.C. and D.R.R. All authors have read and agreed to the published version of the manuscript.

Funding

This project was supported by the Undergraduate Research Opportunity Program (UROP) at Florida State University, NSF DEB 1145987 and NSF DEB 1638902 to DRR, and NSF PRFB 1711141 to MLH.

Acknowledgments

Carl Whittington provided invaluable assistance with HPLC analyses. Special thanks go to the students of the Reptile and Amphibian Evolution Reading Group. The authors thank Dick Barlett, Mark S. Margres, and Jason Strickland for assistance in the field and Megan Lamb, Ethan Bourque, Jennifer Wanat, and Rebecca Bernard with the Florida DEP and Apalachicola River NERR for access to field sites.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cabej, N.R. Epigenetic Principles of Evolution; Elsevier Inc.: Amsterdam, The Netherlands, 2012. [Google Scholar] [CrossRef]
  2. Werner, E.E.; Gilliam, J.F. The Ontogenetic Niche and Species Interactions in Size-Structured Populations. Annu. Rev. Ecol. Syst. 1984, 15, 393–425. [Google Scholar] [CrossRef]
  3. Claessen, D.; Dieckmann, U. Ontogenetic niche shifts and evolutionary branching in size-structured populations. Evol. Ecol. Res. 2002, 4, 189–217. [Google Scholar]
  4. Rowe, L.; Ludwig, D. Size and Timing of Metamorphosis in Complex Life Cycles: Time Constraints and Variation. Ecology 1991, 72, 413–427. [Google Scholar] [CrossRef]
  5. Graham, B.S.; Grubbs, D.; Holland, K.; Popp, B.N. A rapid ontogenetic shift in the diet of juvenile yellowfin tuna from Hawaii. Mar. Biol. 2007, 150, 647–658. [Google Scholar] [CrossRef]
  6. Cipriani, V.; Debono, J.; Goldenberg, J.; Jackson, T.N.; Arbuckle, K.; Dobson, J.; Koludarov, I.; Li, B.; Hay, C.; Dunstan, N.; et al. Correlation between ontogenetic dietary shifts and venom variation in Australian brown snakes (Pseudonaja). Comp. Biochem. Physiol. Part C 2017, 197, 53–60. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Barton, K.E.; Boege, K. Future directions in the ontogeny of plant defence: Understanding the evolutionary causes and consequences. Ecol. Lett. 2017, 20, 403–411. [Google Scholar] [CrossRef]
  8. 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 Boil. Sci. 2009, 276, 2443–2449. [Google Scholar] [CrossRef] [Green Version]
  9. Calvete, J.J.; Sanz, L.; Cid, P.; De La Torre, P.; Flores-Díaz, M.; Dos Santos, M.C.; Borges, A.; Bremo, A.; Angulo, Y.; Lomonte, B.; et al. Snake Venomics of the Central American RattlesnakeCrotalus simusand the South AmericanCrotalus durissusComplex Points to Neurotoxicity as an Adaptive Paedomorphic Trend alongCrotalusDispersal in South America. J. Proteome Res. 2010, 9, 528–544. [Google Scholar] [CrossRef]
  10. Durban, J.; Pérez, A.; Sanz, L.; Gómez, A.; Bonilla, F.; Rodríguez, S.; Chacón, D.; Sasa, M.; Angulo, Y.; Gutiérrez, J.M.; et al. Integrated “omics” profiling indicates that miRNAs are modulators of the ontogenetic venom composition shift in the Central American rattlesnake, Crotalus simus simus. BMC Genom. 2013, 14, 234. [Google Scholar] [CrossRef] [Green Version]
  11. Durban, J.; Sanz, L.; Trevisan-Silva, D.; Neri-Castro, E.; Alagón, A.; Calvete, J.J. Integrated Venomics and Venom Gland Transcriptome Analysis of Juvenile and Adult Mexican Rattlesnakes Crotalus simus, C. tzabcan, and C. culminatus Revealed miRNA-modulated Ontogenetic Shifts. J. Proteome Res. 2017, 16, 3370–3390. [Google Scholar] [CrossRef]
  12. Gibbs, H.L.; Sanz, L.; Chiucchi, J.E.; Farrell, T.M.; Calvete, J.J. Proteomic analysis of ontogenetic and diet-related changes in venom composition of juvenile and adult Dusky Pigmy rattlesnakes (Sistrurus miliarius barbouri). J. Proteom. 2011, 74, 2169–2179. [Google Scholar] [CrossRef] [PubMed]
  13. Jackson, T.N.W.; Koludarov, I.; Ali, S.A.; Dobson, J.; Zdenek, C.N.; Dashevsky, D.; Brouw, B.O.D.; 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. [Google Scholar] [CrossRef] [PubMed]
  14. Mackessy, S.P.; Sixberry, N.M.; Heyborne, W.H.; Fritts, T. Venom of the Brown Treesnake, Boiga irregularis: Ontogenetic shifts and taxa-specific toxicity. Toxicon 2006, 47, 537–548. [Google Scholar] [CrossRef] [PubMed]
  15. Madrigal, M.; Sanz, L.; Flores-Díaz, M.; Sasa, M.; Núñez, V.; Alape-Girón, A.; Calvete, J.J. Snake venomics across genus Lachesis. Ontogenetic changes in the venom composition of Lachesis stenophrys and comparative proteomics of the venoms of adult Lachesis melanocephala and Lachesis acrochorda. J. Proteom. 2012, 77, 280–297. [Google Scholar] [CrossRef]
  16. Rokyta, D.R.; Margres, M.J.; Ward, M.; Sanchez, E.E. The genetics of venom ontogeny in the eastern diamondback rattlesnake (Crotalus adamanteus). PeerJ 2017, 5, e3249. [Google Scholar] [CrossRef] [Green Version]
  17. Underwood, A.H.; Seymour, J. Venom ontogeny, diet and morphology in Carukia barnesi, a species of Australian box jellyfish that causes Irukandji syndrome. Toxicon 2007, 49, 1073–1082. [Google Scholar] [CrossRef]
  18. Wray, K.P.; Margres, M.J.; Seavy, M.; Rokyta, D.R. Early significant ontogenetic changes in snake venoms. Toxicon 2015, 96, 74–81. [Google Scholar] [CrossRef]
  19. Zelanis, A.; Ventura, J.D.S.; Chudzinski-Tavassi, A.M.; Furtado, M.D.F.D. Variability in expression of Bothrops insularis snake venom proteases: An ontogenetic approach. Comp. Biochem. Physiol. Part C 2007, 145, 601–609. [Google Scholar] [CrossRef]
  20. McElroy, T.; McReynolds, C.N.; Gulledge, A.; Knight, K.R.; Smith, W.E.; Albrecht, E.A. Differential toxicity and venom gland gene expression in Centruroides vittatus. PLoS ONE 2017, 12, e0184695. [Google Scholar] [CrossRef] [Green Version]
  21. Casewell, N.R.; Wüster, W.; Vonk, F.; Harrison, R.A.; Fry, B.G. Complex cocktails: The evolutionary novelty of venoms. Trends Ecol. Evol. 2013, 28, 219–229. [Google Scholar] [CrossRef]
  22. Mackessy, S.P. Venom Ontogeny in the Pacific Rattlesnakes Crotalus viridis helleri and C. v. oreganus. Copeia 1988, 1988, 1445927. [Google Scholar] [CrossRef]
  23. Gutiérrez, J.M.; Avila, C.; Camacho, Z.; Lomonte, B. Ontogenetic changes in the venom of the snake Lachesis muta stenophrys (bushmaster) from Costa Rica. Toxicon 1990, 28, 419–426. [Google Scholar] [CrossRef]
  24. Borja, M.; Neri-Castro, E.; Pérez-Morales, R.; Strickland, J.L.; Ponce-López, R.; Parkinson, C.L.; Espinosa-Fematt, J.A.; Sáenz-Mata, J.; Flores-Martínez, E.; Alagón, A.; et al. Ontogenetic Change in the Venom of Mexican Black-Tailed Rattlesnakes (Crotalus molossus nigrescens). Toxins 2018, 10, 501. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Casewell, N.R.; Wagstaff, S.C.; Wüster, W.; Cook, D.A.N.; Bolton, F.M.S.; King, S.I.; Pla, D.; Sanz, L.; Calvete, J.J.; Harrison, R.A. Medically important differences in snake venom composition are dictated by distinct postgenomic mechanisms. Proc. Natl. Acad. Sci. USA 2014, 111, 9205–9210. [Google Scholar] [CrossRef] [Green Version]
  26. 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. [Google Scholar] [CrossRef]
  27. 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. [Google Scholar] [CrossRef] [PubMed]
  28. Calvete, J.J.; Sanz, L.; Pérez, A.; Borges, A.; Vargas, A.M.; Lomonte, B.; Angulo, Y.; Gutiérrez, J.M.; Chalkidis, H.M.; Mourão, R.H.; et al. Snake population venomics and antivenomics of Bothrops atrox: Paedomorphism along its transamazonian dispersal and implications of geographic venom variability on snakebite management. J. Proteom. 2011, 74, 510–527. [Google Scholar] [CrossRef]
  29. Claunch, N.M.; Frazier, J.A.; Escallón, C.; Vernasco, B.J.; Moore, I.T.; Taylor, E.N. Physiological and behavioral effects of exogenous corticosterone in a free-ranging ectotherm. Gen. Comp. Endocrinol. 2017, 248, 87–96. [Google Scholar] [CrossRef]
  30. 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. [Google Scholar] [CrossRef] [Green Version]
  31. Holding, M.L.; Biardi, J.E.; Gibbs, H.L. Coevolution of venom function and venom resistance in a rattlesnake predator and its squirrel prey. Proc. R. Soc. B Boil. Sci. 2016, 283, 20152841. [Google Scholar] [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. [Google Scholar] [CrossRef] [PubMed]
  33. 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 Boil. Sci. 2019, 286, 20190810. [Google Scholar] [CrossRef] [PubMed]
  34. Pla, D.; Sanz, L.; Sasa, M.; Acevedo, M.E.; Dwyer, Q.; Durban, J.; Pérez, A.; Rodriguez, 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. [Google Scholar] [CrossRef] [PubMed]
  35. Minton, S.A.; Weinstein, S.A. Geographic and ontogenic variation in venom of the western diamondback rattlesnake (Crotalus atrox). Toxicon 1986, 24, 71–80. [Google Scholar] [CrossRef]
  36. Smiley-Walters, S.A.; Farrell, T.M.; Gibbs, H.L. High levels of functional divergence in toxicity towards prey among the venoms of individual pigmy rattlesnakes. Biol. Lett. 2019, 15, 20180876. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Wittenberg, R.D.; Beaupre, S.J. Growth of Timber Rattlesnakes (Crotalus horridus) in an Agriculturally Fragmented and a Contiguously Forested Habitat. Herpetologica 2014, 70, 171–183. [Google Scholar] [CrossRef]
  38. Rokyta, D.R.; Wray, K.P.; Margres, M.J. The genesis of an exceptionally lethal venom in the timber rattlesnake (Crotalus horridus) revealed through comparative venom-gland transcriptomics. BMC Genom. 2013, 14, 394. [Google Scholar] [CrossRef] [Green Version]
  39. Taylor, E.N.; DeNardo, D.F. Sexual size dimorphism and growth plasticity in snakes: An experiment on the Western Diamond-backed Rattlesnake (Crotalus atrox). J. Exp. Zool. Part A 2005, 303, 598–607. [Google Scholar] [CrossRef]
  40. Lutterschmidt, W.I.; Reinert, H.K. The Effect of Ingested Transmitters upon the Temperature Preference of the Northern Water Snake, Nerodia S. sipedon. Herpetologica 1990, 46, 39–42. [Google Scholar] [CrossRef]
  41. Margres, M.J.; McGivern, J.J.; Wray, K.P.; Seavy, M.; Calvin, K.; Rokyta, D.R. Linking the transcriptome and proteome to characterize the venom of the eastern diamondback rattlesnake (Crotalus adamanteus). J. Proteom. 2014, 96, 145–158. [Google Scholar] [CrossRef]
  42. Rokyta, D.R.; Margres, M.J.; Calvin, K. Post-transcriptional Mechanisms Contribute Little to Phenotypic Variation in Snake Venoms. G3 Genes Genomes Genet. 2015, 5, 2375–2382. [Google Scholar] [CrossRef] [Green Version]
  43. Schield, D.R.; Card, D.C.; Hales, N.R.; Perry, B.W.; Pasquesi, G.M.; Blackmon, H.; Adams, R.H.; Corbin, A.B.; Smith, C.F.; Ramesh, B.; et al. The origins and evolution of chromosomes, dosage compensation, and mechanisms underlying venom regulation in snakes. Genome Res. 2019, 29, 590–601. [Google Scholar] [CrossRef]
  44. Post, Y.; Puschhof, J.; Beumer, J.; Kerkkamp, H.M.; De Bakker, M.A.; Slagboom, J.; De Barbanson, B.; Wevers, N.R.; Spijkers, X.M.; Olivier, T.; et al. Snake Venom Gland Organoids. Cell 2020, 180, 233–247.e21. [Google Scholar] [CrossRef]
  45. King, R.B.; Cline, J.H.; Hubbard, C.J. Age and Litter Effects on Testosterone Levels in Young Water Snakes. Copeia 2000, 2000, 593–596. [Google Scholar] [CrossRef]
  46. Taylor, E.N.; DeNardo, D.F.; Jennings, D.H. Seasonal steroid hormone levels and their relation to reproduction in the Western Diamond-backed Rattlesnake, Crotalus atrox (Serpentes: Viperidae). Gen. Comp. Endocrinol. 2004, 136, 328–337. [Google Scholar] [CrossRef] [Green Version]
  47. Lind, C.M.; Moore, I.T.; Vernasco, B.J.; Farrell, T.M. Seasonal testosterone and corticosterone patterns in relation to body condition and reproduction in a subtropical pitviper, Sistrurus miliarius. Gen. Comp. Endocrinol. 2018, 267, 51–58. [Google Scholar] [CrossRef]
  48. Lutterschmidt, W.I.; Lutterschmidt, D.I.; Mason, R.T.; Reinert, H.K. Seasonal variation in hormonal responses of timber rattlesnakes (Crotalus horridus) to reproductive and environmental stressors. J. Comp. Physiol. B 2009, 179, 747–757. [Google Scholar] [CrossRef] [PubMed]
  49. Lind, C.M.; Husak, J.F.; Eikenaar, C.; Moore, I.T.; Taylor, E.N.; Lind, C.R.P. The relationship between plasma steroid hormone concentrations and the reproductive cycle in the Northern Pacific rattlesnake, Crotalus oreganus. Gen. Comp. Endocrinol. 2010, 166, 590–599. [Google Scholar] [CrossRef] [Green Version]
  50. Hoss, S.K.; Schuett, G.W.; Earley, R.L.; Smith, L.L. Reproduction in Male Crotalus adamanteus Beauvois (Eastern Diamond-Backed Rattlesnake): Relationship of Plasma Testosterone to Testis and Kidney Dimensions and the Mating Season. Southeast. Nat. 2011, 10, 95. [Google Scholar] [CrossRef]
  51. Claunch, N.M.; Holding, M.L.; Escallón, C.; Vernasco, B.; Moore, I.T.; Taylor, E.N. Good vibrations: Assessing the stability of snake venom composition after researcher-induced disturbance in the laboratory. Toxicon 2017, 133, 127–135. [Google Scholar] [CrossRef] [PubMed]
  52. Martins, M.; Marques, O.a.V.; Sazima, I. Bothrops—Ecological and Phylogenetic Correlates of Feeding Habits. Biol. Vipers. 1996, 592, 307–328. Available online: http://jararacailhoa.org/conservacaoinsularis/bothrops_feeding.pdf (accessed on 13 May 2020).
  53. Glaudas, X.; Glennon, K.L.; Martins, M.; Luiselli, L.; Fearn, S.; Trembath, D.F.; Jelić, D.; Alexander, G.J. Foraging mode, relative prey size and diet breadth: A phylogenetically explicit analysis of snake feeding ecology. J. Anim. Ecol. 2019, 88, 757–767. [Google Scholar] [CrossRef]
  54. Babb, R.D.; Bradley, G.L.; Brennan, T.C.; Holycross, A.T. Preliminary assessment of the diet of gyalopion quadrangulare (serpentes: Colubridae). Southwest. Nat. 2005, 50, 390–392. [Google Scholar] [CrossRef]
  55. Mackessy, S.P.; Williams, K.; Ashton, K.G. Ontogenetic Variation in Venom Composition and Diet of Crotalus oreganus concolor: A Case of Venom Paedomorphosis? Copeia 2003, 2003, 769–782. [Google Scholar] [CrossRef]
  56. Andrade, D.V.; Abe, A.S. Relationship of Venom Ontogeny and Diet in Bothrops; Allen Press: Lawrence, KS, USA, 1999. [Google Scholar]
  57. Amazonas, D.R.; Freitas-De-Sousa, L.A.; Orefice, D.P.; Sousa, L.F.; Martinez, M.G.; Mourão, R.H.V.; Chalkidis, H.D.M.; Camargo, P.; Moura-Da-Silva, A. Evidence for Snake Venom Plasticity in a Long-Term Study with Individual Captive Bothrops atrox. Toxins 2019, 11, 294. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Calvete, J.J. The challenge of integrating proximate and ultimate causes to reconstruct the natural histories of venoms: The evolutionary link. Expert Rev. Proteom. 2016, 13, 1059–1061. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Putman, B.J.; Coss, R.G.; Clark, R.W. The ontogeny of antipredator behavior: Age differences in California ground squirrels (Otospermophilus beecheyi) at multiple stages of rattlesnake encounters. Behav. Ecol. Sociobiol. 2015, 69, 1447–1457. [Google Scholar] [CrossRef]
  60. Holding, M.L. Short-Distance Translocation of the Northern Pacific Rattlesnake (Crotalus o. oreganus): Effects on Volume and Neurogenesis in the Cortical Forebrain, Steroid Hormone Concentrations and Behaviors; California Polytechnic State University: San Luis Obispo, CA, USA, 2011. [Google Scholar] [CrossRef]
  61. Holding, M.L.; Frazier, J.A.; Dorr, S.W.; Henningsen, S.N.; Moore, I.T.; Taylor, E.N. Physiological and Behavioral Effects of Repeated Handling and Short-Distance Translocations on Free-Ranging Northern Pacific Rattlesnakes (Crotalus oreganus oreganus). J. Herpetol. 2014, 48, 233–239. [Google Scholar] [CrossRef]
  62. DetectX® Sample Types Validated: Testosterone Enzyme Immunoassay Kit Species Independent, n.d. Available online: www.ArborAssays.com (accessed on 18 April 2020).
  63. Clark, R.W. Diet of the Timber Rattlesnake, Crotalus horridus. J. Herpetol. 2002, 36, 494. [Google Scholar] [CrossRef]
  64. Means, D.B. Diamonds in the Rough: Natural History of the Eastern Diamondback Rattlesnake; Tall Timbers Press: Tallahassee, FL, USA, 2017. [Google Scholar]
  65. Core, R.; Team, R. A Language and Environment for Statistical Computing; R Foundation: Boston, MA, USA, 2019; Available online: https://www.r-project.org/ (accessed on 18 April 2020).
  66. Pawlowsky-Glahn, V.; Egozcue, J.J.; Tolosana-Delgado, R. Modelling and Analysis of Compositional Data; John Wiley & Sons, Ltd.: Chichester, UK, 2015. [Google Scholar] [CrossRef]
  67. Oksanen, J.; Blanchet, F.G.; Friendly, M.; Kindt, R.; Legendre, P.; Mcglinn, D.; Minchin, P.R.; O’hara, R.B.; Simpson, G.L.; Solymos, P.; et al. Package “Vegan”. Community Ecology Package. R Package Version 2.5-6. 2019. Available online: https://cran.r-project.org/package=vegan (accessed on 27 September 2020).
  68. Palarea-Albaladejo, J.; Martín-Fernández, J.A. zCompositions—R package for multivariate imputation of left-censored data under a compositional approach. Chemom. Intell. Lab. Syst. 2015, 143, 85–96. [Google Scholar] [CrossRef]
  69. Elzhov, T.V.; Mullen, K.M.; Spiess, A.-N.; Bolker, B. minpack lm: R Interface to the Levenberg-Marquardt Nonlinear Least-Squares Algorithm Found in MINPACK, Plus Support for Bounds. 2016. Available online: https://cran.r-project.org/package=minpack.lm (accessed on 27 September 2020).
Figure 1. High performance liquid chromatography (HPLC) profiles of C. adamanteus (n = 3) and C. horridus (n = 5) venom samples at 220 nm. The black line reflected into the negative absorbance range shows the profile for the first sample taken as a juvenile, while the red line shows the last sample taken as a large adult. An asterisk indicates significant ontogenetic variation (p < 0.05). Peak height corresponds to absorbance at 220 nm, but heights were scaled to max peak height = 1 and absorbance not shown for simplicity.
Figure 1. High performance liquid chromatography (HPLC) profiles of C. adamanteus (n = 3) and C. horridus (n = 5) venom samples at 220 nm. The black line reflected into the negative absorbance range shows the profile for the first sample taken as a juvenile, while the red line shows the last sample taken as a large adult. An asterisk indicates significant ontogenetic variation (p < 0.05). Peak height corresponds to absorbance at 220 nm, but heights were scaled to max peak height = 1 and absorbance not shown for simplicity.
Toxins 12 00659 g001
Figure 2. Mode and tempo of venom ontogeny. Line graphs show the change in both venom composition (as Aitchison distance from initial samples) and snout-to-vent length over time for all eight snakes studied. Open circles show venom Aitchison distance from early juvenile average, while filled circles show contemporaneous body size (SVL).
Figure 2. Mode and tempo of venom ontogeny. Line graphs show the change in both venom composition (as Aitchison distance from initial samples) and snout-to-vent length over time for all eight snakes studied. Open circles show venom Aitchison distance from early juvenile average, while filled circles show contemporaneous body size (SVL).
Toxins 12 00659 g002
Figure 3. Non-metric multidimensional scaling of venom phenotypes binned into eleven retention time blocks (R1–R11) to allow comparison across all individuals and samples. Polygons are minimum convex shapes enclosing all samples from a single individual snake. Blue dashed arrows point from the earliest sample taken to the last for each individual. HPLC retention time window centroids are also plotted on the chart to clarify the relative position of the samples in venom space. The retention time windows used are as follows for each region designation: R1 (10–20 min), R2 (20–30 min), R3 (30–40 min), R4 (40–50 min), R5 (50–60 min), R6 (60–70 min), R7 (70–80 min), R98 (80–90 min), R9 (90–100 min), R10 (100–110 min), and R11 (110–125 min).
Figure 3. Non-metric multidimensional scaling of venom phenotypes binned into eleven retention time blocks (R1–R11) to allow comparison across all individuals and samples. Polygons are minimum convex shapes enclosing all samples from a single individual snake. Blue dashed arrows point from the earliest sample taken to the last for each individual. HPLC retention time window centroids are also plotted on the chart to clarify the relative position of the samples in venom space. The retention time windows used are as follows for each region designation: R1 (10–20 min), R2 (20–30 min), R3 (30–40 min), R4 (40–50 min), R5 (50–60 min), R6 (60–70 min), R7 (70–80 min), R98 (80–90 min), R9 (90–100 min), R10 (100–110 min), and R11 (110–125 min).
Toxins 12 00659 g003
Figure 4. Peak-by-peak analysis of relative abundance shifts over time for a Type B Crotalus horridus. Panels with best-fit lines indicate peaks that changed significantly in abundance over time. Black lines indicate best fit of a linear model to the data, while blue lines show cases where a non-linear fit was best.
Figure 4. Peak-by-peak analysis of relative abundance shifts over time for a Type B Crotalus horridus. Panels with best-fit lines indicate peaks that changed significantly in abundance over time. Black lines indicate best fit of a linear model to the data, while blue lines show cases where a non-linear fit was best.
Toxins 12 00659 g004
Figure 5. Peak-by-peak analysis of relative abundance shifts over time for a Type A Crotalus horridus. Panels with best fit lines indicate peaks that changed significantly in abundance over time. Black lines indicate best fit of a linear model to the data, while blue lines show cases where a non-linear fit was best.
Figure 5. Peak-by-peak analysis of relative abundance shifts over time for a Type A Crotalus horridus. Panels with best fit lines indicate peaks that changed significantly in abundance over time. Black lines indicate best fit of a linear model to the data, while blue lines show cases where a non-linear fit was best.
Toxins 12 00659 g005
Figure 6. Peak-by-peak analysis of relative abundance shifts over time for an individual Crotalus adamanteus. Panels with best fit lines indicate peaks that changed significantly in abundance over time. Black lines indicate best fit of a linear model to the data, while blue lines show cases where a non-linear fit was best.
Figure 6. Peak-by-peak analysis of relative abundance shifts over time for an individual Crotalus adamanteus. Panels with best fit lines indicate peaks that changed significantly in abundance over time. Black lines indicate best fit of a linear model to the data, while blue lines show cases where a non-linear fit was best.
Toxins 12 00659 g006
Figure 7. Relationship between venom compositional change (Aitchison distance from juvenile phenotype) and natural log-transformed testosterone concentration for six rattlesnakes raised in controlled laboratory conditions, including two Crotalus adamanteus (A,B), two Type A Crotalus horridus (C,D), and two Type B Crotalus horridus (E,F).
Figure 7. Relationship between venom compositional change (Aitchison distance from juvenile phenotype) and natural log-transformed testosterone concentration for six rattlesnakes raised in controlled laboratory conditions, including two Crotalus adamanteus (A,B), two Type A Crotalus horridus (C,D), and two Type B Crotalus horridus (E,F).
Toxins 12 00659 g007
Figure 8. Size-structured dietary ecology in two rattlesnakes. Both Crotalus horridus (top) and Crotalus adamanteus (bottom) show ontogenetic shifts in the contribution of different mammalian prey items to their diet. Diet items from Clark et al. (2002) and Means et al. (2017) show differences among some body size quartiles. Data in cells are proportions of total diet (darker color = greater proportion), and values above each matrix are p-values for chi-squared tests for differences in counts of prey items between adjacent size class quartiles.
Figure 8. Size-structured dietary ecology in two rattlesnakes. Both Crotalus horridus (top) and Crotalus adamanteus (bottom) show ontogenetic shifts in the contribution of different mammalian prey items to their diet. Diet items from Clark et al. (2002) and Means et al. (2017) show differences among some body size quartiles. Data in cells are proportions of total diet (darker color = greater proportion), and values above each matrix are p-values for chi-squared tests for differences in counts of prey items between adjacent size class quartiles.
Toxins 12 00659 g008
Table 1. PERMANOVA results for tests of ontogenetic shifts between each species, with time as a predictor of venom composition.
Table 1. PERMANOVA results for tests of ontogenetic shifts between each species, with time as a predictor of venom composition.
Snake IDSpeciesDfSums of SqsMean SqsFR2p-Value
KW1425C. adamanteus161.1661.163.810.212 × 10−4
KW1726C. adamanteus1263.08263.0820.580.589.9 × 10−5
MM0114C. adamanteus117.7617.764.100.230.0032
KW1576C. horridus1395.08395.0816.980.499.9 × 10−5
KW1594C. horridus1217.70217.7023.430.639.9 × 10−5
MM0027C. horridus1102.5102.569.60.809.9 × 10−5
KW1753C. horridus1210.43210.4314.460.419.9 × 10−5
MM0071C. horridus167.7767.7715.250.459.9 × 10−5
Table 2. Mode and tempo of the ontogenetic shift. Model summaries for both the linear and sigmoid model fits to the relationship between Aitchison distance and time in each snake. Lowest AIC values are in bold. The tempo of the shift is indicated by the predicted number of days to shift at maximum estimated rate of change. Asterisks indicate ∆AIC values >2. The tempo value with a double asterisk is associated with snakes where the simple linear model was a significantly better fit, and the number of days is simply the total number of days over which that individual was sampled.
Table 2. Mode and tempo of the ontogenetic shift. Model summaries for both the linear and sigmoid model fits to the relationship between Aitchison distance and time in each snake. Lowest AIC values are in bold. The tempo of the shift is indicated by the predicted number of days to shift at maximum estimated rate of change. Asterisks indicate ∆AIC values >2. The tempo value with a double asterisk is associated with snakes where the simple linear model was a significantly better fit, and the number of days is simply the total number of days over which that individual was sampled.
Linear Sigmoid
Snake IDSpeciesAICR2p-ValueAICR2p-ValueTempo (Days)
KW1425C. adamanteus22.70.81<0.000119.1 *0.880.0394612
KW1726C. adamanteus22.80.82<0.000124.30.840.0087811
MM0114C. adamanteus27.70.71<0.000127.10.780.0468213
KW1576C. horridus29.50.76<0.000122.5 *0.860.0005989
KW1594C. horridus−0.7 *0.96<0.00011.70.960.0001677 **
MM0027C. horridus17.60.88<0.0001−16.7 *0.980.0000242
KW1753C. horridus51.70.53<0.000144.2 *0.700.0209309
MM0071C. horridus41.30.65<0.000140.90.710.0254719
Table 3. Number of peaks showing different modes of ontogenetic change in each sampled rattlesnake. The percent gradual is the sum of the “Linear-Gradual” and “Sigmoid-Gradual” count divided by the total number of peaks with statistically significant shifts.
Table 3. Number of peaks showing different modes of ontogenetic change in each sampled rattlesnake. The percent gradual is the sum of the “Linear-Gradual” and “Sigmoid-Gradual” count divided by the total number of peaks with statistically significant shifts.
Snake IDSpeciesNo ShiftLinear-GradualSigmoid-GradualSigmoid-DiscreteProportion ShiftingPercent Gradual
KW1425C. adamanteus14131822/36, 61%64%
KW1726C. adamanteus7115723/30, 77%70%
MM0114C. adamanteus1201719/20, 95%11%
KW1576C. horridus1068519/29, 66%74%
KW1594C. horridus557719/24, 79%63%
MM0027C. horridus237616/18, 89%63%
KW1753C. horridus1075416/26, 62%75%
MM0071C. horridus9101011/20, 55%100%
Table 4. Locality, sex, and size data for the C. adamanteus (n = 3) and C. horridus (n = 5) used in this study. CB = Captive born. WC = Wild caught.
Table 4. Locality, sex, and size data for the C. adamanteus (n = 3) and C. horridus (n = 5) used in this study. CB = Captive born. WC = Wild caught.
Snake IDSpeciesSexSVL (cm)LocalitySample Collection Date Range
KW1425C. adamanteusF61–114WC, Lake Co, FL21 May 2013–16 February 2017
KW1726C. adamanteusM40–110.5CB, Little St. George, FL7 October 2013–28 March 2017 *
MM0114 **C. adamanteusF37–95WC, Leon Co, FL19 November 2014–28 March 2017
KW1576C. horridusM34–104.5CB, Baker Co, GA20 August 2013–5 July 2018 *
KW1594C. horridusM42–89.5CB, Baker Co, GA21 September 2013–30 July 2015
MM0027 **C. horridusM64–113WC, Flint River, GA21 May 2013–1 July 2016
KW1753C. horridusM34–100WC, Baker Co, FL12 September 2013–16 February 2017
MM0071C. horridusM41–100.5WC, Columbia Co, FL14 May 2014–16 February 2017
* KW1726 last blood sample date 1 November 2016; KW1576 last blood sample date 28 March 2017. ** Blood samples not assessed.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Schonour, R.B.; Huff, E.M.; Holding, M.L.; Claunch, N.M.; Ellsworth, S.A.; Hogan, M.P.; Wray, K.; McGivern, J.; Margres, M.J.; Colston, T.J.; et al. Gradual and Discrete Ontogenetic Shifts in Rattlesnake Venom Composition and Assessment of Hormonal and Ecological Correlates. Toxins 2020, 12, 659. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins12100659

AMA Style

Schonour RB, Huff EM, Holding ML, Claunch NM, Ellsworth SA, Hogan MP, Wray K, McGivern J, Margres MJ, Colston TJ, et al. Gradual and Discrete Ontogenetic Shifts in Rattlesnake Venom Composition and Assessment of Hormonal and Ecological Correlates. Toxins. 2020; 12(10):659. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins12100659

Chicago/Turabian Style

Schonour, Richard B., Emma M. Huff, Matthew L. Holding, Natalie M. Claunch, Schyler A. Ellsworth, Michael P. Hogan, Kenneth Wray, James McGivern, Mark J. Margres, Timothy J. Colston, and et al. 2020. "Gradual and Discrete Ontogenetic Shifts in Rattlesnake Venom Composition and Assessment of Hormonal and Ecological Correlates" Toxins 12, no. 10: 659. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins12100659

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop