The genetics of venom ontogeny in the eastern diamondback rattlesnake (Crotalus adamanteus)

Animals and tissues

Adult/juvenile pairs of C. adamanteus were collected from five populations previously shown to have different venom compositions (Fig. 1; Margres et al., 2015a; Margres et al., 2015b). Individuals were collected live from the field and temporarily housed at Florida State University; animal information is provided in Table 1. Maturity was assessed on the basis of snout-vent length (SVL), with individuals ≥102.0 cm classified as adults and individuals <102.0 cm classified as juveniles. Waldron et al. (2013) used SVL data to estimate maturity in C. adamanteus and determined that 102.0 cm SVL was the smallest size of reproductively active individuals, and Margres et al. (2015b) found that the ontogenetic shift in toxin gene expression in C. adamanteus corresponded with sexual maturation (i.e., occurred at approximately 1 m SVL). We followed the approach of Rokyta et al. (2012) for preparation of venom glands. Briefly, we stimulated transcription in the glands by means of venom extraction under anesthesia (McCleary & Heard, 2010). Each individual was anesthetized with a propofol injection (10 mg/kg), and venom expulsion was initiated by means of electrostimulation. After allowing four days for transcription to reach maximum levels (Rotenberg, Bamberger & Kochva, 1971), each individual was euthanized by injection of sodium pentobarbital (100 mg/kg). Left and right venom glands were removed and transferred into separate aliquots of RNAlater. Specimens were collected under the following permits: Florida Fish and Wildlife Conservation Commission (FWC) LSSC-13-00004 and LSSC-09-0399 and Florida Department of Environmental Protection permit #03131424. The above procedures were approved by the Florida State University Institutional Animal Care and Use Committee (IACUC) under protocol #0924.

Venom-gland transcriptome sequencing

Snake venom-gland RNA extraction was performed as previously described (Rokyta et al., 2011). Briefly, venom-gland RNA was extracted by first homogenizing the venom gland tissue and submerging in Trizol (Invitrogen), followed by the addition of 20% chloroform and centrifugation in phase lock heavy gel tubes (5Prime) to separate the RNA from DNA and other cellular debris. Isolated RNA was then pelleted with isopropyl alcohol and washed with 75% ethanol. A secondary ethanol precipitation step was used to further purify and concentrate the RNA by using Pellet Paint Co-Precipitant (EMD Millipore), 10% sodium acetate, and 100% ethanol, followed by centrifugation. Purified RNA was then washed with 70% ethanol and quality checked using a Bioanalyzer with an RNA 6000 Pico Kit (Agilent Technologies) following the manufacturer’s instructions.

Prior to RNA-seq library preparation, 1–2 µg of total RNA was used to isolate mRNA using the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs). A fragmentation time of 15.5 min was used to achieve fragment sizes of approximately 370 nucleotides (adapter-ligated). The purified mRNA was immediately used for cDNA library preparation using the NEBNext Ultra RNA Library Prep Kit and Multiplex Oligos for Illumina (New England Biolabs). PCR was performed using the NEBNext High-Fidelity 2X Hot Start PCR Master Mix and 14 cycles of PCR to achieve the desired DNA concentration for sequencing. DNA was purified using Agencourt AMPure XP PCR Purification Beads throughout and at the end of the protocol. The library quality was assessed using a Bioanalyzer with a High Sensitivity DNA Kit (Agilent Technologies) following the manufacturer’s instructions. To determine the amplifiable concentration, KAPA PCR was performed on individual samples by the Florida State University Molecular Cloning Facility. Using the amplifiable concentrations, individual samples with unique indices were pooled to achieve the desired final concentration for sequencing such that each individual snake venom gland was equally represented. The quality of the pooled DNA samples was assessed using a Bioanalyzer with a High Sensitivity DNA Kit (Agilent Technologies), and an additional round of KAPA PCR was performed to confirm amplifiable concentration of the pooled sample prior to sequencing.

All sequencing was performed on an Illumina HiSeq 2000 by the Florida State University College of Medicine Translational Science Laboratory. Samples were multiplexed and sequenced in a rapid run with 150 nucleotide paired-end reads. Left and right glands were processed and sequenced separately for each individual as technical and pseudobiological replicates. The numbers of read pairs per sample are provided in Table 2.

Transcriptome assembly and analysis

Sequencing reads passing the Illumina quality filter were merged on the basis of their 3′ overlaps with PEAR version 0.9.6 (Zhang et al., 2014). Merged reads from both glands for each individual (≤10 million total) were assembled with SeqMan NGen version 12.3.1 with default de novo transcriptome settings, retaining only contigs with ≥200 assembled reads. In addition, 1,000 merged reads for each individual were used as seeds with our in-house assembler Extender (Rokyta et al., 2012), which simply extends the ends of provided seeds on the basis of reads that overlap and match these ends. We only used merged reads as seeds or for extension if all positions had phred quality scores ≥30. We required an exact match of 120 nucleotides for extension. This second de novo assembler was included because of its superior performance at assembling long toxins with high levels of paralogy such as snake-venom metalloproteinases. Note that more commonly used assemblers such as Trinity perform poorly for snake venom transcriptomes (Rokyta et al., 2012; Archer et al., 2014). Toxin-encoding transcripts were identified by means of blastx (v.2.2.30+) searches with a minimum e-value of 10⁻⁴ against the UniProt animal venom proteins and toxins database (http://www.uniprot.org/program/Toxins) downloaded on November 16, 2015. For each assembly, we retained only full-length putative toxins. We combined the results from both assemblies for each of the ten individuals, eliminated duplicates, and then screened for chimeric sequences by aligning all merged reads from both glands against the unique putative toxin coding sequences with bowtie2 version 2.2.7 (Langmead & Salzberg, 2012). Transcripts showing strongly multimodal or extremely uneven coverage distributions were eliminated, and the remaining transcripts were clustered into groups within individuals showing ≤1% nucleotide divergence. A consensus toxin transcriptome across all ten individuals was generated in Geneious version 8.1.8. Putative toxin-transcript coding sequences were aligned by gene family using the ClustalW algorithm (Thompson, Higgins & Gibson, 1994). Transcripts with ≤1.5% nucleotide divergence were combined into clusters by taking the consensus sequence of all cluster members. This clustering of sequences provided an operational definition of paralogs in our final transcriptome. Nontoxin transcripts were derived from a previously published venom-gland transcriptome assembly for C. adamanteus (Margres et al., 2015a). We used RSEM version 1.2.28 (Li & Dewey, 2011) with bowtie version 1.1.2 (Langmead et al., 2009) as the aligner to estimate transcript abundances.

To determine the repeatability of our RNA-seq protocol and transcript-abundance estimates, we aligned left- and right-gland merged reads of each individual separately against the consensus transcriptome and used the estimates of transcripts per million (TPM) from RSEM as our abundance estimates, following the approach of Rokyta, Margres & Calvin (2015). Any estimates of 0.0 TPM were replaced with a value of 1.0 TPM to allow logratio transforms necessary for comparisons across glands for individuals.

To test for the absence of consensus transcripts in any of the ten sequenced transcriptomes, we aligned all of the merged reads from the left and right venom glands from each individual against the consensus transcriptome with BWA MEM (https://sourceforge.net/projects/bio-bwa/), using the -M option. Individual reads that showed more than two mismatches (gaps or nucleotide differences) were removed from the alignments, and we used picard (http://broadinstitute.github.io/picard/) for sorting and indexing. We then used bedtools (Quinlan & Hall, 2010) to calculate the coverage for each site of each of the 59 consensus toxin-encoding transcripts. Transcripts with more than 10% of the coding sequence showing less than 5 × coverage were considered to be absent from the transcriptome. Note that we applied a strict criterion to consider a transcript as present in a transcriptome for this analysis, which was meant to approximate its potential for being assembled de novo. In all other abundance-based analyses, this criterion was not applied. Some transcripts identified as absent for this analysis, therefore, still have abundance estimates for other analyses. To assess the effects of differences in numbers of reads among samples for this analysis, we repeated the analysis using only 9.5 million merged reads for each sample.

To test for evidence of differential expression for toxin genes and to concomitantly account for geographic expression variation, we implemented a pairwise test on adult/juvenile pairs from each population. Because we focused on toxin expression variation, we used available nontoxin expression levels to estimate a null distribution of adult/juvenile expression divergence, then looked for toxins that were outliers relative to this null distribution. Expression levels were estimated as transcripts per million (TPM) using RSEM version 1.2.28 (Li & Dewey, 2011) with bowtie version 1.1.2 (Langmead et al., 2009) as the aligner and centered logratio (clr) transformed. To generate null distributions, we took the absolute values of the difference in transformed adult and juvenile expression levels for each nontoxin transcript and found the value for the 99th percentile. Toxin transcripts showing differences larger than the nontoxin-based 99th percentile were considered outliers for that particular adult/juvenile pair. Because of the known geographic variation for C. adamanteus, we only considered a toxin to show ontogenetic variation if it was an outlier for a majority of the comparisons (≥3 of 5 comparisons) and the change was in the same direction for a majority (≥3 of 5 comparisons) of comparisons (i.e., either upregulated or downregulated in adults). We also tested for differential expression using DESeq version 1.26.0 (Anders & Huber, 2010) and DESeq2 version 1.14.1 (Love, Huber & Anders, 2014), using a 0.1 false-discovery rate (FDR) threshold for both.

Toxicity assays

To generate a pool of adult venoms, we combined 10 mg per individual for 10 adult snakes from the Apalachicola National Forest (ANF). For the juvenile venom pool, we combined 5 mg per individual for 13 juvenile snakes from the ANF population. Five groups of eight mice per venom were housed in cages and observed throughout the experiments. The median lethal doses (LD₅₀s) of snake venoms were determined in BALB/c mice. Venoms were dissolved in 0.85% saline at the highest concentration of venom that was used for the injection (4.89 mg/kg in adult venom and 15.79 mg/kg in juvenile venom). Two-fold serial dilutions using saline were made to obtain four additional concentrations. All solutions during the experiment were stored at 4 °C and warmed to 37 °C just before being injected into mice. Lethal toxicity was determined by injecting 0.2 mL of venom (at various concentrations) into the tail veins of 18–20 g BALB/c mice. The injections were administered using a 1 mL syringe fitted with a 30-gauge, 0.5-in. needle. A saline control was used. The endpoint of lethality of the mice was determined after 48 h. Calculations of final LD₅₀s were determined by the Spearman-Karber method (Spearman-Karber, 1964).

A complete consensus venom-gland transcriptome for Crotalus  adamanteus

We separately assembled and annotated the toxin-encoding genes from the venom-gland transcriptomes (Table 2) from ten individuals (five adults and five juveniles) of C. adamanteus and merged them into a single representative consensus transcriptome for the species. Our final set of transcripts consisted of 59 putative toxin-encoding transcripts. The venom-gland transcriptome of C. adamanteus was previously characterized by means of 454 (Rokyta et al., 2011) and 100-nucleotide paired-end Illumina sequencing (Rokyta et al., 2012; Rokyta, Margres & Calvin, 2015; Margres et al., 2015a). These previous characterizations, however, were all based on RNA from the venom glands of the same single juvenile female from the Apalachicola National Forest (ANF). The most recent assembly (Margres et al., 2015a) identified 44 toxin clusters for C. adamanteus. We have substantially increased this number through a combination of longer reads, which helps the assembly of longer toxins such as snake-venom metalloproteinases (SVMPs), and inclusion of animals from throughout the species’ range. For example, Margres et al. (2015a) identified only seven SVMP paralogs, and we identified 15. The putative toxins and their abundances are provided in Table 3.

To assess the frequency of presence/absence variation of toxin transcripts in venom-gland transcriptomes for C. adamanteus and to determine the potential for missing toxins by sequencing glands from a single individual as the representative of the entire species, we tested for the presence of each of our 59 consensus transcripts in the ten new transcriptomes (Table 4). Seventeen of the 59 transcripts (29%) were absent in at least one of the ten transcriptomes. Only one transcriptome had all 59 (ENP-J), and the most missing was 13 of 59 (CAL-J). The average number missing per transcriptome was 4.5 of 59 (7.6%). The average numbers missing for adults and juveniles were 3.2 and 5.8, respectively, and these numbers were not significantly different (Welch two-sample t test: p = 0.31). The use of a single representative therefore would likely only result in missing a small number of generally low-expression putative toxins, and adults and juveniles should give equivalent characterizations in terms of the toxins identified, even for a species known to show significant geographic (Margres et al., 2015a) and ontogenetic (Margres et al., 2015b; Wray et al., 2015) expression variation. This result, however, may depend on high sequencing coverage; Durban et al. (2013) found more presence/absence differences between a single adult and juvenile venom-gland transcriptome pair for Crotalus simus, but their 454 sequencing approach yielded more than an order of magnitude fewer sequenced base pairs per individual.

The number of merged reads per sample (Table 2) ranged from 9,555,490 (CAL-J) to 20,447,473 (ENP-J). To determine whether this variation among samples contributed to the presence/absence patterns described above, we repeated the analysis using 9.5 million merged reads per sample (Table 4). This change increased the number of transcripts showing presence/absence variation from 17 to 20 of 59, and eliminated the one transcriptome (ENP-J) that had all 59 transcripts. A single transcriptome (LSG-A) lacked only a single transcript. The average number missing per transcriptome increased from 4.5 to 5.4. The average missing for adults was 3.4, the average missing for juveniles was 7.4, and these values were not significantly different (Welch two-sample t test: p = 0.12). Variation in the number of reads per sample therefore did have minor effects on the detectability of some transcripts in some transcriptomes, but these effects did not change our conclusions above.

Because C. adamanteus undergoes a significant ontogenetic change in venom composition (Margres et al., 2015b; Wray et al., 2015), we looked for evidence that certain toxins were absent from adults or juveniles (Table 4). None of the transcripts were entirely unique to either adults or juveniles, suggesting that the ontogenetic change is more quantitative than qualitative. Of the 17 putative toxin transcripts showing presence/absence variation, six (CTL-11, NGF-1, SVMPII-3, SVMPIII-7, and SVSP-12) were only missing from one of the ten transcriptomes. Four of these six were uniquely missing from the CAL-J transcriptome. Three other transcripts (3FTx-1, 3FTx-2, and CTL-12) showed no particular pattern in terms of ontogeny. Of the remaining eight transcripts, six showed a pattern of being present in all adults, but missing from at least two juveniles. These included two type II metalloproteinases (SVMPII-1 and SVMPII-5), two type III metalloproteinases (SVMPIII-1 and SVMPIII-6), and two serine proteinases (SVSP-7 and SVSP-8). The remaining two of 17, CTL-13 and SVMPIII-2, were present in all juveniles, but absent in two adults. Although this analysis only considered qualitative differences in venom composition, we found some evidence for both adult- and juvenile-associated toxins. Adults appear to generally show a higher complexity of proteases, which may be related to the digestive functions of venoms (Mackessy, 2010).

The repeatability of transcript abundance estimates

Because venom-gland transcriptome sequencing necessitates the sacrifice of animals, we will always be limited in the number of true biological replicates that can be performed, particularly for vertebrate species like C. adamanteus that are of conservation concern and typically have relatively small population sizes. To maximize the information gained for each individual snake and to assess the repeatability of our estimates of transcript abundances, we sequenced the left and right venom glands of each individual separately. We found nearly perfect agreement between toxin-transcript abundance estimates across glands for each of our ten individuals (Fig. 2). We measured abundances in transcripts per million (TPM) using RSEM (Li & Dewey, 2011) and used a centered logratio (clr) transform (Aitchison, 1986). This monotonic transform has no effect on rank-order relationships and is equivalent to a log transform for linear relationships. Spearman’s rank correlation coefficients (ρ) ranged from 0.97 through 1.0 (Fig. 2), with six of ten individuals showing exact agreement (i.e., ρ = 1.0). Pearson’s correlation coefficients range from 0.90 through 1.0 (Fig. 2), with four of ten showing essentially perfect linear correspondence (i.e., R = 1.0). Our transcript abundance estimates for toxin transcripts were therefore highly repeatable, and we found no evidence for different expression patterns between left and right glands for individuals. This analysis demonstrated that, for statistical purposes, left and right glands are similar to technical replicates. For characterizing venom composition of a single individual, however, we found no clear benefit to separately sequencing the two glands.

Ontogenetic expression variation

The presence of an ontogenetic venom change in C. adamanteus is well documented. Venom proteomic composition changes through the first postnatal shed (Wray et al., 2015) and appears to have a defined shift at sexual maturity (Margres et al., 2015b). Margres et al. (2015b) used a reversed-phase high performance liquid chromatography (RP-HPLC) approach to quantify ontogenetic and geographic venom expression differentiation for 25 protein peaks and found that the majority of the ontogenetic change could be attributed to five of 25 peaks. They identified proteins encoded by nine loci from a juvenile venom-gland transcriptome for C. adamanteus (Rokyta et al., 2012) in these five RP-HPLC fractions, but four of the five peaks had three or more detectable proteins, so the specific loci involved in the ontogenetic change could not be unambiguously identified. The candidate loci they identified included C-type lectins (CTLs), L-amino acid oxidase (LAAO), type II and type II snake-venom metalloproteinases (SVMPII and SVMPIII), and serine proteinases (SVSPs). This list of loci probably includes loci not involved in the ontogenetic change but that happen to elute in the same peaks as loci that are involved in the change, and the list is possibly incomplete because of the reliance on a single juvenile transcriptome for the identification of proteins in each peak. As shown above, any single individual transcriptome is likely to be missing loci for the species.

We collected adult/juvenile pairs of C. adamanteus from the five populations shown by Margres et al. (2015a) and Margres et al. (2015b) to have different venom compositions (Fig. 1). For each of our five adult/juvenile venom-gland transcriptome pairs, we identified toxin transcripts outside a 99% confidence interval estimate on the basis of nontoxin transcripts (Fig. 3). Toxins outside these confidence intervals in a majority of comparisons (i.e., three of five) with expression changes in the same direction for a majority of comparisons (e.g., upregulated in adults relative to juveniles) were considered to be involved in the ontogenetic change. This outlier approach, requiring a majority consensus, ensured that interpopulation expression variation not related to ontogeny would not result in false positives. Of 59 total toxin transcripts, we identified 13 toxin loci that were upregulated in adults relative to juveniles and four loci that were downregulated in adults relative to juveniles (Table 5). The juvenile-biased toxins included a C-type lectin (CTL), a type II snake-venom metalloproteinase (SVMPII), and two type III snake-venom metallproteinases (SVMPIIIs). The adult-biased toxins included Bradykinin-potentiating and C-type natriuretic peptides (BPP), nerve growth factor (NGF), a phospholipase A₂ (PLA2), two SVMPIIs, four SVMPIIIs, and four snake-venom serine proteinases (SVSPs).

As further confirmation for the results of our outlier analysis, we applied DESeq (Anders & Huber, 2010) and DESeq2 (Love, Huber & Anders, 2014) to compare our five adult transcriptomes to our five juvenile transcriptomes (Table 5). DESeq identified 35 of 4,177 total transcripts as being differentially expressed with a false-discovery rate (FDR) of 0.1. Of these 35, 24 were nontoxins, and 11 were toxins. DESeq2 identified 182 of 4,177 total transcripts as being differentially expressed between adults and juveniles with an FDR of 0.1, including 166 nontoxins and 16 toxins. The three methods showed good agreement (Table 5). Of the 11 toxins detected by DESeq, only one (PLA2-1) was not detected in the outlier analysis. Of the 16 toxins detected by DESeq2, seven were not detected by the outlier analysis. Three of these seven were outliers in at least one of the five adult/juvenile comparisons but did not meet the stringent (i.e., conservative) criterion of being outliers in a majority of comparisons. Of the 17 toxins detected in the outlier analysis, 12 were confirmed by either DESeq or DESeq2 (Table 5).

The 12 loci detected in the outlier analysis and DESeq and/or DESeq2 represented our best conservative estimate of the loci involved in the ontogenetic change in C. adamanteus (highlighted rows in Table 5). Juveniles expressed SVMPII-4, SVMPIII-2, and SVMPIII-4 at higher levels than adults; all of the juvenile-biased toxins were SVMP paralogs. Relative to juveniles, adults upregulated BPP-1, NGF-1, PLA2-2, SVMPII-3, SVMPIII-1, SVMPIII-6, SVSP-7, SVSP-8, and SVSP-9. Therefore, the primary, range-wide ontogenetic change for C. adamanteus involved down regulation of three SVMPs and upregulation of a diverse array of nine toxins from six toxin classes. Snake-venom metalloproteinases induce the local and systemic hemorrhage characteristic of viperid bites and are classified on the basis of their domain structures (Fox & Serrano, 2005; Fox & Serrano, 2010). All SVMPs have a metalloproteinase domain with a zinc-binding motif. Type II SVMPs (SVMPIIs) have an additional disintegrin domain, which may be proteolytically cleaved posttranslationally to produce a free disintegrin. Type III SVMPs (SVMPIIIs) have additional disintegrin-like and cysteine-rich domains. Adults and juveniles of C. adamanteus express similar total amounts of SVMPs. The average for adults was 65,173.7 TPM, and the average for juveniles was 75,654.6 TPM. These averages were not statistically different (Welch two-sample t test: p = 0.68). The ontogenetic change therefore resulted in a maintenance of overall SVMP expression with a shift in expression from juvenile to adult paralogs. The adult-biased toxins include a diverse set of toxin types. Snake-venom serine proteinases (SVSPs) interfere with blood coagulation and hemostasis and belong to the trypsin family of serine proteases (Serrano & Maroun, 2005; Phillips, Swenson & Francis S. Markland, 2010). Bradykinin-potentiating and C-type natriuretic peptides are presumed to cause a reduction in blood pressure (Pahari, Mackessy & Kini, 2007). Despite being present in diverse snake venoms, nerve growth factor (NGF) has an unknown role as a component of venoms (Kostiza & Meier, 1996; Lavin et al., 2010) but was also more highly expressed in the venom-gland transcriptome of an adult C. simus than a juvenile (Durban et al., 2013). Phospholipases A₂ are among the most functionally diverse classes of snake-venom toxins and have pharmocological effects that include neurotoxicity (presynaptic or postsynaptic), myotoxicity, and cardiotoxicity. Anticoagulant and hemolytic effects due to PLA2s are also known (Lynch, 2007; Doley, Zhou & Kini, 2010). Crotalus adamanteus expresses only one PLA2 transcript at a high level, but a second (PLA2-2) is expressed at low levels (Table 3; Dowell et al., 2016). The low-expression paralog was detected as a component of the ontogenetic change. Although all of the adult-biased toxins were also expressed in juveniles, their upregulation in adults should result in a more functionally complex venom. Adults and juveniles share the same genomes and are therefore confined to the same sets of available toxin genes. The presence of toxins differentially expressed in adults and juveniles may ultimately allow us to assess the relative strength of selection acting on venoms of adults and juveniles by comparing evolutionary patterns in the sequences of juvenile-biased and adult-biased toxin genes.

Durban et al. (2013) conducted a proteomic and transcriptomic characterization of the ontogenetic change in C. simus and showed that the proteome-based venom profiles showed less agreement than the transcriptome-based profiles. In particular, they claimed that the venom-gland transcriptomes were largely indistinguishable between adults and juveniles, despite major divergence in the venom proteomes. They therefore attributed at least some portion of the ontogenetic change to posttranscriptional regulation mediated by microRNAs. They used a single adult and juvenile for transcriptomics and pooled adult and juvenile venoms for proteomics, which leaves no real possiblity of a true statistical comparison between adults and juveniles or between transcriptomes and proteomes. Conclusions as significant as the implication of specific posttranscriptional regulatory mechanisms are best not made on the basis of samples sizes of one. We showed that the ontogenetic change in C. adamanteus has a simple, statistically significant transcriptional basis, which was consistent with the proteomic patterns underlying the ontogenetic change in C. adamanteus described by Margres et al. (2015b) and the functional characterizations described by Margres et al. (2016a).

Ontogenetic toxicity variation

To determine whether the significantly different venom profiles for adults and juveniles of C. adamanteus detected by Margres et al. (2015b), which are underlain by the transcriptomic differences described above, resulted in different effects in model prey, we measured median lethal doses (LD₅₀s) of pooled adult and juvenile venoms in mice. We used venoms from adults and juveniles from the Apalachicola National Forest (ANF) population, because the venom-gland transcriptomic comparison for this population (Fig. 3) showed all of the differences detected in our outlier analysis. We estimated an adult LD₅₀ of 3.46 mg/kg and a juvenile LD₅₀ of 2.79 mg/kg, indicating that the juvenile venom was more potent in mice. The differential venom-transcript expression described above therefore translates to a toxicity difference in the venoms. Although this toxicity difference was detected in a model prey species, the difference clearly indicated large phenotypic differences in the functions and actions of adult and juvenile venoms. The difference therefore implies a difference in the effects of these venom in natural prey species, although the magnitude and direction of the difference is likely to be different in natural prey. In other species with known ontogenetic venom compositional changes, juveniles tend to show more toxic venoms (Mackessy, 1988; Saldarriaga et al., 2003; Mackessy et al., 2006). This result may depend, however, on the choice of prey species used in the toxicity assays. Juveniles of Bothrops jararaca have more toxic venoms than adults in chicks, but adults are more potent against mice (Zelanis et al., 2010), suggesting that the difference between adult and juvenile venoms may simply maintain higher toxicity to the different preferred prey of each of the two age classes.