Next Article in Journal
Insights into the Cardiotoxic Effects of Veratrum Lobelianum Alkaloids: Pilot Study
Previous Article in Journal
Mycotoxins of Concern in Children and Infant Cereal Food at European Level: Incidence and Bioaccessibility
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Divergent Specialization of Simple Venom Gene Profiles among Rear-Fanged Snake Genera (Helicops and Leptodeira, Dipsadinae, Colubridae)

by
Peter A. Cerda
1,2,*,
Jenna M. Crowe-Riddell
1,2,3,
Deise J. P. Gonçalves
1,
Drew A. Larson
1,4,
Thomas F. Duda, Jr.
1,2 and
Alison R. Davis Rabosky
1,2
1
Ecology and Evolutionary Biology, University of Michigan, Ann Arbor, MI 48109, USA
2
Museum of Zoology, University of Michigan, Ann Arbor, MI 48108, USA
3
School of Agriculture, Biomedicine and Environment, La Trobe University, Melbourne, VIC 3086, Australia
4
Department of Biology, Indiana University, Bloomington, IN 47405, USA
*
Author to whom correspondence should be addressed.
Toxins 2022, 14(7), 489; https://doi.org/10.3390/toxins14070489
Submission received: 3 June 2022 / Revised: 11 July 2022 / Accepted: 12 July 2022 / Published: 15 July 2022
(This article belongs to the Section Animal Venoms)

Abstract

:
Many venomous animals express toxins that show extraordinary levels of variation both within and among species. In snakes, most studies of venom variation focus on front-fanged species in the families Viperidae and Elapidae, even though rear-fanged snakes in other families vary along the same ecological axes important to venom evolution. Here we characterized venom gland transcriptomes from 19 snakes across two dipsadine rear-fanged genera (Leptodeira and Helicops, Colubridae) and two front-fanged genera (Bothrops, Viperidae; Micrurus, Elapidae). We compared patterns of composition, variation, and diversity in venom transcripts within and among all four genera. Venom gland transcriptomes of rear-fanged Helicops and Leptodeira and front-fanged Micrurus are each dominated by expression of single toxin families (C-type lectins, snake venom metalloproteinase, and phospholipase A2, respectively), unlike highly diverse front-fanged Bothrops venoms. In addition, expression patterns of congeners are much more similar to each other than they are to species from other genera. These results illustrate the repeatability of simple venom profiles in rear-fanged snakes and the potential for relatively constrained venom composition within genera.
Key Contribution: This paper provides novel comparative data on the diversity of venom gene expression in both well-known front-fanged and understudied rear-fanged snakes from the most biodiverse region on Earth: the Amazon.

1. Introduction

Many animal groups use venoms that are comprised of toxic compounds to subdue prey and defend against predators [1]. Venom composition tends to differ considerably among closely related taxa within these groups, possibly due to differences in the type and diversity of prey species [2,3,4,5,6,7,8]. While venoms of a variety of taxa have been characterized [1,4,9,10], the venom composition of members of related but understudied groups is important for reconstructing the evolution of these complex phenotypes across taxa.
Much of our knowledge of snake venoms is based on information from front-fanged species of the families Viperidae (e.g., rattlesnakes) and Elapidae (e.g., coral snakes) [6,11,12,13]. These snakes inject venom through hypodermic needle-like fangs that are either hinged (viperids) or fixed (elapids) at the front of the mouth (Figure 1B) and utilize a high-pressure venom delivery system (Figure 1A) [14]. While most venom components can be found throughout venomous snake species, the venom composition of these two families tends to be drastically different [14]. Viperid venoms are typically hemorrhagic or cytotoxic and largely contain snake venom metalloproteinases (SVMPs), snake venom serine proteases (SVSPs), and phospholipase A2s (PLA2s) [15]. On the other hand, venoms of elapids are usually neurotoxic and primarily contain either PLA2s, three-finger toxins (3FTxs), or a combination of these two types [15]. While these toxin families dominate the venom profiles of front-fanged snakes, other toxins, such as C-type lectins (CTLs), are also found at varying levels [15]. Front-fanged species can differ in venom complexity [16], ranging from “simple” venoms comprised mostly of one toxin family to “complex” venoms comprised of many toxin families [6] that may be associated with inter- and intraspecific differences in predator–prey interactions [2,5,6,17].
Less of our knowledge on snake venoms comes from rear-fanged snakes of the family Colubridae (sensu Pyron et al. [18]), which makes up approximately half of all snake species [19]. Venomous colubrids produce toxins in the oral Duvernoy’s gland (hereafter referred to as the “venom gland”; Figure 1A), which are delivered through grooved or ungrooved fangs located at the back of the mouth (Figure 1B) via a low-pressure venom delivery system [19,20]. Although the family contains approximately 700 venom-producing paraphyletic rear-fanged species [14], venoms have only been investigated in a few of them (Table 1). Emerging trends seem to suggest that Colubrid venoms tend to have relatively simple compositions and that the subfamilies appear to follow a compositional dichotomy [21]. Venoms of the dipsadine subfamily tend to be dominated by SVMPs, such as some viperids, while venoms of the colubrine subfamily largely contain 3FTxs, such as some elapids (Table 1 and references within). Variation in the venoms of colubrid snakes has not been well studied, but there is evidence of ontogenetic shifts in venom composition of Boigia irregularis [22], as well as some geographic variation in the venom composition of Tantilla nigriceps [23]. Nonetheless, the venom compositions of rear-fanged snakes are still largely unknown, limiting our ability to both describe general patterns in the diversity of colubrid venoms as well as accurately model the evolutionary dynamics of snake venoms more broadly.
To increase our understanding of the venom composition of rear-fanged snakes, we characterized venom gland transcriptomes of members of the dipsadine genera Helicops and Leptodeira. While venoms of members of these genera have not previously been examined through RNA sequencing approaches, they have been subject to functional studies, and so some inferences about their composition are available [32,33,34,35,36,37]. For example, Leptodeira annulata and its subspecies Leptodeira annulata pulchriceps have venoms that differ in terms of serine protease activity, but both appear to be rich in SVMPs and PLA2s as they exhibit proteolytic, hemorrhagic, and neurotoxic activities [33,34,35]. The venom of Helicops angulatus exhibits neurotoxic but not hemorrhagic activity and contains a cysteine-rich secretory protein, termed helicopsin, that causes respiratory paralysis in mice [36,37]. The absence of hemorrhagic activity implies that venoms of H. angulatus do not contain SVMPs, as these metalloproteases tend to induce hemorrhaging [38]. Hence, not all dipsadine snakes apparently possess a viperid-like venom, and some members of this subfamily may instead contain members of other toxin types.
We provide a first characterization of the venom gland transcriptomes of members of the dipsadine genera, Helicops and Leptodeira, and compare them with venom gland transcriptomes of front-fanged snakes of the genera Bothrops (Viperidae) and Micrurus (Elapidae). To do this, we sequenced venom gland transcriptomes from 14 species of these four genera, including multiple individuals of four species. We determined (i) which major toxin families are expressed, (ii) whether these venoms are simple or complex, and (iii) inter- and intraspecific levels of variation in venom composition among members of these genera. Although the venom profiles of Bothrops and Micrurus species have been characterized previously, we generated new sequence data for these genera to enable effective comparisons of transcriptomes generated via the same library preparation methods, sequencing approaches, and bioinformatic procedures.

2. Results

2.1. Venom Gland Gene Family Recovery

We produced and analyzed the venom gland transcriptomes of 19 individuals of six rear-fanged snake species (Helicops n = 5 and Leptodeira n = 4; Table 2) and eight front-fanged snake species (Bothrops n = 3 and Micrurus n = 8; Table 2). The number of paired reads recovered per sample ranged from 12,936,858 to 24,938,732 (Table 3). We found that the total toxin sequence count was associated with sequencing platforms used: transcriptomes that were sequenced on a HiSeq 4000 produced a higher toxin transcript count than those on a NovaSeq 6000, including transcriptomes of conspecifics or congeners. The number of unique toxin transcripts recovered from the venom gland of Helicops species ranged from 47 to 91 (Figure 1C). The most frequently identified toxin families in Helicops were snake venom metalloproteinase (SVMPs) and C-type lectins (CTLs). Several copies of cysteine-rich secretory proteins (CRiSPs) were recovered from two individuals of H. angulatus and one individual of H. leopardinus (Figure 1C). The number of unique toxins recovered from Leptodeira venom gland transcriptomes ranged from 29 to 99 (Figure 1D). Similar to the rear-fanged Helicops, the most common toxin families in Leptodeira were SVMPs, CTLs, and PLA2s (Figure 1D). The number of unique toxins recovered from Micrurus species ranged from 37 to 141. These largely included SVMPs, PLA2s, and CTLs transcripts (Figure 1E). The number of unique toxin transcripts recovered from Bothrops species ranged from 89 to 131, and the most common toxin families recovered were SVMPs, snake venom serine proteinases (SVSPs), CTLs, and PLA2s (Figure 1F).

2.2. Venom Gland Transcriptome Expression

In total, we found that total toxin expression encompassed between 17% and 91% of the total expression in the venom gland transcriptomes studied (Table 3). We found that CTLs were the most highly expressed toxin family of Helicops and comprised between 63 and 99% of all toxins expressed (Figure 2). SVMPs were the second most abundant toxin family expressed (0.6–36%) in all but one individual of H. angulatus, for which CRiSPs toxins made up 24% of the venom expression profile (Figure 2). We found that SVMPs also comprised a considerable portion of the venom gland expression profile of rear-fanged Leptodeira species and encompassed 83–98% of the expressed toxins (Figure 2). The second most highly expressed toxins in Leptodeira were CTLs which ranged from 0.7 to 7% of expressed toxins (Figure 2). In front-fanged Micrurus species, we found that venom gland expression was dominated by one or two toxin families. All Micrurus species expressed PLA2s at a high level (49–99%; Figure 2). Kunitz-type serine protease inhibitors were the second most highly expressed toxin in both individuals of M. lemniscatus that we examined (22–26%; Figure 2), whereas 3FTxs were the second most highly expressed toxins of all other Micrurus species (13–21%; Figure 2). In front-fanged Bothrops atrox and Bothrops brazili, SVMPs had the highest expression, comprising 48 and 50% of the venom profile, respectively. (Figure 2). Bradykinin-potentiating peptides were the second most abundant toxin family found in B. atrox (28%), while CTLs were the second most abundant toxin found in B. brazili (20%; Figure 2). SVMPs and CTLs were expressed at similar levels (~30%) in Bothrops bilineatus (Figure 2).

2.3. Complexity and Variation

We calculated Shannon diversity indices to compare levels of venom complexity across individuals. These values ranged from a low of 0.662 for Helicops leopardinus to a high of 3.251 for Bothrops brazili (Figure 2). Overall, the four genera ranked from lowest to highest levels of venom complexity as follows: Helicops (1.289), Micrurus (1.602), Leptodeira (1.954), and Bothrops (3.014). We further estimated beta diversity statistics to determine levels of intra- and interspecific variation in venom composition. For the four species for which multiple individuals were examined, L. annulata and M. lemniscatus exhibited the lowest values (0.097 in both cases), while H. angulatus (0.239) and M. obscurus (0.497) had higher values. Within-genera comparisons showed that Leptodeira (0.097) and Bothrops (0.355) exhibited the lowest and highest values for genera, respectively, whereas Helicops (0.190) and Micrurus (0.180) had intermediate values. While intraspecific variation was found among members of these genera, we saw a greater similarity in venom gland transcriptome composition among individuals within genera than among genera (Figure 2).

3. Discussion

We used a transcriptomic approach to characterize and compare patterns of toxin variation in the venom gland transcriptomes of members of two rear-fanged (Helicops and Leptodeira) and two front-fanged (Bothrops and Micrurus) snake genera. While snake venom metalloproteinases (SVMPs) transcripts dominate the venom profile of the Leptodeira species, C-type lectins (CTLs) are the most highly expressed toxin family of Helicops species (Figure 2). Venom profiles of the front-fanged species are similar to those reported from these taxa previously [39,40,41]. While we were able to recover toxin transcripts from numerous toxin families across all individuals, we found that toxin expression within genera is typically dominated by only a single toxin class: CTLs in Helicops; phospholipase A2s (PLA2s) in Micrurus; and SVMPs in both Leptodeira and Bothrops. Shannon diversity indices suggest that while complexity may vary among individuals, differences may be due to the relative contribution of each underlying transcript rather than differences in toxin family abundance. We also observed low levels of variation among individuals of the same genera but high levels of variation across genera.
The dipsadine rear-fanged species that we examined mostly exhibited low levels of venom complexity, given that their venom gland expression profiles were largely composed of transcripts of only single toxin classes (Figure 2). In line with previous results, front-fanged snakes tended to exhibit higher levels of venom complexity. Increased venom complexity has been correlated with large dietary breadth in both venomous cone snails and North American pit vipers [4,6]. The venoms of several rear-fanged snakes are known to have a “simpler” venom than their front-fanged relatives [21,26,27,42,43]. The lack of complexity observed for rear-fanged snakes may be due to the highly specialized diets that many colubrid snakes exhibit [14]. Broadly, Leptodeira appear to specialize in frogs while Helicops specialize in fish, though a formal analysis of dietary specialization in these two genera has not been performed [44,45]. However, the toxicological diversity of specific toxin families from the species described here has not been determined. The investigation of toxin function may reveal physiological targets and functions specific to prey types, which have been found in neurotoxins described from several rear-fanged snakes [22,26].
Several members of the colubrid subfamily Dipsadinae have been shown to use a hemorrhagic venom that is largely comprised of SVMPs [21,42]. The Leptodeira species we examined exhibited this type of venom profile. However, Helicops species had venom gland expression profiles that are dominated by transcripts of a single toxin family, CTLs. No previous studies have shown snakes that produce venom dominated by CTLs. While functional attributes of CTLs of rear-fanged snakes are not known, CTLs of front-fanged snakes are multifunctional heterodimers that affect hemostasis by acting as anticoagulants, which can cause hemorrhaging, or as procoagulants, which can cause blood clotting [46]. In addition, CTLs of front-fanged snakes have been shown to evolve rapidly [47]. Recently, Xie et al. [48] found extensive duplication of novel dimeric CTLs genes unique to H. leopardinus. The novel CTLs found in H. leopardinus were shown to be under positive selection, but the distribution of these CTLs across the genus Helicops is currently unknown [48]. A previous study of H. angulatus found its venom lacked hemorrhagic activity, which is typically associated with SVMPs and CTLs [36]. Instead, its venom was shown to exhibit neurotoxic activity, likely due to the presence of a previously uncharacterized cysteine-rich secretory protein (CRiSP) named helicopsin [36]. Of the H. angulatus individuals examined here, only two possess CRiSP transcripts, and only one expressed CRiSP at a considerable level (Figure 2). It is unknown if the CRiSP transcripts found in our Peruvian individuals are similar to helicopsin that was previously isolated from an individual of this species from Venezuela [36].
Our sampling included multiple individuals of some species and multiple species of four genera to evaluate patterns of intra- and interspecific variation of venom profiles as inferred from transcriptome data. The levels of variation within species differed considerably. Intraspecific variation of venom profiles of L. annulata and M. lemniscatus were low, while those of H. angulatus and M. obscurus were quite distinct. For example, individuals of H. angulatus differed in terms of the relative abundance of CRiSPs, SVMPs, and CTLs, despite being from the same locality. Further, the described absence of hemorrhagic activity of venoms of H. angulatus from Venezuela [37] suggests that this species exhibits geographic variation in venom. While venom gland profiles generated here generally match venom profiles that were described previously in other studies, geographic variation in venom composition may occur among some of the species surveyed here as well [12,40,41,49]. For example, while pooled venom of B. brazili from Pará, Brazil, contained 33% PLA2s, 27% SVMPs, and 14% SVSPs [49], SVMPs and CTLs were the predominant toxin components that are represented in the venom gland transcriptome of an individual from the Madre de Dios region of Peru. Geographic variation is rather common in predatory venomous species as populations presumably adapt to local prey assemblages [50,51]. However, caution should be taken when comparing venom gland transcriptomes and venom proteomes as there are cases in which there is a lack of correspondence between the abundance of expressed transcripts and translated proteins [40]. Further, the proportion of toxin expression in the whole transcriptome varied widely among the specimens examined (17.15–91.81%; Table 3). The variation in venom expression can possibly be attributed to differences in toxin expression over time [52]. For example, toxin expression is higher after feeding events to replenish the used venom proteins [52].
The two individuals of M. obscurus that we examined differed considerably in venom complexity, with the venom profile of one individual being nearly completely dominated by PLA2s, while that of the other individual was more complex and contained wapirins, 3FTxs, and Kunitz-type serine protease inhibitors (Figure 2; note that wapirins are included in the “Other” category in the figure). These two individuals differed in age, as one was a juvenile (M. obscurus 0665), while the other was an adult (M. obscurus 1054). Thus, the difference observed in M. obscurus may be due to an ontogenic shift from a more “complex” venom that is expressed by juveniles to a “simple” venom that is expressed by adults. Ontogenetic shifts in venom composition have been observed in many snake species, and this change is potentially due to differences in the diet as individuals age [17,22,53,54]. However, more intensive sampling is needed to determine if this is indeed the case in M. obscurus.
The species studied here displayed varying levels of interspecific variation in venom gland transcriptomes, which is a common pattern observed across venomous taxa [1,4,6,55]. Barua and Mikheyev [56] found that while many combinations of venom components were possible, different species tend to show similar venom profiles despite phylogenetic relatedness, albeit with different proportions of respective venom compositions. While interspecific variation was observed among species examined here, we do note that most of our taxa had venom gland profiles similar to those expected given their respective families [21,56], except for species of Helicops. It is not clear how or why the Helicops species examined here arrived at a potentially novel venom phenotype. The use of CTLs by Helicops species may be a more efficient strategy to capture their aquatic prey [44,45]. Proteomic and functional studies should be performed on Helicops venoms to determine the abundance of CTL proteins in these venoms and how they might be used in prey capture. Future exploration of colubrid venoms will help us further our understanding of how convergent and novel venom phenotypes have evolved across all venomous snakes.

4. Materials and Methods

4.1. Sampling

We collected individuals of Bothrops, Helicops, Leptodeira, and Micrurus species from multiple localities in Nicaragua and Peru between 2016 and 2018 (Table 2). Within ten minutes of euthanasia, we excised venom gland tissues and stored them in RNALater (Invitrogen, Carlsbad, CA, USA). While at remote locations, we stored samples at room temperature for up to three weeks during active collection expeditions. Then we stored the samples at −20 °C prior to export to the University of Michigan and subsequently at −80 °C until RNA extraction. We deposited whole voucher specimens in the Museo de Historia Natural, Universidad Nacional Mayor de San Marcos (MUSM) in Lima, Peru, and the University of Michigan Museum of Zoology (UMMZ; Table 2).

4.2. Extraction, Library Preparation, and Sequencing

We extracted the total RNA from venom glands using the PureLink RNA mini kit (Life Technologies, Carlsbad, CA, USA), following the recommended protocol for animal tissue. The total RNA was stored at −80 °C until submission to the University of Michigan’s Advanced Genomics Core, where the RNA was quantified with a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA), and the size was visualized with an Aligent TapeStation (Santa Clara, CA, USA). Poly-A tail selected libraries were constructed with Illumina TruSeq RNASeq (San Diego, CA, USA) and NEBNext Ultra II RNA (Ipswich, MA, USA) library kits, and 150 bp paired-end sequencing was conducted on Illumina HiSeq4000 or Illumina NovaSeq6000 machines (Table 3).

4.3. Bioinformatics

We assessed the raw read quality for each individual using FastQC v0.11.6 [57] and used Trimmomatic v0.38 to trim adapters and remove low-quality reads [58] for downstream phylogenetic and transcriptomic analysis. We assembled the transcriptome of each individual using two methods, Trinity v2.6.6 [59] and Extender [60], a seed-and-extend assembler that has been shown to recover a high level of isoform diversity in snake venom families [61]. The extender assemblies were generated with combined trimmed forward and reverse paired reads that were merged with PEAR v0.9.6 [62].
We merged the assembled transcriptomes produced by Trinity and Extender to construct a single transcriptome for each individual. For each transcriptome, we followed a previously published protocol to filter out low-quality transcripts and chimeras [63,64] using transRate v1.0.3 [65] and a BLAST-based method [66], respectively. We used Corset v1.07 [67] to remove duplicate transcripts and select a single representative transcript (the longest) for each putative gene, using SALMON v0.11.2 [68] as our aligner option. To find open reading frames, we used transDecoder [69] and BLAST [70]; CD-HIT was used to reduce redundancy at 95% similarity [71]. We estimated RNA abundance (i.e., transcripts per million [TPM]) with the align_and_estimate_abundance.pl script in Trinity [59], which calls on RSEM v1.2.28 [72] and Bowtie2 v2.3.4.1 [73].
We created a custom database of non-toxin and toxin nucleotide sequences by downloading venom gland transcriptomes of Crotalus adamanteus [60], Crotalus horridus [74], Micrurus fulvius [75], Boiga irregularis [21], Hypsiglena sp. [21], and Spilotes sulphureus [26], the annotated genome of Ophiophagus hannah [76], and sequences of snake venom matrix metalloproteinase [77] from NCBI GenBank. A combination of these species has been used previously to identify putative rear-fanged snake venom components and represent a diversity of venomous snakes [26]. We used BLASTn [70] to identify known toxins in our transcriptomes using this custom database. We wrote a custom script in R [78] to annotate our nucleotide sequences, sort non-toxin and toxin transcripts by identifying annotations that matched known toxins from the custom database using the ‘grep’ function, and associate sequence identity with transcript abundance estimates. To determine the composition of toxin transcripts in the venom gland transcriptomes, we wrote a custom R script that counted the number of venom gene contigs and total TPM for each snake toxin gene family across individuals. Toxin sequences with a TPM of 0 were removed. We divided the total TPM of each toxin family over the total TPM of all toxin genes to give a proportion of each toxin family present in the venom gland transcriptome.

4.4. Assembly of Mitochondrial Sequences and Phylogenetic Tree Estimation

We assembled mitochondrial sequences for our samples using HybPiper v1.3.1 [79] and a target file consisting of sequences from 16 fully annotated, publicly available mitochondrial genomes downloaded from GenBank [80,81,82,83,84,85] (Table S1). To prepare the target file, all sequences annotated as rRNA, tRNA, and protein-coding genes were extracted from the mitogenome using Geneious Prime 2020.2.3 (https://www.geneious.com (accessed on 5 January 2022)). We then manually curated sequence label formatting and combined all sequences into a single HybPiper target file. We processed the raw sequence reads for all of the transcrioptomes with trimmomatic [58] and used the options “ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:2:TRUE SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:36” to trim adapters and reads with a PHRED score of less than 20. We spot-checked the trimmed reads with FASTQC [57]. We combined reads that became unpaired during trimming into a single file per sample and assembled target mitochondrial genes with HybPiper using forward, reverse, and unpaired reads and default settings. The HybPiper pipeline calls on Exonerate [86], BLAST+ [70], Biopython [87], BWA [88], SAMtools [89], GNU Parallel [90], and SPAdes [91].
The assemblies for protein coding, rRNA, and informative tRNA sequences were aligned with MAFFT v7.271 [92] and the options “–maxiterate 1000”, “–ep 0.123”, and “–genafpair”. Columns with >70% gaps were removed, and alignments were concatenated into a supermatrix with the pxclsq and pxcat commands in phyx, respectively [93]. We estimated a phylogenetic tree using IQ-TREE v 2.1.3 [94] and the options “–m TEST” and “–mset raxml” to determine the best-fitting model, which had 19 partitions, variously assigned GTR + F, GTR + F + I, GTR + F + G4, and GTR + F + I + G4 models. The maximum likelihood tree was visualized with Figtree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/ (accessed on 10 January 2022)).

4.5. Complexity and Variation

To estimate the levels of venom complexity, we calculated Shannon indices (H′) [95] based on each unique toxin and their respective TPM expressed in the transcriptomes of the individuals examined. For cases in which multiple individuals of a species were examined, we averaged the H′ values from these individuals to estimate the venom complexity of the species. We also averaged the values across species to estimate the relative levels of venom complexity of the genera. We quantified the extent that samples differ in levels of intra- and interspecific variation of venom composition with calculations of pairwise proportional dissimilarity (PPD) values (i.e., Brays–Curtis distances, [96]) based on proportions of toxin families recovered. We used the PPD values to estimate and compare levels of intraspecific variation in venom composition for species from which multiple individuals were examined (i.e., Helicops angulatus, Leptodeira annulata, Micrurus lemniscatus, and Micrurus obscurus). We further averaged PPD values that were calculated among species of genera to evaluate levels of interspecific variation in venom composition.

4.6. diceCT and Segmentation

We used diffusible iodine contrast-enhanced computed tomography (diceCT) to scan a representative specimen from each genus using 1.25% Lugol’s iodine solution and a Nikon Metrology XTH 225ST microCTscanner (Xtect, Tring, UK) at the UMMZ, following protocols outlined in Callahan et al. [97]. We segmented the maxillary bone and venom gland using the ‘draw’ and ‘thresholding’ tools in Volume Graphics Studio Max version 3.2 (Volume Graphics, Heidelberg, Germany).

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/toxins14070489/s1, Table S1: GenBank accession numbers for mitochondrial sequences by species.

Author Contributions

Conceptualization, P.A.C., T.F.D.J. and A.R.D.R.; Methodology, P.A.C., J.M.C.-R., D.J.P.G. and D.A.L.; Formal Analysis, P.A.C., D.J.P.G. and D.A.L.; Investigation, P.A.C., D.J.P.G. and D.A.L.; Resources, P.A.C. and A.R.D.R.; Data Curation, P.A.C., D.J.P.G. and D.A.L.; Writing—Original Draft Preparation, P.A.C.; Writing—Review & Editing, P.A.C., J.M.C.-R., D.J.P.G., D.A.L., T.F.D.J. and A.R.D.R.; Visualization, P.A.C., J.M.C.-R. and A.R.D.R.; Supervision, T.F.D.J. and A.R.D.R.; Project Administration, A.R.D.R.; Funding Acquisition, P.A.C. and A.R.D.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the University of Michigan to A.D.R. through startup funds and to P.A.C. through the C.F. Walker Scholarship and Graduate Research Awards. This research was also supported by a Theodore Roosevelt Memorial Fund award from the American Museum of Natural History to P.A.C. Field research was partially supported through funding from the David and Lucile Packard Foundation to Dan Rabosky and the University of Michigan to Iris Holmes.

Institutional Review Board Statement

The animal study protocol was approved by the Institutional Animal Care and Use Committee of University of Michigan (protocols PRO00006234, approved 4/22/2015, and PRO00008306, approved 11/08/2018).

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw sequence data for each individual are available on NCBI SRA under BioProject PRJNA843733 and BioSample accession SAMN28768753-SAMN28768771.

Acknowledgments

We thank Dan Rabosky, Rudi von May, Joanna Larson, Jose Martínez-Fonseca, and Iris Holmes for assistance in organizing and running field collection in Peru and Nicaragua, along with the many people who helped with collection in the field. We thank Conservación Amazónica (ACCA) and Project Amazonas for support at field stations. We also thank SERFOR (Servicio Nacional Forestal y de Fauna Silvestre) in Peru, MARENA (Ministerio del Ambiente y los Recursos Naturales) in Nicaragua, and the U.S. Fish and Wildlife Service for collection, export, and import permits. We thank Matt Holding and Darin Rokyta for access to the Extender program, and Ramon Nagesan, Sean Callahan, and Greg Schneider for assistance with CT scanning. We also thank the Davis Rabosky and Duda lab members and anonymous reviewers for helpful comments that strengthened this manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Casewell, N.R.; Wüster, W.; Vonk, F.J.; Harrison, R.A.; Fry, B.G. Complex Cocktails: The Evolutionary Novelty of Venoms. Trends Ecol. Evol. 2013, 28, 219–229. [Google Scholar] [CrossRef]
  2. Daltry, J.C.; Wüster, W.; Thorpe, R.S. Diet and Snake Venom Evolution. Nature 1996, 379, 537–542. [Google Scholar] [CrossRef]
  3. Sanz, L.; Gibbs, H.L.; Mackessy, S.P.; Calvete, J.J. Venom Proteomes of Closely Related Sistrurus Rattlesnakes with Divergent Diets. J. Proteome Res. 2006, 5, 2098–2112. [Google Scholar] [CrossRef]
  4. Phuong, M.A.; Mahardika, G.N.; Alfaro, M.E. Dietary Breadth Is Positively Correlated with Venom Complexity in Cone Snails. BMC Genom. 2016, 17, 1–15. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Lyons, K.; Dugon, M.M.; Healy, K. Diet Breadth Mediates the Prey Specificity of Venom Potency in Snakes. Toxins 2020, 12, 74. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Holding, M.L.; Strickland, J.L.; Rautsaw, R.M.; Hofmann, E.P.; Mason, A.J.; Hogan, M.P.; Nystrom, G.S.; Ellsworth, S.A.; Colston, T.J.; Borja, M.; et al. Phylogenetically Diverse Diets Favor More Complex Venoms in North American Pitvipers. Proc. Natl. Acad. Sci. USA 2021, 118, 2015579118. [Google Scholar] [CrossRef]
  7. 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 Biol. Sci. 2009, 276, 2443–2449. [Google Scholar] [CrossRef] [Green Version]
  8. Remigio, E.A.; Duda, T.F. Evolution of Ecological Specialization and Venom of a Predatory Marine Gastropod. Mol. Ecol. 2008, 17, 1156–1162. [Google Scholar] [CrossRef] [Green Version]
  9. Ligabue-Braun, R.; Verli, H.; Carlini, C.R. Venomous Mammals: A Review. Toxicon 2012, 59, 680–695. [Google Scholar] [CrossRef]
  10. Pekár, S.; Bočánek, O.; Michálek, O.; Petráková, L.; Haddad, C.R.; Šedo, O.; Zdráhal, Z. Venom Gland Size and Venom Complexity—Essential Trophic Adaptations of Venomous Predators: A Case Study Using Spiders. Mol. Ecol. 2018, 27, 4257–4269. [Google Scholar] [CrossRef]
  11. Giorgianni, M.W.; Dowell, N.L.; Griffin, S.; Kassner, V.A.; Selegue, J.E.; Carroll, S.B. The Origin and Diversification of a Novel Protein Family in Venomous Snakes. Proc. Natl. Acad. Sci. USA 2020, 117, 10911–10920. [Google Scholar] [CrossRef] [PubMed]
  12. Sanz, L.; Quesada-Bernat, S.; Ramos, T.; Casais-e-Silva, L.L.; Corrêa-Netto, C.; Silva-Haad, J.J.; Sasa, M.; Lomonte, B.; Calvete, J.J. New Insights into the Phylogeographic Distribution of the 3FTx/PLA2 Venom Dichotomy across Genus Micrurus in South America. J. Proteom. 2019, 200, 90–101. [Google Scholar] [CrossRef] [PubMed]
  13. Facente, J.; McGivern, J.J.; Seavy, M.; Wray, K.P.; 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]
  14. Mackessy, S.P.; Saviola, A.J. Understanding Biological Roles of Venoms among the Caenophidia: The Importance of Rear-Fanged Snakes. Integr. Comp. Biol. 2016, 56, 1004–1021. [Google Scholar] [CrossRef] [PubMed]
  15. Mackessy, S.P. Handbook of Venoms and Toxins of Reptiles, 1st ed.; Mackessy, S.P., Ed.; CRC Press: Boca Raton, FL, USA, 2016; ISBN 9780429186394. [Google Scholar]
  16. Casewell, N.R.; Jackson, T.N.W.; Laustsen, A.H.; Sunagar, K. Causes and Consequences of Snake Venom Variation. Trends Pharmacol. Sci. 2020, 41, 570–581. [Google Scholar] [CrossRef] [PubMed]
  17. Jackson, T.N.W.; Koludarov, I.; Ali, S.A.; Dobson, J.; Zdenek, C.N.; Dashevsky, D.; Op Den Brouw, B.; Masci, P.P.; Nouwens, A.; Josh, P.; et al. Rapid Radiations and the Race to Redundancy: An Investigation of the Evolution of Australian Elapid Snake Venoms. Toxins 2016, 8, 309. [Google Scholar] [CrossRef]
  18. Pyron, R.A.; Burbrink, F.T.; Wiens, J.J. A Phylogeny and Revised Classification of Squamata, Including 4161 Species of Lizards and Snakes. BMC Evol. Biol. 2013, 13, 93. [Google Scholar] [CrossRef] [Green Version]
  19. Junqueira-de-Azevedo, I.L.M.; Campos, P.F.; Ching, A.T.C.; Mackessy, S.P. Colubrid Venom Composition: An Omics Perspective. Toxins 2016, 8, 230. [Google Scholar] [CrossRef] [Green Version]
  20. Jackson, T.N.W.; Young, B.; Underwood, G.; McCarthy, C.J.; Kochva, E.; Vidal, N.; van der Weerd, L.; Nabuurs, R.; Dobson, J.; Whitehead, D.; et al. Endless Forms Most Beautiful: The Evolution of Ophidian Oral Glands, Including the Venom System, and the Use of Appropriate Terminology for Homologous Structures. Zoomorphology 2017, 136, 107–130. [Google Scholar] [CrossRef]
  21. McGivern, J.J.; Wray, K.P.; Margres, M.J.; Couch, M.E.; Mackessy, S.P.; Rokyta, D.R. RNA-Seq and High-Definition Mass Spectrometry Reveal the Complex and Divergent Venoms of Two Rear-Fanged Colubrid Snakes. BMC Genom. 2014, 15, 1061. [Google Scholar] [CrossRef] [Green Version]
  22. 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]
  23. Hofmann, E.P.; Rautsaw, R.M.; Mason, A.J.; Strickland, J.L.; Parkinson, C.L. Duvernoy’s Gland Transcriptomics of the Plains Black-Headed Snake, Tantilla Nigriceps (Squamata, Colubridae): Unearthing the Venom of Small Rear-Fanged Snakes. Toxins 2021, 13, 336. [Google Scholar] [CrossRef] [PubMed]
  24. Modahl, C.M.; Frietze, S.; Mackessy, S.P. Transcriptome-Facilitated Proteomic Characterization of Rear-Fanged Snake Venoms Reveal Abundant Metalloproteinases with Enhanced Activity. J. Proteom. 2018, 187, 223–234. [Google Scholar] [CrossRef] [PubMed]
  25. Pla, D.; Sanz, L.; Whiteley, G.; Wagstaff, S.C.; Harrison, R.A.; Casewell, N.R.; Calvete, J.J. What Killed Karl Patterson Schmidt? Combined Venom Gland Transcriptomic, Venomic and Antivenomic Analysis of the South African Green Tree Snake (the Boomslang), Dispholidus Typus. Biochim. Biophys. Acta Gen. Subj. 2017, 1861, 814–823. [Google Scholar] [CrossRef] [Green Version]
  26. Modahl, C.M.; Mrinalini; Frietze, S.; Mackessy, S.P. Adaptive Evolution of Distinct Prey-Specific Toxin Genes in Rear-Fanged Snake Venom. Proc. R. Soc. B Biol. Sci. 2018, 285, 20181003. [Google Scholar] [CrossRef]
  27. Mackessy, S.P.; Bryan, W.; Smith, C.F.; Lopez, K.; Fernández, J.; Bonilla, F.; Camacho, E.; Sasa, M.; Lomonte, B. Venomics of the Central American Lyre Snake Trimorphodon Quadruplex (Colubridae: Smith, 1941) from Costa Rica. J. Proteom. 2020, 220, 103778. [Google Scholar] [CrossRef]
  28. Schramer, T.D.; Rautsaw, R.M.; Bayona-Serrano, J.D.; Nystrom, G.S.; West, T.R.; Ortiz-Medina, J.A.; Sabido-Alpuche, B.; Meneses-Millán, M.; Borja, M.; Junqueira-de-Azevedo, I.L.M.; et al. An Integrative View of the Toxic Potential of Conophis Lineatus (Dipsadidae: Xenodontinae), a Medically Relevant Rear-Fanged Snake. Toxicon 2022, 205, 38–52. [Google Scholar] [CrossRef]
  29. Campos, P.F.; Andrade-Silva, D.; Zelanis, A.; Leme, A.F.P.; Rocha, M.M.T.; Menezes, M.C.; Serrano, S.M.T.; Junqueira-De-Azevedo, I.D.L.M. Trends in the Evolution of Snake Toxins Underscored by an Integrative Omics Approach to Profile the Venom of the Colubrid Phalotris Mertensi. Genome Biol. Evol. 2016, 8, 2266–2287. [Google Scholar] [CrossRef] [Green Version]
  30. Ching, A.T.C.; Rocha, M.M.T.; Leme, A.F.P.; Pimenta, D.C.; de Fátima, D.; Furtado, M.; Serrano, S.M.T.; Ho, P.L.; Junqueira-de-Azevedo, I.L.M. Some Aspects of the Venom Proteome of the Colubridae Snake Philodryas Olfersii Revealed from a Duvernoy’s (Venom) Gland Transcriptome. FEBS Lett. 2006, 580, 4417–4422. [Google Scholar] [CrossRef] [Green Version]
  31. Ching, A.T.C.; Leme, A.F.P.; Zelanis, A.; Rocha, M.M.T.; Furtado, M.D.F.D.; Silva, D.A.; Trugilho, M.R.O.; Da Rocha, S.L.G.; Perales, J.; Ho, P.L.; et al. Venomics Profiling of Thamnodynastes Strigatus Unveils Matrix Metalloproteinases and Other Novel Proteins Recruited to the Toxin Arsenal of Rear-Fanged Snakes. J. Proteome Res. 2012, 11, 1152–1162. [Google Scholar] [CrossRef]
  32. Lemoine, K.; Girón, M.E.; Aguilar, I.; Navarrete, L.F.; Rodríguez-Acosta, A. Proteolytic, Hemorrhagic, and Neurotoxic Activities Caused by Leptodeira annulata ashmeadii (Serpentes: Colubridae) Duvernoy’s Gland Secretion. Wilderness Environ. Med. 2004, 15, 82–89. [Google Scholar] [CrossRef]
  33. Sánchez, M.N.; Gonzalez, K.Y.; Sciani, J.M.; Gritti, M.A.; Maruñak, S.L.; Tavares, F.L.; Teibler, G.P.; Peichoto, M.E. First Insights into the Biochemical and Toxicological Characterization of Venom from the Banded Cat-Eyed Snake Leptodeira annulata pulchriceps. Comp. Biochem. Physiol. Part C Toxicol. Pharmacol. 2021, 239, 108897. [Google Scholar] [CrossRef] [PubMed]
  34. Torres-Bonilla, K.A.; Schezaro-Ramos, R.; Floriano, R.S.; Rodrigues-Simioni, L.; Bernal-Bautista, M.H.; da Cruz-Höfling, M.A. Biological Activities of Leptodeira annulata (Banded Cat-Eyed Snake) Venom on Vertebrate Neuromuscular Preparations. Toxicon 2016, 119, 345–351. [Google Scholar] [CrossRef] [PubMed]
  35. Torres-Bonilla, K.A.; Panunto, P.C.; Pereira, B.B.; Zambrano, D.F.; Herrán-Medina, J.; Bernal, M.H.; Hyslop, S. Toxinological Characterization of Venom from Leptodeira annulata (Banded Cat-Eyed Snake; Dipsadidae, Imantodini). Biochimie 2020, 174, 171–188. [Google Scholar] [CrossRef]
  36. Estrella, A.; Sánchez, E.E.; Galán, J.A.; Tao, W.A.; Guerrero, B.; Navarrete, L.F.; Rodríguez-Acosta, A. Characterization of Toxins from the Broad-Banded Water Snake Helicops angulatus (Linnaeus, 1758): Isolation of a Cysteine-Rich Secretory Protein, Helicopsin. Arch. Toxicol. 2011, 85, 305–313. [Google Scholar] [CrossRef]
  37. Estrella, A.; Rodríguez-Torres, A.; Serna, L.; Navarrete, L.F. Is the South American Water Snake Helicops angulatus (Linnaeus, 1758) (Ddipsadidae:Xenodontinae) Venomous? Herpetotropicos 2012, 5, 79–84. [Google Scholar]
  38. Gutiérrez, J.M.; Calvete, J.J.; Habib, A.G.; Harrison, R.A.; Williams, D.J.; Warrell, D.A. Snakebite Envenoming. Nat. Rev. Dis. Prim. 2017, 3, 17063. [Google Scholar] [CrossRef]
  39. Lomonte, B.; Rey-Suárez, P.; Fernández, J.; Sasa, M.; Pla, D.; Vargas, N.; Bénard-Valle, M.; Sanz, L.; Corrêa-Netto, C.; Núñez, V.; et al. Venoms of Micrurus Coral Snakes: Evolutionary Trends in Compositional Patterns Emerging from Proteomic Analyses. Toxicon 2016, 122, 7–25. [Google Scholar] [CrossRef]
  40. Amazonas, D.R.; Portes-Junior, J.A.; Nishiyama, M.Y., Jr.; Nicolau, C.A.; Chalkidis, H.M.; Mourão, R.H.V.; Grazziotin, F.G.; Rokyta, D.R.; Gibbs, H.L.; Valente, R.H.; et al. Molecular Mechanisms Underlying Intraspecific Variation in Snake Venom. J. Proteom. 2018, 181, 60–72. [Google Scholar] [CrossRef]
  41. Sanz, L.; Quesada-Bernat, S.; Pérez, A.; De Morais-Zani, K.; SantEänna, S.S.; Hatakeyama, D.M.; Tasima, L.J.; De Souza, M.B.; Kayano, A.M.; Zavaleta, A.; et al. Danger in the Canopy. Comparative Proteomics and Bioactivities of the Venoms of the South American Palm Pit Viper Bothrops bilineatus Subspecies bilineatus and smaragdinus and Antivenomics of B. b. bilineatus (Rondônia) Venom against the Brazilian Pentabothropic Antivenom. J. Proteome Res. 2020, 19, 3518–3532. [Google Scholar] [CrossRef]
  42. Peichoto, M.E.; Tavares, F.L.; Santoro, M.L.; MacKessy, S.P. Venom Proteomes of South and North American Opisthoglyphous (Colubridae and Dipsadidae) Snake Species: A Preliminary Approach to Understanding Their Biological Roles. Comp. Biochem. Physiol. Part D Genom. Proteom. 2012, 7, 361–369. [Google Scholar] [CrossRef] [PubMed]
  43. Modahl, C.M.; Mackessy, S.P. Venoms of Rear-Fanged Snakes: New Proteins and Novel Activities. Front. Ecol. Evol. 2019, 7, 279. [Google Scholar] [CrossRef] [Green Version]
  44. Grundler, M.C. SquamataBase: A Natural History Database and R Package for Comparative Biology of Snake Feeding Habits. Biodivers. Data J. 2020, 8, e49943. [Google Scholar] [CrossRef] [PubMed]
  45. Duellman, W.E. Cusco Amazonico: The Lives of Amphibians and Reptiles in an Amazonian Rainforest; Comstock Pub. Associates: Ithaca, NY, USA, 2005; Volume 43, ISBN 9780801439971. [Google Scholar]
  46. Arlinghaus, F.T.; Eble, J.A. C-Type Lectin-like Proteins from Snake Venoms. Toxicon 2012, 60, 512–519. [Google Scholar] [CrossRef]
  47. Ogawa, T.; Chijiwa, T.; Oda-Ueda, N.; Ohno, M. Molecular Diversity and Accelerated Evolution of C-Type Lectin-like Proteins from Snake Venom. Toxicon 2005, 45, 1–14. [Google Scholar] [CrossRef]
  48. Xie, B.; Dashevsky, D.; Rokyta, D.; Ghezellou, P.; Fathinia, B.; Shi, Q.; Richardson, M.K.; Fry, B.G. Dynamic Genetic Differentiation Drives the Widespread Structural and Functional Convergent Evolution of Snake Venom Proteinaceous Toxins. BMC Biol. 2022, 20, 4. [Google Scholar] [CrossRef]
  49. Sanz, L.; Pérez, A.; Quesada-Bernat, S.; Diniz-Sousa, R.; Calderón, L.A.; Soares, A.M.; Calvete, J.J.; Caldeira, C.A.S. Venomics and Antivenomics of the Poorly Studied Brazil’s Lancehead, Bothrops brazili (Hoge, 1954), from the Brazilian State of Pará. J. Venom. Anim. Toxins Incl. Trop. Dis. 2020, 26, 20190103. [Google Scholar] [CrossRef]
  50. Strickland, J.L.; Smith, C.F.; Mason, A.J.; Schield, D.R.; Borja, M.; Castañeda-Gaytán, G.; Spencer, C.L.; Smith, L.L.; Trápaga, A.; Bouzid, N.M.; et al. Evidence for Divergent Patterns of Local Selection Driving Venom Variation in Mojave Rattlesnakes (Crotalus scutulatus). Sci. Rep. 2018, 8, 17622. [Google Scholar] [CrossRef]
  51. Weese, D.A.; Duda, T.F. Effects of Predator-Prey Interactions on Predator Traits: Differentiation of Diets and Venoms of a Marine Snail. Toxins 2019, 11, 299. [Google Scholar] [CrossRef] [Green Version]
  52. Rotenberg, D.; Bamberger, E.S.; Kochva, E. Studies on Ribonucleic Acid Synthesis in the Venom Glands of Vipera Palaestinae (Ophidia, Reptilia). Biochem. J. 1971, 121, 609–612. [Google Scholar] [CrossRef] [Green Version]
  53. Mora-Obando, D.; Salazar-Valenzuela, D.; Pla, D.; Lomonte, B.; Guerrero-Vargas, J.A.; Ayerbe, S.; Gibbs, H.L.; Calvete, J.J. Venom Variation in Bothrops asper Lineages from North-Western South America. J. Proteom. 2020, 229, 103945. [Google Scholar] [CrossRef] [PubMed]
  54. Mackessy, S.P. Venom Ontogeny in the Pacific Rattlesnakes Crotalus Viridis Helleri and C. v. Oreganus. Copeia 1988, 1988, 92. [Google Scholar] [CrossRef]
  55. Haney, R.A.; Clarke, T.H.; Gadgil, R.; Fitzpatrick, R.; Hayashi, C.Y.; Ayoub, N.A.; Garb, J.E. Effects of Gene Duplication, Positive Selection, and Shifts in Gene Expression on the Evolution of the Venom Gland Transcriptome in Widow Spiders. Genome Biol. Evol. 2016, 8, 228–242. [Google Scholar] [CrossRef] [Green Version]
  56. Barua, A.; Mikheyev, A.S. Many Options, Few Solutions: Over 60 My Snakes Converged on a Few Optimal Venom Formulations. Mol. Biol. Evol. 2019, 36, 1964–1974. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data. Available online: http://www.Bioinformatics.Babraham.Ac.Uk/Projects/Fastqc/ (accessed on 10 March 2021).
  58. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A Flexible Trimmer for Illumina Sequence Data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef] [Green Version]
  59. Haas, B.J.; Papanicolaou, A.; Yassour, M.; Grabherr, M.; Blood, P.D.; Bowden, J.; Couger, M.B.; Eccles, D.; Li, B.; Lieber, M.; et al. De Novo Transcript Sequence Reconstruction from RNA-Seq Using the Trinity Platform for Reference Generation and Analysis. Nat. Protoc. 2013, 8, 1494–1512. [Google Scholar] [CrossRef]
  60. Rokyta, D.R.; Lemmon, A.R.; Margres, M.J.; Aronow, K. The Venom-Gland Transcriptome of the Eastern Diamondback Rattlesnake (Crotalus adamanteus). BMC Genom. 2012, 13, 312. [Google Scholar] [CrossRef] [Green Version]
  61. Holding, M.L.; Margres, M.J.; Mason, A.J.; Parkinson, C.L.; Rokyta, D.R. Evaluating the Performance of de Novo Assembly Methods for Venom-Gland Transcriptomics. Toxins 2018, 10, 249. [Google Scholar] [CrossRef] [Green Version]
  62. Zhang, J.; Kobert, K.; Flouri, T.; Stamatakis, A. PEAR: A Fast and Accurate Illumina Paired-End ReAd MergeR. Bioinformatics 2014, 30, 614–620. [Google Scholar] [CrossRef] [Green Version]
  63. Yang, Y.; Smith, S.A. Orthology Inference in Nonmodel Organisms Using Transcriptomes and Low-Coverage Genomes: Improving Accuracy and Matrix Occupancy for Phylogenomics. Mol. Biol. Evol. 2014, 31, 3081–3092. [Google Scholar] [CrossRef]
  64. Morales-Briones, D.F.; Kadereit, G.; Tefarikis, D.T.; Moore, M.J.; Smith, S.A.; Brockington, S.F.; Timoneda, A.; Yim, W.C.; Cushman, J.C.; Yang, Y. Disentangling Sources of Gene Tree Discordance in Phylogenomic Data Sets: Testing Ancient Hybridizations in Amaranthaceae s.l. Syst. Biol. 2021, 70, 219–235. [Google Scholar] [CrossRef] [PubMed]
  65. Smith-Unna, R.; Boursnell, C.; Patro, R.; Hibberd, J.M.; Kelly, S. TransRate: Reference-Free Quality Assessment of de Novo Transcriptome Assemblies. Genome Res. 2016, 26, 1134–1144. [Google Scholar] [CrossRef] [Green Version]
  66. Yang, Y.; Smith, S.A. Optimizing de Novo Assembly of Short-Read RNA-Seq Data for Phylogenomics. BMC Genom. 2013, 14, 328. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  67. Davidson, N.M.; Oshlack, A. Corset: Enabling Differential Gene Expression Analysis for de Novo Assembled Transcriptomes. Genome Biol. 2014, 15, 410. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  68. Patro, R.; Duggal, G.; Love, M.I.; Irizarry, R.A.; Kingsford, C. Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression. Nat. Methods 2017, 14, 417–419. [Google Scholar] [CrossRef] [Green Version]
  69. Available online: http://Transdecoder.Github.Io/ (accessed on 12 March 2021).
  70. Camacho, C.; Coulouris, G.; Avagyan, V.; Ma, N.; Papadopoulos, J.; Bealer, K.; Madden, T.L. BLAST+: Architecture and Applications. BMC Bioinform. 2009, 10, 421. [Google Scholar] [CrossRef] [Green Version]
  71. Fu, L.; Niu, B.; Zhu, Z.; Wu, S.; Li, W. CD-HIT: Accelerated for Clustering the next-Generation Sequencing Data. Bioinformatics 2012, 28, 3150–3152. [Google Scholar] [CrossRef]
  72. Li, B.; Ruotti, V.; Stewart, R.M.; Thomson, J.A.; Dewey, C.N. RNA-Seq Gene Expression Estimation with Read Mapping Uncertainty. Bioinformatics 2009, 26, 493–500. [Google Scholar] [CrossRef] [Green Version]
  73. Langmead, B.; Salzberg, S.L. Fast Gapped-Read Alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef] [Green Version]
  74. 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]
  75. Margres, M.J.; Aronow, K.; Loyacano, J.; Rokyta, D.R. The Venom-Gland Transcriptome of the Eastern Coral Snake (Micrurus fulvius) Reveals High Venom Complexity in the Intragenomic Evolution of Venoms. BMC Genom. 2013, 14, 531. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  76. Vonk, F.J.; Casewell, N.R.; Henkel, C.V.; Heimberg, A.M.; Jansen, H.J.; McCleary, R.J.R.; Kerkkamp, H.M.E.; Vos, R.A.; Guerreiro, I.; Calvete, J.J.; et al. The King Cobra Genome Reveals Dynamic Gene Evolution and Adaptation in the Snake Venom System. Proc. Natl. Acad. Sci. USA 2013, 110, 20651–20656. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Bayona-Serrano, J.D.; Viala, V.L.; Rautsaw, R.M.; Schramer, T.D.; Barros-Carvalho, G.A.; Nishiyama, M.Y.; Freitas-De-Sousa, L.A.; Moura-Da-Silva, A.M.; Parkinson, C.L.; Grazziotin, F.G.; et al. Replacement and Parallel Simplification of Nonhomologous Proteinases Maintain Venom Phenotypes in Rear-Fanged Snakes. Mol. Biol. Evol. 2020, 37, 3563–3575. [Google Scholar] [CrossRef] [PubMed]
  78. R Core Team. R: A Language and Environment for Statistical Computing. Available online: https://www.r-project.org/ (accessed on 3 June 2020).
  79. Johnson, M.G.; Gardner, E.M.; Liu, Y.; Medina, R.; Goffinet, B.; Shaw, A.J.; Zerega, N.J.C.; Wickett, N.J. HybPiper: Extracting Coding Sequence and Introns for Phylogenetics from High-throughput Sequencing Reads Using Target Enrichment. Appl. Plant Sci. 2016, 4, 1600016. [Google Scholar] [CrossRef] [Green Version]
  80. Dong, S.; Kumazawa, Y. Complete Mitochondrial DNA Sequences of Six Snakes: Phylogenetic Relationships and Molecular Evolution of Genomic Features. J. Mol. Evol. 2005, 61, 12–22. [Google Scholar] [CrossRef]
  81. Kumazawa, Y.; Ota, H.; Nishida, M.; Ozawa, T. The Complete Nucleotide Sequence of a Snake (Dinodon Semicarinatus) Mitochondrial Genome with Two Identical Control Regions. Genetics 1998, 150, 313–329. [Google Scholar] [CrossRef]
  82. Jiang, Z.J.; Castoe, T.A.; Austin, C.C.; Burbrink, F.T.; Herron, M.D.; McGuire, J.A.; Parkinson, C.L.; Pollock, D.D. Comparative Mitochondrial Genomics of Snakes: Extraordinary Substitution Rate Dynamics and Functionality of the Duplicate Control Region. BMC Evol. Biol. 2007, 7, 123. [Google Scholar] [CrossRef] [Green Version]
  83. Yan, J.; Li, H.; Zhou, K. Evolution of the Mitochondrial Genome in Snakes: Gene Rearrangements and Phylogenetic Relationships. BMC Genom. 2008, 9, 569. [Google Scholar] [CrossRef] [Green Version]
  84. Kumazawa, Y.; Nishida, M. Complete Mitochondrial DNA Sequences of the Green Turtle and Blue- Tailed Mole Skink: Statistical Evidence for Archosaurian Affinity of Turtles. Mol. Biol. Evol. 1999, 16, 784–792. [Google Scholar] [CrossRef] [Green Version]
  85. Macey, J.R.; Papenfuss, T.J.; Kuehl, J.V.; Fourcade, H.M.; Boore, J.L. Phylogenetic Relationships among Amphisbaenian Reptiles Based on Complete Mitochondrial Genomic Sequences. Mol. Phylogenet. Evol. 2004, 33, 22–31. [Google Scholar] [CrossRef] [Green Version]
  86. Slater, G.S.C.; Birney, E. Automated Generation of Heuristics for Biological Sequence Comparison. BMC Bioinform. 2005, 6, 31. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Cock, P.J.A.; Antao, T.; Chang, J.T.; Chapman, B.A.; Cox, C.J.; Dalke, A.; Friedberg, I.; Hamelryck, T.; Kauff, F.; Wilczynski, B.; et al. Biopython: Freely Available Python Tools for Computational Molecular Biology and Bioinformatics. Bioinformatics 2009, 25, 1422–1423. [Google Scholar] [CrossRef] [PubMed]
  88. Li, H.; Durbin, R. Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform. Bioinformatics 2009, 25, 1754–1760. [Google Scholar] [CrossRef] [Green Version]
  89. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R. The Sequence Alignment/Map Format and SAMtools. Bioinformatics 2009, 25, 2078–2079. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  90. Tange, O. GNU Parallel 20150322 (‘Hellwig’). Login USENIX Mag. 2015, 436, 42–47. [Google Scholar] [CrossRef]
  91. Bankevich, A.; Nurk, S.; Antipov, D.; Gurevich, A.A.; Dvorkin, M.; Kulikov, A.S.; Lesin, V.M.; Nikolenko, S.I.; Pham, S.; Prjibelski, A.D.; et al. SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing. J. Comput. Biol. 2012, 19, 455–477. [Google Scholar] [CrossRef] [Green Version]
  92. Katoh, K.; Standley, D.M. MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability. Mol. Biol. Evol. 2013, 30, 772–780. [Google Scholar] [CrossRef] [Green Version]
  93. Brown, J.W.; Walker, J.F.; Smith, S.A. Phyx: Phylogenetic Tools for Unix. Bioinformatics 2017, 33, 1886–1888. [Google Scholar] [CrossRef] [Green Version]
  94. Nguyen, L.T.; Schmidt, H.A.; Von Haeseler, A.; Minh, B.Q. IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies. Mol. Biol. Evol. 2015, 32, 268–274. [Google Scholar] [CrossRef]
  95. Shannon, C.E. A Mathematical Theory of Communication. Bell Syst. Tech. J. 1948, 27, 379–423. [Google Scholar] [CrossRef] [Green Version]
  96. Bray, J.R.; Curtis, J.T. An Ordination of the Upland Forest Communities of Southern Wisconsin. Ecol. Monogr. 1957, 27, 325–349. [Google Scholar] [CrossRef]
  97. Callahan, S.; Crowe-Riddell, J.M.; Nagesan, R.S.; Gray, J.A.; Rabosky, A.R.D. A Guide for Optimal Iodine Staining and High-Throughput DiceCT Scanning in Snakes. Ecol. Evol. 2021, 11, 11587–11603. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Representative venom delivery system for each genus and venom transcript counts of individuals. (A) Soft-tissue scan of the venom delivery system for each genus mapped onto a simplified phylogeny developed in this study. (B) Isolated maxilla bone from each representative genus showing position of fangs (arrow). (C) Unique toxin transcripts recovered from each Helicops individual. (D) Unique toxin transcripts recovered from each Leptodeira individual. (E) Unique toxin transcripts recovered from each Micrurus individual. (F) Unique toxin transcripts recovered from each Bothrops individual. MUSM = Museo de Historia Natural, Universidad Nacional Mayor de San Marco, UMMZ = University of Michigan Museum of Zoology, SVMP = snake venom metalloproteinase, PLA2 = phospholipase A2, 3FTx = three-finger toxin, CTL = C-type lectin, SVSP = snake venom serine proteinase, BPP = bradykinin-potentiating peptides, LAAO = L-amino acid oxidase, CRiSP = cystine rich secretory protein, Kunitz = Kunitz-type serine protease. Micro-CT scans of specimens vouchered in the University of Michigan Museum of Zoology (UMMZ) and Museo de Historia Natural de la Universidad Nacional Major de San Marcos (MUSM). We deposited these micro-CT scans used for venom and fang morphology for public access in the Morphosource ‘Scan All Snakes’ Project ID 00000C374 (https://www.morphosource.org/Detail/ProjectDetail/Show/project_id/374).
Figure 1. Representative venom delivery system for each genus and venom transcript counts of individuals. (A) Soft-tissue scan of the venom delivery system for each genus mapped onto a simplified phylogeny developed in this study. (B) Isolated maxilla bone from each representative genus showing position of fangs (arrow). (C) Unique toxin transcripts recovered from each Helicops individual. (D) Unique toxin transcripts recovered from each Leptodeira individual. (E) Unique toxin transcripts recovered from each Micrurus individual. (F) Unique toxin transcripts recovered from each Bothrops individual. MUSM = Museo de Historia Natural, Universidad Nacional Mayor de San Marco, UMMZ = University of Michigan Museum of Zoology, SVMP = snake venom metalloproteinase, PLA2 = phospholipase A2, 3FTx = three-finger toxin, CTL = C-type lectin, SVSP = snake venom serine proteinase, BPP = bradykinin-potentiating peptides, LAAO = L-amino acid oxidase, CRiSP = cystine rich secretory protein, Kunitz = Kunitz-type serine protease. Micro-CT scans of specimens vouchered in the University of Michigan Museum of Zoology (UMMZ) and Museo de Historia Natural de la Universidad Nacional Major de San Marcos (MUSM). We deposited these micro-CT scans used for venom and fang morphology for public access in the Morphosource ‘Scan All Snakes’ Project ID 00000C374 (https://www.morphosource.org/Detail/ProjectDetail/Show/project_id/374).
Toxins 14 00489 g001
Figure 2. Expression of toxin gene families from venom gland transcriptomes (center) and overall venom transcriptome diversity (right) mapped to phylogeny inferred from mitochondrial gene sequences (left). Note that these data suggest there is generally higher similarity in venom profiles among individuals within genera than among genera for both metrics. Bothrops and Micrurus species are front-fanged, while Leptodeira and Helicops species are rear-fanged. SVMP = snake venom metalloproteinase, PLA2 = phospholipase A2, 3FTx = three-finger toxin, CTL = C-type lectin, SVSP = snake venom serine proteinase, BPP = bradykinin-potentiating peptides, LAAO = L-amino acid oxidase, CRiSP = cystine-rich secretory protein, Kunitz = Kunitz-type serine protease.
Figure 2. Expression of toxin gene families from venom gland transcriptomes (center) and overall venom transcriptome diversity (right) mapped to phylogeny inferred from mitochondrial gene sequences (left). Note that these data suggest there is generally higher similarity in venom profiles among individuals within genera than among genera for both metrics. Bothrops and Micrurus species are front-fanged, while Leptodeira and Helicops species are rear-fanged. SVMP = snake venom metalloproteinase, PLA2 = phospholipase A2, 3FTx = three-finger toxin, CTL = C-type lectin, SVSP = snake venom serine proteinase, BPP = bradykinin-potentiating peptides, LAAO = L-amino acid oxidase, CRiSP = cystine-rich secretory protein, Kunitz = Kunitz-type serine protease.
Toxins 14 00489 g002
Table 1. Summary of venom gland transcriptome toxin profiles of several species of rear-fanged colubrid snakes. SVMP = snake venom metalloproteinase, 3FTx = three-finger toxin, CRiSP = cystine rich secretory protein, svMMP = snake venom matrix metalloproteinase, CTL = C-type lectins, Kunitz = Kunitz-type serine protease.
Table 1. Summary of venom gland transcriptome toxin profiles of several species of rear-fanged colubrid snakes. SVMP = snake venom metalloproteinase, 3FTx = three-finger toxin, CRiSP = cystine rich secretory protein, svMMP = snake venom matrix metalloproteinase, CTL = C-type lectins, Kunitz = Kunitz-type serine protease.
SubfamilySpeciesMajor Venom Component(s)Reference
ColubrinaeAhaetulla prasinaSVMPs, 3FTxs[24]
Boiga irregularis3FTxs, SVMPs[21]
Dispholidus typusSVMPs[25]
Spilotes sulphureus3FTxs[26]
Tantilla nigriceps3FTxs, CRiSPs, SVMPs[23]
Trimorphodon quadruplex3FTxs, SVMPs[27]
DipsadinaeBorikenophis portoricensisSVMPs[24]
Conophis lineatussvMMPs[28]
Hypsiglena sp. SVMPs, CRiSPs[21]
Phalotris mertensiKunitzs, SVMPs, CTLs[29]
Philodryas olfersiiSVMPs, CNPs[30]
Thamnodynastes strigatussvMMPs[31]
Table 2. Specimen information for samples sequenced in this study. Numbers at the end of species names are field codes used to identify individuals. MUSM = Museo de Historia Natural, Universidad Nacional Mayor de San Marcos, UMMZ = University of Michigan Museum of Zoology, EBLA = Estación Biológica Los Amigos, EBMS = Estación Biológica Madre Selva, LBM = Las Brisas del Mogotón, EBVC = Estación Biológica Villa Carmen, RB = Refugio Bartola, mm = millimeters, g = grams, F = female, M = male, J = juvenile, A = adult.
Table 2. Specimen information for samples sequenced in this study. Numbers at the end of species names are field codes used to identify individuals. MUSM = Museo de Historia Natural, Universidad Nacional Mayor de San Marcos, UMMZ = University of Michigan Museum of Zoology, EBLA = Estación Biológica Los Amigos, EBMS = Estación Biológica Madre Selva, LBM = Las Brisas del Mogotón, EBVC = Estación Biológica Villa Carmen, RB = Refugio Bartola, mm = millimeters, g = grams, F = female, M = male, J = juvenile, A = adult.
FamilyTaxonMuseum Accession No.Date CapturedStation, CountrySVL (mm)Mass (g)SexAge
Viperidae Bothrops atrox 0365MUSM 3572121 March 2016EBLA, Peru58981FJ
Bothrops bilineatus 0065UMMZ 24508411 March 2016EBLA, Peru74485FA
Bothrops brazili 1278MUSM 369221 December 2016EBLA, Peru60676MA
ElapidaeMicrurus annellatus 3275UMMZ 24845026 November 2018EBLA, Peru49718.11FA
Micrurus hemprichii 1810UMMZ 24685718 January 2017EBMS, Peru74086MA
Micrurus lemniscatus 0249UMMZ 24508216 March 2016EBLA, Peru71565MA
Micrurus lemniscatus 0336MUSM 3590521 March 2016EBLA, Peru72550FA
Micrurus nigrocinctus 3053UMMZ 24714222 May 2018LBM, Nicaragua71764.8FA
Micrurus obscurus 0665UMMZ 2468597 November 2016EBVC, Peru2615.19MJ
Micrurus obscurus 1054UMMZ 24686022 November 2016EBLA, Peru77581MA
ColubridaeHelicops angulatus 0143UMMZ 24505313 March 2016EBLA, Peru37348.36FA
Helicops angulatus 3440UMMZ 2488792 December 2018EBLA, Peru41160FA
Helicops angulatus 3559MUSM 398269 December 2018EBLA, Peru30724.19FA
Helicops leopardinus 1812UMMZ 24680818 January 2017EBMS, Peru685220FA
Helicops polylepis 1932UMMZ 24680918 January 2017EBMS, Peru823600FA
Leptodeira annulata 0468UMMZ 24505924 March 2016EBLA, Peru46318.56MA
Leptodeira annulata 0497UMMZ 24506027 March 2016EBLA, Peru59038.02FA
Leptodeira rhombifera 3241UMMZ 24709812 June 2018Tecomapa, Nicaragua665169.5FA
Leptodeira septentrionalis 3176UMMZ 2470993 June 2018RB, Nicaragua654113.2FA
Table 3. Library preparation, sequencing platform, sequencing output, and percent of the transcriptome that was comprised of toxin transcripts for individuals used in the study.
Table 3. Library preparation, sequencing platform, sequencing output, and percent of the transcriptome that was comprised of toxin transcripts for individuals used in the study.
FamilyTaxonLibrary PreparationIllumina PlatformReads PairsPercent Toxin Expression
Viperidae Bothrops atrox 0365TruSeq RNASeqHiSeq 400022,392,18247.97%
Bothrops bilineatus 0065NEBNext Ultra IINovaSeq 600017,985,94574.53%
Bothrops brazili 1278NEBNext Ultra IINovaSeq 600014,816,99847.25%
ElapidaeMicrurus annellatus 3275NEBNext Ultra IINovaSeq 600016,398,04639.12%
Micrurus hemprichii 1810NEBNext Ultra IINovaSeq 600019,149,98633.68%
Micrurus lemniscatus 0249NEBNext Ultra IINovaSeq 600016,413,19019.58%
Micrurus lemniscatus 0336TruSeq RNASeqHiSeq 400024,938,73253.04%
Micrurus nigrocinctus 3053NEBNext Ultra IINovaSeq 600018,981,69263.38%
Micrurus obscurus 0665NEBNext Ultra IINovaSeq 600015,955,90440.44%
Micrurus obscurus 1054NEBNext Ultra IINovaSeq 600016,791,60148.82%
ColubridaeHelicops angulatus 0143TruSeq RNASeqHiSeq 400023,374,95817.15%
Helicops angulatus 3440NEBNext Ultra IINovaSeq 600017,274,73531.03%
Helicops angulatus 3559NEBNext Ultra IINovaSeq 600015,797,92119.57%
Helicops leopardinus 1812NEBNext Ultra IINovaSeq 600017,894,32252.98%
Helicops polylepis 1932NEBNext Ultra IINovaSeq 600012,936,85891.81%
Leptodeira annulata 0468NEBNext Ultra IINovaSeq 600019,191,19327.03%
Leptodeira annulata 0497TruSeq RNASeqHiSeq 400020,844,57936.83%
Leptodeira rhombifera 3241NEBNext Ultra IINovaSeq 600017,838,00041.59%
Leptodeira septentrionalis 3176NEBNext Ultra IINovaSeq 600017,147,03146.39%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Cerda, P.A.; Crowe-Riddell, J.M.; Gonçalves, D.J.P.; Larson, D.A.; Duda, T.F., Jr.; Davis Rabosky, A.R. Divergent Specialization of Simple Venom Gene Profiles among Rear-Fanged Snake Genera (Helicops and Leptodeira, Dipsadinae, Colubridae). Toxins 2022, 14, 489. https://doi.org/10.3390/toxins14070489

AMA Style

Cerda PA, Crowe-Riddell JM, Gonçalves DJP, Larson DA, Duda TF Jr., Davis Rabosky AR. Divergent Specialization of Simple Venom Gene Profiles among Rear-Fanged Snake Genera (Helicops and Leptodeira, Dipsadinae, Colubridae). Toxins. 2022; 14(7):489. https://doi.org/10.3390/toxins14070489

Chicago/Turabian Style

Cerda, Peter A., Jenna M. Crowe-Riddell, Deise J. P. Gonçalves, Drew A. Larson, Thomas F. Duda, Jr., and Alison R. Davis Rabosky. 2022. "Divergent Specialization of Simple Venom Gene Profiles among Rear-Fanged Snake Genera (Helicops and Leptodeira, Dipsadinae, Colubridae)" Toxins 14, no. 7: 489. https://doi.org/10.3390/toxins14070489

APA Style

Cerda, P. A., Crowe-Riddell, J. M., Gonçalves, D. J. P., Larson, D. A., Duda, T. F., Jr., & Davis Rabosky, A. R. (2022). Divergent Specialization of Simple Venom Gene Profiles among Rear-Fanged Snake Genera (Helicops and Leptodeira, Dipsadinae, Colubridae). Toxins, 14(7), 489. https://doi.org/10.3390/toxins14070489

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