Next Article in Journal
The Deciphered Genome of Mesobuthus martensii Uncovers the Resistance Mysteries of Scorpion to Its Own Venom and Toxins at the Ion Channel Level
Next Article in Special Issue
Evolution Stings: The Origin and Diversification of Scorpion Toxin Peptide Scaffolds
Previous Article in Journal
Potential Degradation of Swainsonine by Intracellular Enzymes of Arthrobacter sp. HW08
Previous Article in Special Issue
Atractaspis aterrima Toxins: The First Insight into the Molecular Evolution of Venom in Side-Stabbers
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Three-Fingered RAVERs: Rapid Accumulation of Variations in Exposed Residues of Snake Venom Toxins

1
CIMAR/CIIMAR, Centro Interdisciplinar de Investigação Marinha e Ambiental, Universidade do Porto, Rua dos Bragas, 177, Porto 4050-123, Portugal
2
Departamento de Biologia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, Porto 4169-007, Portugal
3
Venom Evolution Lab, School of Biological Sciences, The University of Queensland, St. Lucia, Queensland 4072, Australia
4
Institute for Molecular Bioscience, The University of Queensland, St. Lucia, Queensland 4072, Australia
5
HEJ Research Institute of Chemistry, International Center for Chemical and Biological Sciences (ICCBS), University of Karachi, Karachi-75270, Pakistan
*
Author to whom correspondence should be addressed.
Submission received: 2 October 2013 / Revised: 8 November 2013 / Accepted: 11 November 2013 / Published: 18 November 2013
(This article belongs to the Collection Evolution of Venom Systems)

Abstract

:
Three-finger toxins (3FTx) represent one of the most abundantly secreted and potently toxic components of colubrid (Colubridae), elapid (Elapidae) and psammophid (Psammophiinae subfamily of the Lamprophidae) snake venom arsenal. Despite their conserved structural similarity, they perform a diversity of biological functions. Although they are theorised to undergo adaptive evolution, the underlying diversification mechanisms remain elusive. Here, we report the molecular evolution of different 3FTx functional forms and show that positively selected point mutations have driven the rapid evolution and diversification of 3FTx. These diversification events not only correlate with the evolution of advanced venom delivery systems (VDS) in Caenophidia, but in particular the explosive diversification of the clade subsequent to the evolution of a high pressure, hollow-fanged VDS in elapids, highlighting the significant role of these toxins in the evolution of advanced snakes. We show that Type I, II and III α-neurotoxins have evolved with extreme rapidity under the influence of positive selection. We also show that novel Oxyuranus/Pseudonaja Type II forms lacking the apotypic loop-2 stabilising cysteine doublet characteristic of Type II forms are not phylogenetically basal in relation to other Type IIs as previously thought, but are the result of secondary loss of these apotypic cysteines on at least three separate occasions. Not all 3FTxs have evolved rapidly: κ-neurotoxins, which form non-covalently associated heterodimers, have experienced a relatively weaker influence of diversifying selection; while cytotoxic 3FTx, with their functional sites, dispersed over 40% of the molecular surface, have been extremely constrained by negative selection. We show that the a previous theory of 3FTx molecular evolution (termed ASSET) is evolutionarily implausible and cannot account for the considerable variation observed in very short segments of 3FTx. Instead, we propose a theory of Rapid Accumulation of Variations in Exposed Residues (RAVER) to illustrate the significance of point mutations, guided by focal mutagenesis and positive selection in the evolution and diversification of 3FTx.

1. Introduction

Venoms are key evolutionary innovations in the Kingdom Animalia and are complex concoctions of biologically active proteins (from polypeptide globular enzymes to small peptides), salts, and organic molecules such as polyamines, amino acids and neurotransmitters [1]. Venom components originate via toxin recruitment events during which ordinary protein-encoding genes, typically those involved in key regulatory processes (such as hemostasis or neurotransmission) are duplicated, and the new copies are selectively expressed in the venom gland [1,2,3,4,5,6,7,8,9,10,11,12,13]. These novel paralogs can further duplicate and give rise to multigene families, following the “birth and death” mode of evolution, where the rapid evolution of these families results in extensive neofunctionalization of some copies, while the other non-functional forms are lost through degradation or get transformed into pseudogenes [14]. Research has shown that despite the extraordinary diversity of animal toxins, most belong to a limited number of enzymatic (e.g., phospholipases, serine proteases, metalloproteinases) and non-enzymatic (e.g., three-finger toxins, natriuretic peptides, Kunitz peptides, lectins) protein superfamilies which have been convergently recruited in various organisms to perform similar functions [1,3].
Three-finger toxins (3FTx) are one of the most abundantly secreted non-enzymatic components of elapid (Elapidae), colubrid (Colubridae) and psammophiide (Psammophiinae subfamily withing the Lamprophidae) snake venom. They are characterised by a broad diversity of functional forms (Table 1). In the past, 3FTx were considered to be exclusive to elapid snake venoms [6]. The discovery of α-colubritoxin, however, revealed this potent toxin type to be widespread in “non-front-fanged” (NFF) Caenophidia snake lineages [15]. Subsequent studies revealed the broad taxonomic distribution of this toxin type [16,17,18,19,20,21,22], which appears to have been recruited into the snake venom arsenal near the base of the snake tree [9].
3FTx are characterised by the presence of three β-loops that extend from the toxin’s small, hydrophobic core, giving them the “three-fingered” appearance and their name [23]. The plesiotypic (ancestral character state of a molecular scaffold) 3FTx form, such as that found in NFF advanced snakes, has ten cysteines in a distinctive pattern, reflective of its molecular origin from the recruitment of a LYNX/SLUR nicotinic receptor binding neuromodulation peptide [3,6,15,17]. All described 3FTx types from Henophidia, NFF and viperid (Viperidae) snakes contain these ten cysteines, and all of these forms that have been functionally investigated are α-neurotoxic, with a potency much greater against birds and/or reptiles than mammals [6,9,17,19,20,21,24]. This taxon specificity in action led to the plesiotypic α-neurotoxins being mistakenly referred to as “weak neurotoxins” (c.f. [25]). Similarly, this taxon-specific toxicity led to misinterpretation of prey-handling behaviour in an experiment investigating the role of venom in the feeding ecology of NFF snakes [26].
Table 1. Bioactivities of 3FTx types with characterised toxicities.
Table 1. Bioactivities of 3FTx types with characterised toxicities.
Functional ClassMode of Action
Basal-type α-neurotoxinsAntagonists of α1 nicotinic acetylcholine receptors, with a 100-fold greater potency to avians/reptiles than mammals. Produces flaccid paralysis [21]
Type I α-neurotoxinsAntagonists of α1 nicotinic acetylcholine receptors. Produces flaccid paralysis [27]
Type II α-neurotoxinsAntagonists of α1 and α7 nicotinic acetylcholine receptors. Produces flaccid paralysis [27]
Type III α-neurotoxinsAntagonists of α1 nicotinic acetylcholine receptors. Produces flaccid paralysis [28]
κ-neurotoxinsAntagonists of α3β2 neuronal nicotinic acetylcholine receptor subtype. Produces flaccid paralysis [29]
Adrenergic/Muscarinic neurotoxinsAntagonists of a wide variety of adrenergic and muscarinic subtypes with extreme specificity for receptor subtypes [30,31,32,33,34,35,36,37,38,39,40,41]
Type B Muscarinic toxinsAntagonists of M2 muscarinic acetylcholine receptors [42]
ASIC channel blockersActs as a reversible gating modifier toxin by antagonistically binding to closed/inactivated ASIC1a-ASIC2a (ACCN2-ACCN1) channels in central neurons and ASIC1b-containing channels in nociceptors [43]
Calcium channel blockersAntagonists of L-type calcium channels, thus inhibiting the transmission of the action potential [44]
Acetylcholinesterase inhibitorsInhibit acetylcholinesterase through competitive binding [45]
Platelet inhibitorsCompetitively bind to platelet GPIIb/IIIa receptor utilising the RGD functional motif, thus blocking platelet aggregation [46]
CytotoxinsCell-damaging activity mediated by hydrophobic-patch on molecular surface that interacts non-specifically with hydrophobic aspects of the cell phospholipid bilayer [47]
SynergisticAlone are non-toxic but form complexes with α-neurotoxins to dramatically enhance neurotoxicity [48]
Figure 1. Bayesian molecular phylogeny of representative three-finger toxins. Uniprot [49] accession numbers are given for each. Cysteine framework variation is displayed, with ancestral cysteines in black and newly evolved cysteines in red.
Figure 1. Bayesian molecular phylogeny of representative three-finger toxins. Uniprot [49] accession numbers are given for each. Cysteine framework variation is displayed, with ancestral cysteines in black and newly evolved cysteines in red.
Toxins 05 02172 g001
Subsequent to elapid snakes evolving a high pressure, syringe-like delivery system including venom gland compressor musculature and hollow front fangs, apotypic (derived character state of a molecular scaffold) forms of 3FTx emerged, characterised by a loss of plesiotypic cysteines 2 and 3 (Figure 1), a change that probably resulted in the dramatic potentiation of α-neurotoxicity through the uncoupling of loop-1, with these apotypic forms becoming much more potent upon mammalian receptors than the more constrained plesiotypic forms [6]. This increased toxicity likely enhanced the role of these proteins in prey capture and resulted in a high level expression of α-neurotoxins (α-ntx) lacking the second and third plesiotypic cysteines (Type I (aka: short-chain) and Type II (aka: long-chain) α-ntx) in the venom glands of these snakes. The increased level of expression was accompanied by punctuated molecular evolution, resulting in the emergence of a myriad of structurally and functionally novel forms [6]. Structurally novel forms included α-ntx with newly evolved cysteines that stabilised the second loop (Type II α-ntx). Venoms of snakes from the Oxyuranus and Pseudonaja clade are unusual in containing significant amounts of 3FTx that display all the features of Type II α-ntx (including lacking the 2nd and 3rd plesiotypic cysteines) but lack the apotypic loop-2 stabilising cysteine pair characteristic of this type [50]. These toxins would thus be expected to be phylogenetically basal to the Type II α-ntx which contain the apotypic loop-2 stabilising cysteine pair. Another derived 3FTx structural variation is represented by the κ-neurotoxins (aka: κ-bungarotoxins), which form non-covalently associated heterodimers. A number of novel functions also emerged ([6]; Table 1; Figure 1), such as κ-ntx specifically targeting the neuronal nicotinic receptors. The most extreme neofunctionalisation is represented by the cytotoxins, which have deviated from the highly focused ion channel targeting of the α-ntx [51]. Instead, they exhibit a number of novel biological activities, including lysis of various types of cells (including erythrocytes and epithelial cells), enzyme inhibition (protein kinase C: [52]; Na+/K+ ATPase: [53]), depolarization and contraction of muscle cells and prevention of platelet aggregation. All structurally and functionally apotypic 3FTxs lack plesiotypic cysteines 2 and 3 (Figure 1), and are expressed in the venoms of elapid snakes in much higher levels than the α-ntxs that contain all ten plesiotypic cysteines.
Although it has been hypothesised and demonstrated through (now obsolete) selection assessments that snake venom three-finger toxins have evolved under the influence of positive Darwinian selection [6,54,55], the underlying mechanism of evolution driving the diversification of functional forms remains unclear. The molecular evolution of several 3FTx forms, such as κ-ntx, plesiotypic 3FTxs from Henophidia, NFF and elapid snakes and Type III α-neurotoxins from Australian elapids, remain unstudied to date. Molecular mechanisms underlying the evolution of the 3FTx gene in viperid snakes, which have independently evolved a sophisticated high pressure, hollow-fanged venom delivery system (VDS), also remain elusive. Moreover, interpretations regarding the evolution of certain 3FTx types, such as the suggestion that cytotoxins evolve under the influence of positive selection (with an ω value of 19.5) [55], have been questionable. Recently, the molecular evolution of 3FTx homologues in Atractaspis sp. was evaluated and the influence of positive Darwinian selection on genes encoding them was highlighted [56].
As the general organization of the 3FTx gene is highly conserved and ordered, Accelerated Segment Switch in Exons to alter Targeting (ASSET) has been postulated as the mechanism driving the molecular evolution of 3FTx [16,57]. This theory suggests that during the evolution of the 3FTx gene, segments in exonic regions have been exchanged with distinctly different ones and the resultant “switching” of segments has generated the observed sequence variation and the functional diversity of 3FTx. The authors further speculated that point mutations alone could not account for the diversity of 3FTx functional forms and that they could only be helpful in fine tuning receptor binding capabilities, which originally arise through ASSET [16,57]. However, this study attempted to classify regions in 3FTx based on simplistic “degree of identity” comparisons. Crucially, such analysis does not take into account the fact that proteins adopt regionally differential rates of evolution [58,59,60]. It has been previously demonstrated in other toxin types that structurally important residues are constrained by negative selection, while regions responsible for biological function and/or those forming the molecular surface accumulate variations under an arms race scenario [58,59,61]. This not only facilitates functional diversification, but also increases the number of active residues on the surface of venom components that can non-specifically interact with novel receptors in prey and induce a plethora of pharmacological effects.
In the present study, we test the following hypotheses: (i) that 3FTx with specific sites of action are involved in a coevolutionary arms race with receptors of prey animals and thus are evolving under the influence of positive selection; (ii) that cytotoxic 3FTx, which interact non-specifically with cell membranes, do not experience an arms race and evolve under the constraints of negative selection pressure; (iii) that, in venoms in which 3FTx are a major component (Elapidae, Colubridae and Psammophiinae), 3FTx are rapidly evolving under positive selection (with the exception of cytotoxic 3FTx in the venom of elapid snakes); (iv) that in venoms in which 3FTx are a minor component (Viperidae), 3FTx evolve under a neutral selection regime and do not experience positive selection; and (v) the evolution of the advanced venom delivery system in elapid snakes led to an increase in diversifying pressure on 3FTx, resulting in the diversity of functional forms present in the venoms of snakes from this family. In order to test these hypotheses and provide further insight into the evolution and diversification of this toxin superfamily, we reconstructed the complex molecular evolutionary history of 3FTx. In particular, this study examined: the relative rate of evolution of the plesiotypic α-ntx (i) before and after the evolution of the advanced venom delivery apparatus in Caenophidia (advanced snakes), in order to test whether the evolution of a sophisticated VDS influenced the regime of 3FTx evolution; (ii) before and after the evolution of the sophisticated, high pressure and hollow-fanged venom delivery system in elapids and viperids to test whether the refining of VDS influenced selection pressures on 3FTx genes; (iii) the relative rate of 3FTx evolution subsequent to the loss of two plesiotypic cysteines and resultant potentiation of α-neurotoxicity; (iv) the relative rate of evolution subsequent to the apotyposis (or derivation) of non-covalently associated dimeric forms, in order to find out whether these new structural constraints affected natural selection pressures on this gene; and finally; (v) the relative rate of evolution subsequent to the apotyposis of cytotoxins was examined to understand whether the novel activity recruited by cytotoxins affected their rate of evolution.

2. Results

Bayesian and maximum-likelihood molecular phylogenetic analyses retrieved phylogenetic trees with similar topologies (Figure 1; Supplementary Figure 1), with different toxin sequences forming distinct phylogenetic clades as shown previously [3,6,17]. The Bayesian tree (convergence diagnostics: Average standard deviation of split frequencies: 0.009; PSRF = 1.00) was built using amino acid sequences, since some 3FTxs are only known from their amino acid sequences, while the maximum-likelihood tree (1,000 bootstrap replicates; GTR + I + G) was built using all the nucleotide sequences analysed in this study (n = 457). The overall topology of these phylogenetic trees was largely in concordance with the previously published 3FTx phylogeny [3,6,17].
A particularly notable finding was that the novel Type II α-ntx sequences from Oxyuranus and Pseudonaja, which lack the Type II α-ntx characteristic cysteine doublet (−2C) between plesiotypic cysteines 5 and 6, were not phylogenetically basal to the other Type II, nor were they monophyletic, but rather were nested within the regular Type II forms (+2C) (Figure 2). This suggests that the apotypic cysteine doublet characteristic of Type II α-ntx was secondarily lost on at least three occasions.
Figure 2. Bayesian molecular phylogeny and structural and functional evolution of +2C/−2C Type II (long-chain) α-neurotoxins. Pseudonaja/Oxyuranus −2C and +2C sequences are coloured purple and green, respectively. Sequences presented (uniprot): (1) R4FIT5 Pseudonaja modesta, (2) R4FIU6 Pseudonaja modesta, (3) A8HDK6 Pseudonaja textilis, (4) R4FK68 Pseudonaja modesta, (5) A8HDK8 Oxyuranus microlepidotus, (6) A7X4Q3 Oxyuranus microlepidotus, (7) A8HDK7 Oxyuranus microlepidotus, (8) A7X4R0 Oxyuranus microlepidotus, (9) A8HDK9 Oxyuranus scutellatus, (10) Q9W7J5 Pseudonaja textilis, (11) R4FIT0 Pseudonaja modesta, (12) R4G7K3 Pseudonaja modesta, (13) R4G321 Pseudonaja modesta, (14) R4G2J4 Pseudonaja modesta, (15) R4G319 Pseudonaja modesta and (16) R4FK75 Pseudonaja modesta.
Figure 2. Bayesian molecular phylogeny and structural and functional evolution of +2C/−2C Type II (long-chain) α-neurotoxins. Pseudonaja/Oxyuranus −2C and +2C sequences are coloured purple and green, respectively. Sequences presented (uniprot): (1) R4FIT5 Pseudonaja modesta, (2) R4FIU6 Pseudonaja modesta, (3) A8HDK6 Pseudonaja textilis, (4) R4FK68 Pseudonaja modesta, (5) A8HDK8 Oxyuranus microlepidotus, (6) A7X4Q3 Oxyuranus microlepidotus, (7) A8HDK7 Oxyuranus microlepidotus, (8) A7X4R0 Oxyuranus microlepidotus, (9) A8HDK9 Oxyuranus scutellatus, (10) Q9W7J5 Pseudonaja textilis, (11) R4FIT0 Pseudonaja modesta, (12) R4G7K3 Pseudonaja modesta, (13) R4G321 Pseudonaja modesta, (14) R4G2J4 Pseudonaja modesta, (15) R4G319 Pseudonaja modesta and (16) R4FK75 Pseudonaja modesta.
Toxins 05 02172 g002
Our analysis of the venom gland transcriptome of the unique viperid species Azemiops feae recovered a 3FTx transcript that was found in the same clade as the representative sequences from the pit-viper Sistrurus catenatus (Figure 1; Supplementary Figure 1). This is the first non-Sistrurus 3FTx to be recovered from a viperid snake. Based on an unresolved neighbour joining tree, viperid 3FTx were previously reported as polyphyletic, with some sequences nested within elapid 3FTx clade, and others in the “non-front-fanged” advanced snake 3FTx clade [24]. However, both Bayesian and maximum-likelihood phylogenetic analyses in this study place all viperid sequences in a monophyletic clade outside elapid 3FTx clade with strong node support (Figure 1; Supplementary Figure 1). This is consistent with the viperid snakes having diverged from the remaining advanced snakes nearly 54 million years ago [62]. viperid 3FTx homologs retrieved to date indicate that apotypic structural and functional forms, such as Type I, II and III α-ntxs, κ-ntxs and cytotoxins, are the result of apotyposis in elapid snakes after this split (Figure 1; Supplementary Figure 1).
The one-ratio model (ORM), the simplest of the lineage-specific models, computed a wide range of ω (dN/dS) values for various types of 3FTx (Supplementary Tables 1–10). This highly conservative model can only detect positive selection when the ω ratio, averaged over all sites along the lineages in a phylogenetic tree, is significantly greater than one. Despite this, the ω value for most three-finger toxins was significantly greater than 1, indicating the strong influence of positive selection in shaping the evolution of 3FTx: plesiotypic α-ntx from viperids, NFF and elapids: 1.79, 1.29 and 1.30, respectively; Type I, Type II and Type III α-ntx: 1.92, 2.01 and 2.59, respectively; and κ-ntx: 1.64. Oxyuranus/Pseudonaja Type II α-ntx with (+2C) and without the cys doublet (−2C) were significantly different, with ω values of 0.97 and 2.69, respectively. In contrast to all other 3FTx types, ORM estimated an ω value of 0.32 for cytotoxins, indicating an unprecedented lack of variation in this unique three-finger toxin.
The segregation of seromucous mixed maxillary glands into discrete protein and mucus glands at the base of the advanced snake tree and the subsequent evolution of a high pressure, hollow-fanged venom delivery apparatus independently in Atractaspis/Homoroselaps, elapids and viperids were major evolutionary advances in snake-feeding ecology [7,17]. In order to assess the effect of apomorphic (derived state of a morphological character) venom delivery systems on the evolution of 3FTx, we employed the three-ratio model (3RM) and the branch-site test A (BST) on the plesiotypic α-ntx from the advanced snakes (NFF, viperids and elapids). While 3RM—which assesses selection pressure only across lineages—estimated ω of 1.13, 1.37 and 2.21 for the plesiotypic α-ntx from NFF, elapids and viperids, respectively; BST—which assess selection pressures across the sites and along the lineages in a tree—estimated ω of 2.16 (22% sites), 4.09 (23% sites) and 6.54 (31% sites), respectively (Table 2). NFF and elapid comparisons were insignificant (p > 0.05) for 3RM, while all other comparisons were significant (significant at 0.001 even after Bonferroni corrections). 3RM and BST indicated a greater influence of positive selection on the plesiotypic α-ntxs in elapid and viperid snakes, which have independently evolved sophisticated high pressure and hollow-fanged venom delivery systems.
Table 2. Lineage-specific analyses of plesiotypic Caenophidia three-finger toxins (3FTxs).
Table 2. Lineage-specific analyses of plesiotypic Caenophidia three-finger toxins (3FTxs).
Modelω aLikelihood (ι)Prop. of Sites with ω > 1 bSignificance c
Plesiotypic α-neurotoxins from ‘non-front-fanged’ advanced snakes
Three-ratio model1.13−5982.232579-P > 0.05 NS
Branch-site model A2.16−5767.71045422.1%* P << 0.001
Clade model C2.20−5670.14924649.9%* P << 0.001
Plesiotypic α-neurotoxins from Elapidae
Three-ratio model1.37−5984.727017-P > 0.05NS
Branch-site model A4.09−5773.84670722.6%* P << 0.001
Clade model C4.48−5670.14924649.9%* P << 0.001
Plesiotypic α-neurotoxins from Viperidae
Three-ratio model2.21−5981.804559-* P << 0.001
Branch-site model A6.54−5762.85319431.1%* P << 0.001
Clade model C6.23−5670.14924649.9%* P << 0.001
Legend: a: mean dN/dS; b: Proportion of sites with ω > 1; c: Significance of the model in comparison with the null model; * : Bonferroni corrected values; NS: Not significant (at 0.05).
We also employed these tests to assess selection pressures along the lineages of the three types of elapid-specific apotypic α-ntx that lack plesiotypic cysteines 2 and 3: Type I, II and III α-ntxs (Table 3). The three-ratio model (3RM) estimated ω of 1.89, 1.94 and 2.39 for Type I, Type II and Type III α-ntxs, respectively. BST estimated ω of 4.39 (10% sites), 4.09 (18% sites) and 4.77 (24% sites), respectively for the three types. All these comparisons were significant at 0.001 even after Bonferroni corrections.
Table 3. Lineage-specificanalyses of α-neurotoxins.
Table 3. Lineage-specificanalyses of α-neurotoxins.
Modelω aLikelihood (ι)Prop. of Sites with ω > 1 bSignificance c
Type I α-neurotoxins
Three-ratio model1.89−12,395.465265-* P << 0.001
Branch-site model A4.39−11,868.1980459.8%* P << 0.001
Clade model C3.36−11,602.54705044.4%* P << 0.001
Type II α-neurotoxins
Three-ratio model1.94−12,395.367372-* P << 0.001
Branch-site model A4.09−11,751.04562118.3%* P << 0.001
Clade model C4.14−11,602.54705044.4%* P << 0.001
Type III α-neurotoxins
Three-ratio model2.39−12,394.644636-* P << 0.001
Branch-site model A4.77−11,812.91289023.7* P << 0.001
Clade model C3.26−11,602.54705044.4%* P << 0.001
Legend: a: mean dN/dS; b: Proportion of sites with ω > 1; c: Significance of the model in comparison with the null model; *: Bonferroni corrected values; NS: Not significant (at 0.05).
Since the aforementioned lineage-specific models often fail to detect episodic diversifying selection that only affects certain sites in proteins, we also employed site-specific models. Similar to ORM, model 8 estimated ω that ranged widely (Table 4; Figure 3 and Figure 4; Supplementary Tables 1–10). The computed ω and the number of positively selected sites detected by model 8 suggest that most three-finger toxins have evolved under the strong influence of positive selection; viperids: 3.28 (30 positively selected (PS) sites; 55% of sites)); plesiotypic α-ntx from NFF: 1.63 (37 PS; 47% of sites); plesiotypic α-ntx from elapid snakes: 1.75 (28 PS; 41% of sites); Type I α-ntx: 1.72 (19 PS; 30.0% of sites); Type II α-ntx: 1.45 (21 PS; 26% of sites); Oxyuranus/Pseudonaja Type II α-ntx with cys doublet (+2C): 3.67 (18 PS; 33% of sites); Oxyuranus/Pseudonaja Type II α-ntx without cys doublet (−2C): 3.41 (20 PS; 58% of sites); Type III α-ntx 2.61 (30 PS; 42% of sites); κ-bungarotoxins: 2.11 (5 PS; 9.3% of sites) (Table 4; Figure 3 and Figure 4; Supplementary Tables 1–10). FUBAR and MEME identified several episodically diversifying sites in various forms of 3FTxs, highlighting the rapid evolution of these toxins (Table 4). In contrast, M8 estimated ω of 0.53 (2 PS and 3% of sites) for cytotoxic 3FTx, suggesting that they have evolved under a regime of strong negative selection. Compared to the α-ntxs, κ-bungarotoxins have experienced a weaker influence of positive diversifying selection, as reflected by the very small number of positively selected sites detected in this lineage by various analyses (5PS and 9% of sites; Table 4).
Table 4. Dynamic molecularevolution of three-finger toxins (3FTx).
Table 4. Dynamic molecularevolution of three-finger toxins (3FTx).
FUBAR aMEME bBSR cPAML d
#sites#branchesM8M2a
Plesiotypic α-neurotoxins
Henophidia 3FTxω > 1e : 0
ω < 1f : 0
0000
ω: 0.78ω: 0.78
Non-front-fanged advanced snake 3FTxω > 1e : 18
ω < 1f : 7
1893726
(27 + 10)(18 + 8)
ω: 1.63ω: 1.71
Elapidae 3FTxω > 1e : 13
ω < 1f : 3
682824
(22 + 6)(15 + 9)
ω: 1.75ω: 1.79
Viperidae 3FTxω > 1e : 7
ω < 1f : 1
063024
(18 + 12)(13 + 11)
ω: 3.28ω: 3.28
Apotypic or Derived α-neurotoxins
Type I α-neurotoxinω > 1e : 22
ω < 1f : 6
1391913
(13 + 6)(11 + 2)
ω: 1.72ω: 1.61
Type II α-neurotoxinω > 1e : 21
ω < 1f : 11
30242119
(19 + 2)(17 + 2)
ω: 1.45ω: 1.43
Type III α-neurotoxinω > 1e : 29
ω < 1f : 4
24133027
(26 + 4)(26 + 1)
ω: 2.61ω: 2.59
Oxyuranus/Pseudonaja derived Type II α-neurotoxins
With-Cysteine- Doubletω > 1e : 11
ω < 1f : 2
671817
(11 + 7)(10 + 7)
ω: 3.67ω: 3.57
Without-Cysteine-Doubletω > 1e : 5
ω < 1f : 1
122014
(9 + 11)(2 + 12)
ω: 3.41ω: 3.41
Cytotoxinsω > 1e : 1
ω < 1f : 2
7022
(2 + 0)(1 + 1)
ω: 0.53ω: 0.57
Kappa neurotoxinsω > 1e : 6
ω < 1f : 0
3255
(3 + 2)(3 + 2)
ω: 2.11ω: 2.11
Legend: a: Fast Unconstrained Bayesian AppRoximation; b: Sites detected as experiencing episodic diversifying selection (0.05 significance) by the Mixed Effects Model Evolution (MEME); c: Number of branches detected by the branch-site REL (Random effects likelihood) test as episodically diversifying (not comparable to positively selected sites); d: Positively selected sites detected by the Bayes Empirical Bayes approach implemented in M8 and M2a. Sites detected at 0.99 and 0.95 significance are indicated in the parenthesis (total number of positively selected sites is also indicated along with the estimated ω); e: number of sites under pervasive diversifying selection at the posterior probability ≥0.9 (FUBAR); f: Number of sites under pervasive purifying selection at the posterior probability ≥0.9 (FUBAR); ω: mean dN/dS.
Figure 3. Molecular evolution of plesiotypic 3FTxs. Three-dimensional homology models of various three-finger toxins (top: surface-fill; bottom: wireframe), depicting the locations of positively selected sites are presented. Site-model 8 computed omega and the total number of positively selected sites detected by its Bayes Empirical Bayes (BEB) approach (PP ≥ 0.95) are also indicated.
Figure 3. Molecular evolution of plesiotypic 3FTxs. Three-dimensional homology models of various three-finger toxins (top: surface-fill; bottom: wireframe), depicting the locations of positively selected sites are presented. Site-model 8 computed omega and the total number of positively selected sites detected by its Bayes Empirical Bayes (BEB) approach (PP ≥ 0.95) are also indicated.
Toxins 05 02172 g003
Figure 4. Molecular evolution of derived 3FTx. Three-dimensional homology models of various three-finger toxins (top: surface-fill; bottom: wireframe), depicting the locations of positively selected sites are presented. Site-model 8 computed omega and the total number of positively selected sites detected by its Bayes Empirical Bayes (BEB) approach (PP ≥ 0.95) are also indicated.
Figure 4. Molecular evolution of derived 3FTx. Three-dimensional homology models of various three-finger toxins (top: surface-fill; bottom: wireframe), depicting the locations of positively selected sites are presented. Site-model 8 computed omega and the total number of positively selected sites detected by its Bayes Empirical Bayes (BEB) approach (PP ≥ 0.95) are also indicated.
Toxins 05 02172 g004
Since lineages may have different proportions of sites under the influence of selection, the direct comparison of ω values computed by the aforementioned analyses could be misleading. Hence, we employed the clade model c approach to simultaneously compute ω values for different 3FTx clades. Clade-model analyses computed ω of 2.20, 4.48 and 6.23 for the plesiotypic α-ntxs from NFF, elapid and viperid snakes, respectively (Table 2). This highlights a greater influence of positive selection on 3FTx from elapid and viperid snake lineages (significant at 0.001 in comparison with site model M1a, even after Bonferroni correction), relative to NFF and Henophidia 3FTxs. Within elapid snake specific 3FTx lineages, clade analyses estimated ω of 3.36, 4.14 and 3.26 for Type I, Type II and Type III α-ntxs, respectively, indicating a greater influence of positive Darwinian selection on Type II α-ntx compared to the rest (Table 3).
To provide further support to the positively detected sites by the nucleotide-level selection analyses (M8), we employed a complementary protein-level approach implemented in TreeSAAP. With this combined approach we were able to identify 18/19 positively selected (PS) sites in Type I, all 21 PS sites in Type II, 13/30 PS sites in Type III α-ntxs; 21/30 sites in viperid 3FTx; 11/37 PS sites in NFF and 2/2 PS sites in cytotoxic 3FTx (Supplementary Table 11: highlighted in bold). However, this approach failed to identify positively selected sites in κ-bungarotoxins and elapid plesiotypic 3FTx lineages and +2C and −2C type II α-ntx lineages, suggesting that mutations in these toxin types may not be as radical as the ones in the aforementioned types.
Evolutionary fingerprint analyses clearly depict several residues in most 3FTx as evolving under the influence of positive selection, while conversely the majority of sites in cytotoxic 3FTx remain constrained under negative selection (Supplementary Figure 2). Branch-site tests assume that a phylogenetic tree can be divided into foreground (branches experiencing positive selection) and background (branches evolving under neutral or negative selection) lineages. It has been demonstrated that these tests can suffer from high rates of false-positives and false-negatives when such modeling assumptions are violated [63]. Hence, we employed the branch-site REL (BSR) test, which does not make such assumptions and outperforms branch-site tests. This test identified a number of branches in the phylogenetic trees of various 3FTx types as episodically diversifying: plesiotypic α-ntx from NFF, elapid and viperid snakes: 9, 8 and 6 branches, respectively; Type I, Type II and Type III α-ntx 9, 24 and 13 branches, respectively; +2C and −2C Oxyuranus/Pseudonaja Type II α-ntx: 7 and 2 branches, respectively (Supplementary Figures 3–5). This test detected two and three branches, respectively in κ-ntxs and cytotoxins, indicating that even these toxin types have experienced the influence of episodic bursts of selection, probably during their early stages of evolution (Supplementary Figure 6). Computation of accessible surface area ratio indicated that a majority of positively selected sites detected by M8 in 3FTx were exposed (completely or partially) to the surrounding medium (Figure 5; Supplementary Tables 11). The likelihood of a hypermutable site to be exposed to the surrounding medium was 2 to 3.5 times greater, compared to the likelihood of it being buried (Table 5). Mapping of positively selected sites onto the alignment of various 3FTx genes revealed that nearly all the hypermutable sites are located outside the structurally and functionally important regions, suggesting that 3FTx genes are subjected to a very sophisticated focal mutagenesis (Figure 2, Figure 6, Figure 7 and Figure 8; Supplementary Table 12).
Figure 5. Surface accessibility of three-finger toxins. A plot of amino acid positions (x-axis) against accessible surface area (ASA) ratio (y-axis) indicating the locations of amino acids (exposed or buried) in the crystal structure of various three-finger toxins is presented. Positively selected residues are presented as large dots, while the remaining sites are presented as small dots in the plot. Residues with an ASA ratio greater than 50% are considered to be exposed (ASA of 40%–50% are likely to have exposed side chains), while those with a ratio lesser than 20% are considered to be buried to the surrounding medium. Three-dimensional structures of various 3FTx types, depicting the locations of positively selected (PS) sites along with model 8 omega values and the number of exposed and buried positively selected sites are also presented (for details see Table 5).
Figure 5. Surface accessibility of three-finger toxins. A plot of amino acid positions (x-axis) against accessible surface area (ASA) ratio (y-axis) indicating the locations of amino acids (exposed or buried) in the crystal structure of various three-finger toxins is presented. Positively selected residues are presented as large dots, while the remaining sites are presented as small dots in the plot. Residues with an ASA ratio greater than 50% are considered to be exposed (ASA of 40%–50% are likely to have exposed side chains), while those with a ratio lesser than 20% are considered to be buried to the surrounding medium. Three-dimensional structures of various 3FTx types, depicting the locations of positively selected (PS) sites along with model 8 omega values and the number of exposed and buried positively selected sites are also presented (for details see Table 5).
Toxins 05 02172 g005
Table 5. Relative frequency of positively selected sites.
Table 5. Relative frequency of positively selected sites.
Toxin TypeTotal ResiduesExposed Total *Exposed PSExposed PS %Buried Total *Buried PSBuried PS %Ratio Exposed % PS:Buried % PS
Plesiotypic NFF74402562.517317.63.5
Plesiotyic Viperidae72402255%204202.8
Plesiotypic Elapidae6539235914321.42.8
Type I α-ntx62381539.5102202
Type II α-ntx74381436.8203152.5
Type III α-ntx57342058.814321.42.7
κ-ntx643538.5170-N.A.
Cytotoxic 3FTx603425.9150-N.A.
PS = positively selected; N.A. = not applicable; NFF: ‘non-front-fanged’ advanced snakes; Ntx: Neurotoxin; * Some residues were ambiguous as being exposed or buried and thus were not assigned to either category.
Figure 6. Structural and functional evolution of plesiotypic 3FTxs. Sequence alignment and homology models depicting structurally and functionally important residues and hypermutable sites of various three-finger toxins are shown. Extremely well-conserved residues implicated in structural/functional roles, and hypermutable sites are shaded. Sequences presented (uniprot): A. (1) M9T1L2 Aspidites melanocephalus, (2) M9T271 Aspidites melanocephalus, (3) M9SZR1 Cylindrophis ruffus and (4) M9SZV9 Cylindrophis ruffus; B. (1) A5X2W8 Sistrurus catenatus edwardsii, (2) B4Y146 Sistrurus catenatus edwardsii, (3) A5X2W7 Sistrurus catenatus edwardsii, (4) B4Y143 Sistrurus catenatus edwardsii, (5) B4Y144 Sistrurus catenatus edwardsii and (6) M9T2J4 Azemiops feae; C. (1) Q06ZW0 Boiga dendrophila, (2) A0S864 Boiga irregularis, (3) A0S865 Boiga irregularis, 4) A7X3V0 Telescopus dhara and (5) A7X3S0 Trimorphodon biscutatus; D. (1) Q6IZ95 Bungarus candidus, (2) Q9PW19 Bungarus multicinctus, (3) Q9YGH9 Bungarus multicinctus, (4) Q8AY51 Bungarus candidus, (5) Q2VBN2 Ophiophagus hannah and 6) Q9YGI2 Naja atra.
Figure 6. Structural and functional evolution of plesiotypic 3FTxs. Sequence alignment and homology models depicting structurally and functionally important residues and hypermutable sites of various three-finger toxins are shown. Extremely well-conserved residues implicated in structural/functional roles, and hypermutable sites are shaded. Sequences presented (uniprot): A. (1) M9T1L2 Aspidites melanocephalus, (2) M9T271 Aspidites melanocephalus, (3) M9SZR1 Cylindrophis ruffus and (4) M9SZV9 Cylindrophis ruffus; B. (1) A5X2W8 Sistrurus catenatus edwardsii, (2) B4Y146 Sistrurus catenatus edwardsii, (3) A5X2W7 Sistrurus catenatus edwardsii, (4) B4Y143 Sistrurus catenatus edwardsii, (5) B4Y144 Sistrurus catenatus edwardsii and (6) M9T2J4 Azemiops feae; C. (1) Q06ZW0 Boiga dendrophila, (2) A0S864 Boiga irregularis, (3) A0S865 Boiga irregularis, 4) A7X3V0 Telescopus dhara and (5) A7X3S0 Trimorphodon biscutatus; D. (1) Q6IZ95 Bungarus candidus, (2) Q9PW19 Bungarus multicinctus, (3) Q9YGH9 Bungarus multicinctus, (4) Q8AY51 Bungarus candidus, (5) Q2VBN2 Ophiophagus hannah and 6) Q9YGI2 Naja atra.
Toxins 05 02172 g006
Figure 7. Structural and functional evolution of derived α-neurotoxins and κ-neurotoxins. Sequence alignment and homology models depicting structurally and functionally important residues and hypermutable sites of various three-finger toxins are shown. Extremely well-conserved residues implicated in structural/functional roles, and hypermutable sites are shaded. Sites that are shown to be structurally/functionally important but exhibit less than 70% identity are indicated by arrow heads in the alignment. Sequences presented (uniprot): (A). (1) B5G6F6 Oxyuranus microlepidotus, (2) H8PG58 Suta nigriceps, (3) P60770 Naja atra, (4) P10455 Laticauda colubrina, (5) B2BRS2 Austrelaps labialis and (6) F8J2H2 Drysdalia coronoides; (B). (1) F8J2E0 Drysdalia coronoides, (2) Q8UW29 Hydrophis hardwickii, (3) P01384 Notechis scutatus, (4) P82662 Ophiophagus hannah, (5) R4G2L3 Suta fasciata and (6) R4G2E5 Brachyurophis roperi; C. (1) R4G2S0 Cacophis squamulosus, (2) R4G7H6 Furina ornata, (3) R4G332 Pseudonaja modesta, (4) Q9W7K2 Pseudonaja textilis, (5) R4G7F3 Brachyurophis roperi and (6) R4G7M0 Vermicella annulata; D. (1) Q8AY56 Bungarus candidus, (2) Q8AY55 Bungarus candidus, (3) P15817 Bungarus multicinctus, (4) P01398 Bungarus multicinctus, (5) O12962 Bungarus multicinctus and (6) Q9W729 Bungarus multicinctus.
Figure 7. Structural and functional evolution of derived α-neurotoxins and κ-neurotoxins. Sequence alignment and homology models depicting structurally and functionally important residues and hypermutable sites of various three-finger toxins are shown. Extremely well-conserved residues implicated in structural/functional roles, and hypermutable sites are shaded. Sites that are shown to be structurally/functionally important but exhibit less than 70% identity are indicated by arrow heads in the alignment. Sequences presented (uniprot): (A). (1) B5G6F6 Oxyuranus microlepidotus, (2) H8PG58 Suta nigriceps, (3) P60770 Naja atra, (4) P10455 Laticauda colubrina, (5) B2BRS2 Austrelaps labialis and (6) F8J2H2 Drysdalia coronoides; (B). (1) F8J2E0 Drysdalia coronoides, (2) Q8UW29 Hydrophis hardwickii, (3) P01384 Notechis scutatus, (4) P82662 Ophiophagus hannah, (5) R4G2L3 Suta fasciata and (6) R4G2E5 Brachyurophis roperi; C. (1) R4G2S0 Cacophis squamulosus, (2) R4G7H6 Furina ornata, (3) R4G332 Pseudonaja modesta, (4) Q9W7K2 Pseudonaja textilis, (5) R4G7F3 Brachyurophis roperi and (6) R4G7M0 Vermicella annulata; D. (1) Q8AY56 Bungarus candidus, (2) Q8AY55 Bungarus candidus, (3) P15817 Bungarus multicinctus, (4) P01398 Bungarus multicinctus, (5) O12962 Bungarus multicinctus and (6) Q9W729 Bungarus multicinctus.
Toxins 05 02172 g007
Figure 8. Structural and functional evolution of cytotoxins. Sequence alignment and homology models depicting A. the locations of putative functional residues, and B. hydrophobic regions of cytotoxins are presented. Extremely well-conserved residues implicated in structural/functional roles, and hypermutable sites are shaded. Sequences presented (uniprot): (1) P60301 Naja atra, (2) P60303 Naja kaouthia, (3) Q9DGH9 Naja kaouthia, (4) P60301 Naja sputatrix, (5) Q02454 Naja sputatrix and (6) O93471 Naja sputatrix.
Figure 8. Structural and functional evolution of cytotoxins. Sequence alignment and homology models depicting A. the locations of putative functional residues, and B. hydrophobic regions of cytotoxins are presented. Extremely well-conserved residues implicated in structural/functional roles, and hypermutable sites are shaded. Sequences presented (uniprot): (1) P60301 Naja atra, (2) P60303 Naja kaouthia, (3) Q9DGH9 Naja kaouthia, (4) P60301 Naja sputatrix, (5) Q02454 Naja sputatrix and (6) O93471 Naja sputatrix.
Toxins 05 02172 g008

3. Discussion

3.1. Evolution of an Advanced Venom Delivery Apparatus and Subsequent Evolution of a High Pressure, Hollow-Fanged Venom Delivery System May Have Facilitated Rapid Evolution of Three-Finger Toxins.

The early evolution of a venom system has been shown to have played a significant role in the diversification of Toxicofera reptiles [9,17,64,65,66,67] and the morphology and development of the venom delivery system has been shown to be intimately linked to feeding ecology/prey-handling behaviour [68,69,70,71,72,73,74,75,76,77]. This relationship between venom system morphology and feeding ecology has been discussed in detail elsewhere [7,78,79,80,81]. While some Henophidia do not use constriction for prey capture (e.g., Cylindrophis) and have mandibular and maxillary venom glands that are predominantly protein secreting, others (such as Boidae and Pythonidae) are not dependent on envenoming, and instead use an apomorphic strategy (constriction) for prey subjugation [9]. In contrast to Henophidia snakes, elapid and viperid snakes possess sophisticated venom delivery apparatus that facilitate efficient injection of venom into the prey. Various selection analyses employed in this study indicated that these contrasting prey capture strategies and the evolution of the advanced venom-delivery apparatus in Caenophidia have significantly influenced the evolution of the plesiotypic α-ntxs (Figure 3 and Figure 6; Table 2 and Table 4). Not surprisingly, plesiotypic α-ntxs in NFF (ω = 1.63 and 7 Positively Selected (PS) sites), viperid (ω = 3.28 and 30 PS) and elapid (ω = 1.75 and 28 PS) accumulated tremendous amount of variations, while their Henophidia homologs (ω = 0.77 and 0 PS) lacked variation completely (Table 2 and Table 4; Figure 3 and Figure 6). This is consistent with 3FTx variation being unnecessary in powerfully constricting snakes, which are not known to use venom for capturing prey [9]. Interestingly, the Henophidiae genera Aspidites (Pythonidae) and Cylindrophis (Cylindrophiidae) are separated by ~92 million years of evolution, almost as long as their common ancestor has been separated from the advanced snakes (approximate date of divergence 103 mya: [62]). Extreme conservation of 3FTx in the largely mucus-secreting maxillary and mandibular glands of these basally divergent snakes [9], despite their 90+ million years of independent evolution, is noteworthy. Due to the paucity of 3FTx sequences (primarily a result of their extremely low levels of venom gland specific transcription), the precise role of these proteins in the venom glands of Henophidia and the reason for the lack of variation in their coding sequences remains unclear.
The significant amount of variation observed in the plesiotypic α-ntxs from the venoms of NFF is consistent with high expression levels in the Colubridae and Psammophiinae [17,18]. 3FTx from several NFF snakes have been demonstrated to be potently toxic, and in some cases exhibit prey-specific toxicity [6,20,21,82]. However, various selection analyses presented in this study indicate a greater accumulation of variations in the plesiotypic α-ntxs from viperids and elapids, relative to the NFF snake α-ntx homologs: Table 2 and Table 4; Figure 3 and Figure 6; Supplementary Tables 1–3; Supplementary Figures 2 and 3). The sophisticated venom delivery apparatus of elapids facilitates efficient injection of the toxin cocktail into prey animals. Thus, the evolution of the high pressure, hollow-fanged venom delivery apparatus appears to have had a considerable influence on the rapid diversification of elapid and viperid 3FTxs in comparison with NFF snakes. The polyphyletic NFF snakes have been demonstrated to employ a characteristically different venom arsenal from that of elapid snakes. The venom in most elapid snakes is rich in three-finger toxins, while only the venoms of gracile snakes in Colubridae or Psammophiinae have been shown to contain large amounts of this toxin type [6,17,18]. Hence, the differences in overall venom composition could have also influenced the differential rate of 3FTx evolution in these lineages. It has already been demonstrated that some venom components that are secreted in large amounts by NFF snakes experience a greater influence of positive selection relative to those of their elapid snake counterparts [58]. Since most viperid snakes are known to employ a haemorrhagic, haemolytic toxin cocktail, the neurotoxic 3FTx may only play an ancillary role in envenoming by these snakes. Evidently, viperid snake venoms are dominated by hydrolytic enzymes (such as snake venom metalloproteinases and serine proteases), while 3FTx expression appears to occur only in insignificant amounts. Hence, the presence of considerable coding sequence variation in the viperid 3FTx gene is puzzling. Further sequencing and functional assessment of 3FTx from this lineage is warranted to precisely understand the participation of these proteins in envenoming.

3.2. The Deletion of Plesiotypic Cysteine Resides (2 and 3) May Have Triggered an Explosive Diversification

Comparative assessments revealed that although all three forms of apotypic α-ntxs that lack the 2nd and 3rd plesiotypic cysteines have evolved particularly rapidly, Type II forms have experienced a greater influence of positive selection than Type I or Type III α-ntxs (Table 3 and Table 4; Figure 4 and Figure 7; Supplementary Tables 4–6; Supplementary Figure 4A,B,C). Unlike Type I α-ntx that can only antagonise the muscle nicotinic acetylcholine receptors (α1-nAChR) (the plesiotypic functional activity), Type II α-ntxs can strongly bind to neuronal α7-nACh receptors, as well. This expanded target potential makes them extremely useful for prey paralysis and probably explains the greater influence of positive selection on this type of 3FTx. Type III α-ntx has been demonstrated to exhibit activities similar to that of Type I α-ntxs, but at a lower potency, in experiments conducted on mice [28]. Prey-specific toxicity has been demonstrated in a variety of α-ntxs [20,21,82]. Hence, Type III α-ntx in the venoms of Australian elapids could be less potent against mammals, but extremely potent against the organisms on which they specialise, such as lizards [83,84,85,86,87]. The majority of Type III α-ntxs have been recovered from venom gland transcriptomes of small Australian elapid snakes that feed exclusively on reptiles. The possibility that this toxin type may exhibit taxon-specific toxicity should be investigated further, as it may account for the tremendous influence of positive selection on this toxin type.

3.3. The Dramatic Rate of Evolution of Type II α-Neurotoxins May Have Resulted in the Loss of Characteristic Cysteine Doublets on Multiple Occasions in Oxyuranus/Pseudonaja

Oxyuranus and Pseudonaja, two closely related genera of Australian elapid snakes with unusually potent venom, are known to express the regular Type II α-ntx (+2C) and an additional type that lacks the two apotypic cysteine residues involved in stabilising loop 2 (−2C). It was previously suggested that the absence of the cysteine pair (−2C) was the phylogenetically basal form and that the Type II α-ntx form with the cysteine pairs (+2) represented the more derived form [88]. However, our molecular phylogenetic analyses revealed that the −2C originate from within the larger Type II α-ntx clade and thus represent secondary losses of the apotypic loop-2 stabilising cysteine pair (Figure 1 and Figure 2; Supplementary Figure 1). Type II α-ntxs have experienced a strong influence of positive selection pressure and some of these in Oxyuranus and Pseudonaja species have accumulated mutations even in structurally important regions that are typically constrained by negative-selection. As a result, on at least three occasions, they have lost the two apotypic cysteine residues that are characteristic of Type II α-ntxs (Table 3 and Table 4; Figure 2 and Figure 4; Supplementary Figures 4B and 5). Taxon-specific variations in potency of the +2C versus −2C have not been investigated but such studies would be illuminating.

3.4. Not All Three-Finger Toxin Clades Evolve via Pervasive Positive Darwinian Selection

Evidence provided by various analyses [Codeml site, branch, branch-site and clade-specific models: M8 and M2a, three-ratio model, branch-site test A, clade model c; HyPhy: MEME, FUBAR, Branch-site REL; TreeSAAP: protein-level selection analyses and evolutionary fingerprint analyses (Table 2, Table 3, and Table 4; Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, and Figure 7; Supplementary Tables 1–9; Supplementary Figures 2–6)], clearly demonstrates the substantial influence of positive diversifying selection on the evolution of various three-finger toxin types. Conspicuously, while most 3FTx lineages evolved rapidly, cytotoxins (ω = 0.53; 2 PS-sites) remain constrained by negative selection (Table 4; Figure 4 and Figure 8; Supplementary Figures 2 and 6; Supplementary Table 10). The precise molecular mechanism by which cytotoxic 3FTx exert their cell-damaging effects remains unclear, but it has been proposed that a highly conserved hydrophobic patch on the cytotoxin’s molecular surface interacts with hydrophobic regions of the cell’s phospholipid bilayer [89,90]. Hydrophobic residues that facilitate such non-specific interactions are spread over the three loops of the cytotoxins (Figure 8: Panel B) and cover nearly 40% of the molecular surface. Additionally, cytotoxins are known to exhibit a diversity of biological functions and the residues implicated in these activities are distributed all along the length of the toxin sequence (Figure 8; Supplementary Table 10). In contrast to α-ntx, which bind to specific receptor sites, cytotoxins bind in a much more non-specific manner and thus may escape the classic predator–prey chemical arms race that influences the diversification of other toxin types. Thus, a lack of variation within the cytotoxins is reflective of the tremendous negative selection pressure on their structurally and functionally important residues, which further prevents the loss of activity and/or potency. In contrast, α-neurotoxins and an array of the ion channels they target, likely exert reciprocal selection pressure upon each other which results in the diversification of the 3FTx gene, in a classic “Red Queen” coevolutionary arms race scenario.
Conversely, κ-ntxs (ω = 2.11; 5 PS-sites) have accumulated fewer variations than the α-ntxs (Table 4; Figure 4 and Figure 7: Panel D; Supplementary Figures 2 and 6), despite also targeting a narrowly defined ion-channel type [91]. The lack of variations in the κ-ntxs could be attributed to their unique ability amongst 3FTx to form non-covalent dimers [91]. Thus, most residues in these toxins evolve under the influence of negative selection, preventing the accumulation of mutations that might affect the three-dimensional structure of the toxin, and thus the ability to non-covalently form dimers. It is noteworthy that κ-ntxs are characterised by two stabilising cysteines in loop 2, which maintain a more rigid three dimensional structure and are homologous to the same apotypic disulfide-bond in loop II of the Type II α-ntxs. Besides the van der Waal’s bond between Phe47 and Leu55 (Figure 7: Panel D), as many as six main chain and three side chain hydrogen bonds are known to stabilise this toxin [92]. The necessity to form these bonds exerts a strong structural constraint against the accumulation of variation in κ-ntxs. Thus variation in κ-ntx is largely limited to five positively selected sites, most of which are found in loop II. Loop III and the region between C9–C10, along with the functionally/structurally important residues evolve under the constraints of negative selection (Figure 7: Panel D).
α-ntxs and κ-ntxs have been theorised to originate from a common ancestral 3FTx subtype gene [6,93]. They are characterised by the highly conserved signal peptide and 3’ non-coding regions and share an apotypic loop-2 stabilising cysteine pair [93,94,95,96]. Moreover, comparison of the structural motifs in type II α-ntxs and the κ-ntxs, particularly the residues around the shared ten cysteines (eight plesiotypic plus the apotypic loop-2 stabilising pair), indicates how closely related the two genes are (Figure 7: Panels B and D). Hence, it is apparent that they have originated from a common ancestral 3FTx subtype gene and the differential rate of evolution has resulted in the observed functional and structural differences between the two homologs. In contrast to α-ntx, κ-bungarotoxin sequences are very well conserved, as a result of the aforementioned structural constraints (Figure 7: Panel D).

3.5. Rapid Accumulation of Variations in Exposed Residues (RAVER) and Focal Mutagenesis in the Evolution of Snake Venom Three-Finger Toxins

Considering the highly conserved organization of the exonic regions of 3FTx, Doley et al. 2008 [16] developed a hypothesis of “Accelerated Segment Switch in Exons to alter Targeting” (ASSET), to explain the molecular mechanism of 3FTx evolution. They postulated that during the evolution of the 3FTx gene, segments in exonic regions have been exchanged with distinctly different ones and the resultant “switching” of segments has generated the coding sequence variation and the diversity of 3FTx. The molecular mechanism of segment “switching” was attributed to a complex scenario involving genetic recombination, splicing variation, and independent recruitment events. We argue that these suggested mechanisms are unlikely to have played a major role in the evolution of three-finger toxins. Recombination cannot account for the tremendous variations observed in extremely short segments of 3FTx; all the observed hypermutations were concentrated in 1–6 adjacent amino acid sites (Figure 6, Figure 7, and Figure 8). Furthermore, recombination would not discriminate between functionally/structurally important and unimportant residues. All 3FTx have highly conserved residues (with eight of the original ten plesiotypic cysteines conserved in all forms), involved in maintaining functional and/or structural stability, and none of the hypermutable sites coincide with these (Figure 6, Figure 7, and Figure 8; Supplementary Table 12). In contrast, evidence supporting the hypothesis that all 3FTx genes originated from the same ancestral gene is strong (Figure 1; Supplementary Figure 1; [15]). Hence, the suggestion of multiple independent recruitment of 3FTx as means to explain the observed variations is implausible. Though alternate splicing or the errors associated with this process can introduce variation in toxin length, like recombination, they cannot account for the considerable variations observed in very short segments of 3FTx. In addition, insertions and/or deletions in shared molecular scaffolds are more plausible explanations for the variations in sequence length.
The role of point mutations in shaping the diversity of 3FTx gene was underestimated by the aforementioned authors [57]. They speculated that point mutations alone could not account for the diversity of 3FTx functional forms and could have only been helpful in fine-tuning receptor binding capabilities, which arise through ASSET [16,57]. A prominent role for point mutations in the evolution of 3FTx gene was discounted based on the following reasoning: (i) point mutations would require numerous generations to account for the observed variations; (ii) all intermediates would have to be selected for via positive selection; and (iii) point mutations would have to occur independently in unrelated lineages to produce the same structurally/functionally important conserved motif.
We argue that these reasons for rejection are not valid. Predatory venom-encoding genes have been demonstrated to evolve under the significant influence of positive Darwinian selection, coupled with focal mutagenesis [6,59,97]. They are likely involved in a coevolutionary arms race, where evolving venom resistance in prey and the evolution of novel venom components exert reciprocal selective pressures on one another [98,99]. The recruitment of the 3FTx gene to reptilian venom systems was previously thought to have occurred within Caenophidia snakes between the divergence of the pareatid snakes approximately 65 million years ago (MYA) and that of the viperid snakes approximately 54 MYA [62]. We have recently recovered 3FTx transcripts from the venom glands of Henophidia snakes, demonstrating that this recruitment event in fact occurred within the ancient Afrophidia snake clade circa 104 MYA [9]. 3FTx genes explosively diversified subsequent to the evolution of the advanced venom delivery apparatus in elapids, 40 MYA (Table 2, Table 3, and Table 4; Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, and Figure 8; Supplementary Figures 2–8). Therefore, an appropriately large number of generations have passed for point mutations, under the influence of positive Darwinian selection to shape the diversity of 3FTxs and account for the tremendous variation observed.
The aforementioned study [16] also failed to take into account the essence of focal mutagenesis, in which point mutations in certain regions of the toxins, such as the molecular surface, may confer an adaptive advantage. In order to conserve catalytic functions, most mutations in proteins concentrate in structurally and/or functionally unimportant regions. Since mutations occur randomly, they are equally likely to occur even within important regions of proteins. However, mutations in structurally or functionally important regions could disrupt the normal functioning of the protein and hence they are filtered out of the population by natural selection. Synthesis and secretion of venom proteins is an energetically expensive process [100,101,102] and hence, natural selection will only permit the rapid accumulation of mutations in certain regions of the toxin, such as the molecular surface, where accumulation of variation could confer an adaptive advantage. Evidently, venom proteins have been demonstrated to evolve via mutation of the surface chemistry, which introduces novel residues (and side-chains) on the molecular surface that could non-specifically interact with the novel targets in prey animals and cause a myriad of pharmacological effects [58,59,103]. In certain organisms like vampire bats, focal mutagenesis in venom proteins can be advantageous as it is likely to delay/prevent the development of immunological resistance in the prey by introducing extremely variable protein surface chemistry in the vampire bat population [61]. Not surprisingly, in vampire bats it has been previously demonstrated that prolonged feeding can result in mounting of immunological resistance in the prey [104], further evidence of the predator–prey (or parasite–host) arms race in which toxins are involved. Mapping of mutations onto three-dimensional homology models and computation of accessible surface area ratios for various three-finger toxins, suggested that a greater proportion of mutations have accumulated on the molecular surface and in the loops of these toxins (Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, and Figure 8; Supplementary Tables 11). We found that the likelihood of hypermutable sites in three-finger toxins occurring on the molecular surface was greater than the likelihood of the site being buried (Table 5).
Hence, we propose a new theory called: “Rapid Accumulation of Variations in Exposed residues” (RAVER) to explain the molecular evolution of the 3FTx gene. We theorise that most predatory venom proteins experience RAVER or focal mutagenesis, which preserves the structural and functional integrity of the toxin, while permitting the rapid accumulation of variations on the molecular surface through positive selection. Thus, while variations accumulate on the molecular surface of the toxin under a coevolutionary arms race scenario, rare mutations in structural and functional residues of rapidly evolving types generate novel molecular scaffolds that get selected only if they considerably improve the potency of the toxin in comparison with the plesiotypic form. Otherwise, such mutations may get eliminated from the population by negative selection. The frequency at which such “rare” mutations accumulate within a functional type likely depends on the strength of positive Darwinian selection. Indeed, rapid accumulation of point mutations under the influence of positive selection in Type II α-ntx gene has led to the origination of −2C Type II α-ntxs on multiple occasions in Oxyuranus and Pseudonaja, two closely related genera of Australian elapid snakes with unusually potent venom (Figure 2).

3.6. Structure–Function Relationship in Three-Finger Toxins

Mapping positively selected sites onto alignments of various three-finger toxins revealed that a large number of residues implicated in structural/functional roles were extremely well conserved (80% ≥ identity; Figure 7 and Figure 8; Supplementary Table 12). However, very few functionally/structurally important residues (A28, F29, R36 and F65) in Type II α-ntxs have been poorly conserved. This is not surprising given the rapid evolution of Type II α-ntxs, which has even resulted in the replacement of two cysteine residues that are usually invariant, and led to the formation of new toxin type (−2C Type II α-ntxs) in Oxyurnaus and Pseudonaja species. It has been demonstrated in erabutoxin (a Type I α-ntx; uniprot: P60775) that mutations at position 29 have a drastic effect on the binding affinity of this toxin to the Torpedo nAChR [105]. Evaluation of mutations at the homologous site in α-cobratoxin (a Type II α-ntx; uniprot: P01391) were shown to have varied effects on binding affinities [106]. Although most sequences had a conserved phenylalanine residue at this position (F29; 49%), a number of others (33%) had tryptophan and histidine residues. It has been demonstrated that these substitutions (F29W and F29H) do not affect the binding affinity of Type II α-ntxs in any way [106]. Two other sequences (1%) had a leucine (F29L) residue in this position, which has been shown to substantially reduce the binding affinity of this toxin [106]. However, lysine residues in loop II (K23) and loop III (K49) have been shown to mediate the binding of Type II α-ntxs to α1-nAChRs and both these residues remain extremely well conserved in these toxins (Figure 7: Panel B; Supplementary Table 12). Thus, α1-nAChR targeting capability of Type II α-ntxs remains intact, despite the mutation at the 29th position. Similarly, it has been demonstrated that A28R mutation does not affect the function of Type II α-ntxs [107] and a large number of sequences (77%) at this site either had alanine or arginine residues. Mutations at R36 and F65 have shown to only slightly affect the binding capabilities of Type II α-ntxs at both neuronal and muscle nAChR. It should be noted that α-ntx possess several extremely well-conserved residues that mediate the specific binding of these toxins to either muscular or neuronal receptors (Figure 7). Thus, despite the mutations at R36 and F65, these toxins can still bind to muscular and neuronal receptors. These examples clearly depict how Type II α-ntxs cope with rapid mutation rate. R33 in Type II α-ntx has been shown to be the most important site for toxicity, as it caused the greatest effect on the function of Type II α-ntx in mutational studies. We noted that this site was extremely well conserved (86% identity; Supplementary Table 12). Mutational studies have also shown that T47 and D53, despite being very closely located to K49 on loop III, are unlikely to play a significant functional role [106]. Interestingly, our selection analyses have detected these sites as hypermutable, among 21 other hypermutable sites in Type II α-ntx (Figure 7: Panel B). Rare mutations at functional sites, such as I26R in erabutoxin (P60775), have been demonstrated to increase the affinity of α-ntx to their receptors [108]. Hence, rapidly evolving α-ntx could derive extremely useful functional variations by conserving most functional residues.
Most hypervariable regions (73%) in Type III α-ntx were concentrated in loop I and loop II (Figure 7), regions that harbour residues implicated in binding of Type I α-ntx to α1-nAChR (Supplementary Table 12). However, unlike Type I α-ntxs, Type III seem to have very well-conserved region between plesiotypic cysteines 9 and 10 [6], which form the 4th disulfide bond in both these toxins. The corresponding region in Type II α-ntx is also very well conserved, and shares with Type III α-ntx three highly conserved amino acids: S54, T55 and D56 (Figure 7). It is very likely that this region, being the only well-conserved region in Type III α-ntx, plays an important structural or functional role. The invariant Asp 56 in this region seems particularly important for the activity of this toxin clade. The observed pattern of accumulation of hypermutations (Figure 7) can also be explained by prey-specific and intraspecific conservation of residues in Type III α-ntx.
κ-bungarotoxins and cytotoxins were unique among all 3FTx lineages examined in this study, in having extremely conserved amino acid sequences (Figure 7 and Figure 8; Supplementary Figure 6). 52% (33 residues) and 65% (53 residues) of sites in κ-bungarotoxins and cytotoxins, respectively, were invariant. Although toxins from these two clades accumulated a few positively selected sites (5 and 2, respectively), they did not overlap with residues implicated in structural/functional stability.
Mapping of hypermutable sites onto the alignments of the plesiotypic 3FTxs reveals a lack of conservation of most residues implicated in receptor binding of Type I and Type II elapid α-neurotoxins, which may be reflective of differences in taxon specificity. However, the residues implicated for activity in Type I and II may be specific to the assay utilised, as plesiotypic 3FTx have shown to be dramatically more potent against birds and/or lizards than mammals, but the specific receptors mediating the activity have not yet been elucidated. The residues implicated in Type I or Type II binding may also be inaccurate as a result of the former prevalence of fish assays (Torpedo species), which may not be reflective of residues actually important in binding to the receptors of natural prey items.
The impact of the development of the venom-delivery system (VDS) and its subsequent specialisation in Caenophidia (advanced snakes) on the evolution of 3FTxs is intriguing. While 3FTx was recruited at the base of the snake evolutionary tree [9], only two NFF lineages (Colubridae and Psammophiinae) and one front-fanged lineage (Elapidae) express this toxin type as the dominant class [18]. As many ‘non-front-fanged’ advanced snakes use venom as part of a combined arsenal that includes variable use of constriction and bite pressure [67] there may therefore be less selective pressure on the diversification of 3FTx within most members of the NFF lineages. Viperids, on the other hand, rely on venom for prey subjugation, and possess arguably the most sophisticated venom-delivery apparatus found amongst snakes (hinged hollow fangs). Viperid snakes utilise a primarily haemolytic and necrotic venom arsenal for prey subjugation and thus 3FTx have not diversified greatly within this family. This may be linked to prey preference, particularly in pit vipers, which possess heat-sensing pits for targeting endothermic (“warm-blooded”) prey. As plesiotypic 3FTx have been shown to be more toxic to birds and reptiles than mammals [21,109], snakes that specialise in feeding on mammals may have less use for this toxin type. Conversely, the earliest elapid snake lineages, with hollow, fixed front fangs, may have fed primarily on reptiles as plesiomorphic elapid snakes do today (e.g., [83,85]), and employed 3FTx as their primary weapon. Hence, three-finger toxins in these snakes underwent rapid and explosive diversification under the significant influence of positive Darwinian selection, resulting in the formation of a myriad of three-finger toxin forms which appear to be exclusive to the elapid snakes (Figure 1). Interestingly, this diversification of 3FTx within elapids led to the origin of new forms that were highly potent against mammals, facilitating more efficient predation on a wider range of prey species including the rodent prey specialised upon by some extant lineages of elapids such as Oxyuranus and Pseudonaja. The plesiotypic 3FTx, which are expressed in relatively lower amounts in elapids, viperids and a majority of “non-front-fanged” advanced snakes (except in Colubridae and Psammophiinae venoms where they are the dominant toxin types), could still play an ancillary envenoming function and thus still be subject to RAVER, potentially facilitating the apotyposis of novel functional forms within this clade. Further sequencing and biochemical activity testing of plesiotypic 3FTx genes will unravel their molecular evolutionary history and the underlying mechanisms of action.

4. Conclusion

In this study, we investigated the following hypotheses: (i) that 3FTx with specific sites of action are involved in coevolutionary arms races with receptors of prey species and thus are evolving under the influence of positive selection; (ii) that cytotoxic 3FTx, which interact non-specifically with cell membranes, are not involved in chemical arms race, and thus are constrained by negative selection; (iii) that, in venoms in which 3FTx are a major component (Elapidae, Colubridae and Psammophiinae), 3FTx are rapidly evolving under positive selection (with the exception of cytotoxic 3FTx in the venom of elapid snakes); (iv) that in venoms in which 3FTx are a minor component (Viperidae), 3FTx evolve under a neutral selection regime and do not experience positive selection; and (v) the evolution of the advanced venom delivery system in elapid snakes led to an increase in positive selection pressure on 3FTx, resulting in the diversity of functional forms present in the venoms of snakes from this family. We were able to confirm hypotheses i, ii, iii and v, but surprisingly the viperid 3FTx were unravelled to have evolved under the influence of positive selection, resulting in the rejection of hypothesis iv.
In addition, we found that point mutations, coupled with focal mutagenesis and positive selection, have played a primary role in the evolution and diversification of snake venom three-finger toxins. We thus reject ASSET as a viable theory for explaining the molecular evolution of three-finger toxins. Instead, we propose “Rapid Accumulation of Variations in Exposed Residues” (RAVER) to explain the dynamic evolution of 3FTx. We elucidate an exquisite molecular mechanism in which the molecular surface, and the regions not implicated in structural/functional stability rapidly diversify, while the core structural/functional residues remain largely constrained by negative selection. However, even these highly conserved residues occasionally accumulate point mutations, especially in toxins with multiple targets such as Type II α-ntx. Rare mutations in these highly conserved residues can generate unique toxin forms with novel folds, which might be more potent than the plesiotypic form.
We have also shown for the first time that not all 3FTx forms are influenced by positive selection and highlighted the differential rate of evolution of various 3FTx types. As well as confirming hypothesis ii, by showing that cytotoxins have been constrained by negative selection, we found that unique forms like the κ-ntx have experienced weaker influence of positive selection than the α-ntx.
These results underscore the dynamic nature of venom evolution, where a myriad of forces, such as morphological refinement of the venom delivery system, shifts in feeding ecology, coevolutionary chemical arms races driving toxin neofunctionalisation via changes in targeting and evolution of toxin complexes with multiple subunits. all influence the molecular evolutionary history of toxins. This generates a molecular complexity that is both a blessing and a curse. On one hand, this complexity complicates snakebite treatment because of highly variable clinical sequelae and relative cross-reactivity of available antivenoms. On the other hand, it provides a rich pool of novel investigational ligands with tremendous potential for biodiscovery.
It should be emphasised that this study only focused upon structural/functional 3FTx variants for which there were sufficient nucleotides to analyse. As shown in Figure 1 and Table 1, there are numerous functional variants that were not studied due to insufficient nucleotide sequences available (with some known only from peptide sequences). In addition, as is clearly evident in Figure 1 and has been covered in-depth by us previously [6], there is a tremendous diversity of 3FTxs for which no functional activity has been determined. Thus, 3FTx represent a tremendous resource not only for investigating the fundamentals of molecular evolution but also as a vast pool of novel peptides for biodiscovery. It is hoped that this study will stimulate further work on this very interesting class of snake toxin.

5. Materials and Methods

5.1. Specimens

Azemiops feae specimen was from Hunan, China.

5.2. Azemiops Feae Library Construction

Total RNA was extracted from the dissected venom gland using the standard TRIzol Plus method (Invitrogen). Extracts were enriched for mRNA using standard RNeasy mRNA mini kit (Qiagen) protocol. mRNA was reverse transcribed, fragmented and ligated to a unique 10-base multiplex identifier (MID) tag prepared using standard protocols and applied to one PicoTitrePlate (PTP) for simultaneous amplification and sequencing on a Roche 454 GS FLX+ Titanium platform (Australian Genome Research Facility). Automated grouping and analysis of sample-specific MID reads informatically separated sequences from the other transcriptomes on the plates, which were then post-processed to remove low quality sequences before de novo assembly into contiguous sequences (contigs) using v 3.4.0.1 of the MIRA software program, resulting in 1656 contigs from 811,967 assembled bases (490 average per contig). Assembled contigs were processed using CLC Main Work Bench (CLC-Bio) and Blast2GO bioinformatic suite [110] to provide Gene Ontology, BLAST and domain/Interpro annotation.

5.3. Sequence Retrieval and Alignment

To assess the molecular evolutionary history of various three-finger toxins, we retrieved nucleotide sequences from National Center for Biotechnology Information [111]. The new toxin sequence recovered from Azemiops feae in this study was identified by comparison of the translated DNA sequence with those of previously characterised toxins using a BLAST search [112] implemented in the UniProtKB protein database. Translated nucleotide sequences were aligned using MUSCLE 3.8 [113] and the alignments were manually inspected to rectify errors. Sequences were refined to remove regions with gaps in more than 50% of the sequences. All alignments are available as supplementary file 1, while all the nucleotide sequences (with their NCBI accession codes) analysed in this study are deposited as supplementary file 2.

5.4. Phylogenetics

Phylogenetic analyses were performed to allow the reconstruction of the molecular evolutionary history of various types of three-finger toxins. Since many three-finger toxin forms are only known from their amino acid sequences, amino acid datasets were used for the Bayesian inference. Trees were generated using both Bayesian and maximum-likelihood inferences. MrBayes 3.2.1 [114] was employed for the Bayesian inference, using lset rates = invgamma with prset aamodelpr = mixed command, which enables the program to optimise between nine different amino acid substitution matrices implemented in MrBayes. WAG [115] was chosen as the best substitution matrix by the program. The analysis was performed by running a minimum of 1×107 generations in four chains, and saving every 100th tree. The log-likelihood score of each saved tree was plotted against the number of generations to establish the point at which the log likelihood scores reached their asymptote, and the posterior probabilities for clades established by constructing a majority-rule consensus tree for all trees generated after completion of the burn-in phase. Maximum-likelihood method implemented in PhyML [116] was also employed on the nucleotide datasets and node support was evaluated with 1,000 bootstrapping replicates.

5.5. Test for Recombination

To overcome the effects of recombination on the phylogenetic and evolutionary interpretations [117], we employed Single Breakpoint algorithms implemented in the HyPhy package and assessed recombination on all the toxin forms examined in this study [118,119]. When potential breakpoints were detected using the small sample Akaike information Criterion (AICc), the multiple sequence alignments were compartmentalised before conducting selection analyses, so as to allow the different recombinant regions to have distinct phylogenetic topology.

5.6. Selection Analyses

We evaluated the influence of natural selection on various types of three-finger toxins using maximum-likelihood models [120,121] implemented in Codeml of the PAML package [122]. We utilised the three-ratio model as well as the optimised branch-site test [123,124] to detect evolutionary selection pressures acting upon various toxin lineages. The three-ratio model evaluates selection across the lineages alone, while the branch-site model allows ω to vary across the sites of the protein and along the branches in the tree. The latter is known for its reasonable power and accuracy to detect short bursts of episodic adaptations [123]. However, both the three-ratio and branch-site models require the foreground (lineages under positive selection) and background lineages (lineages lacking positive selection) to be defined a priori. Such predefined biological hypotheses are often unavailable and it becomes difficult to define the foreground lineages. Therefore, we treated each 3FTx type being compared as a foreground branch alternatively and tested multiple hypotheses. A likelihood-ratio test was then conducted by comparing the model that allows ω to be greater than 1 in the foreground branch, with a null model that does not (ω constrained 1). It has been suggested that while implementing multiple hypotheses using branch and branch-site models, it is necessary to control the family-wise error rate (FWER or Type I error) [123]. We used Bonferroni correction to account for such errors. It uses α/n as the significance level to test each hypothesis, where “α” is the significance level and “n” is the number of independent true null hypotheses. We further utilised branch-site Random Effects Likelihood (REL) [63] to identify lineages evolving under the influence of episodic diversifying selection pressures. Unlike the aforementioned lineage-specific models, branch-site REL does not require a priori identification of foreground and background branches. Moreover, branch-site tests have been shown to be prone to false positives when the background lineages are also under the influence of positive selection and branch-site test was designed to address these shortcomings [63].
Since lineage-specific models assess the effects of selection only along the branches, they often fail to identify regions in proteins that might be affected by episodic selection pressures and ultimately underestimate the strength of selection. Hence, we employed site-specific models that estimate positive selection statistically as a non-synonymous-to-synonymous nucleotide-substitution rate ratio (ω) significantly greater than 1. We compared likelihood values for three pairs of models with different assumed ω distributions as no a priori expectation exists for the same: M0 (constant ω rates across all sites) versus M3 (allows ω to vary across sites within “n” discrete categories, n ≥ 3); M1a (a model of neutral evolution) where all sites are assumed to be either under negative (ω < 1) or neutral selection (ω = 1) versus M2a (a model of positive selection) which in addition to the site classes mentioned for M1a, assumes a third category of sites; sites with ω > 1 (positive selection) and M7 (Beta) versus M8 (Beta and ω), and models that mirror the evolutionary constraints of M1 and M2 but assume that ω values are drawn from a beta distribution [125]. Only if the alternative models (M3, M2a and M8: allow sites with ω > 1) show a better fit in Likelihood Ratio Test (LRT) relative to their null models (M0, M1a and M7: do not allow sites ω > 1), are their results considered significant. It should be noted that although we have reported the results of all site-specific models (M2a, M3 and M8), only the results of M8, the most accurate model among these, were considered for downstream analyses (Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, and Figure 8; Table 4 and Table 5). LRT is estimated as twice the difference in maximum likelihood values between nested models and compared with the χ2 distribution with the appropriate degree of freedom—the difference in the number of parameters between the two models. The Bayes empirical Bayes (BEB) approach [126] was used to identify amino acids under positive selection by calculating the posterior probabilities that a particular amino acid belongs to a given selection class (neutral, conserved or highly variable). Sites with greater posterior probability (PP ≥ 95%) of belonging to the “ω > 1 class” were inferred to be positively selected.
FUBAR [127] implemented in HyPhy [128] was employed to detect codon sites evolving under the influence of pervasive diversifying and purifying selection pressures. The Mixed Effects Model of Evolution (MEME) [129] was also employed to efficiently detect sites that experience diversifying selection for a short period in the evolutionary timescale. Non-synonymous mutations that introduce extremely variant amino acids are more likely to influence the structure-function, and hence the fitness of the organism, relative to mutations that introduce amino acids with similar side-chains as the ancestral residues. Further support for the results of the nucleotide-level selection analyses was obtained and the radicalness of mutations were assessed using a complementary protein-level approach implemented in TreeSAAP [130].
Direct comparison of ω values computed from different datasets can be misleading, as they can have different proportions of sites under selection. Hence, we assessed the selection pressures shaping different 3FTx clades by employing clade model analyses implemented in Codeml and simultaneously estimated ω values [131]. The significance of the analysis was tested by comparing the likelihood of this model with that of model M1a. To clearly depict the proportion of sites under different regimes of selection, an evolutionary fingerprint analysis was carried out using the evolutionary selection distance (ESD) algorithm implemented in datamonkey [132].

5.7. Structural Analyses

To depict the natural selection pressures influencing the evolution of various three-finger toxins, we mapped the sites under positive selection on the homology models created using the Phyre 2 webserver [133]. Pymol 1.3 [134] was used to visualise and generate the images of homology models. The Consurf webserver [135] was used for mapping the evolutionary selection pressures on the three-dimensional homology models. GETAREA [136] was used to calculate the Accessible Surface Area (ASA) or the solvent exposure of amino-acid side chains. It uses the atom co-ordinates of the PDB file and indicates if a residue is buried or exposed to the surrounding medium by comparing the ratio between side chain ASA and the “random coil” values per residue. An amino acid is considered to be buried if it has an ASA less than 20% and exposed if the ASA is more than or equal to 50%. When ASA ratio lies between 40% and 50%, it is highly likely that the residues have their side chains exposed to the surrounding medium.

Acknowledgements

KS was funded by the PhD grant (SFRH/BD/61959/2009) from F.C.T (Fundação para a Ciência e a Tecnologia). TNWJ was funded by a PhD scholarship from the University of Queensland. AA was funded by the European Regional Development Fund (ERDF) through the COMPETE - Operational Competitiveness Programme and national funds through F.C.T under the projects PEst-C/MAR/LA0015/2013 and PTDC/AAC-AMB/121301/2010 (FCOMP-01-0124-FEDER-019490), and BGF was funded by the Australian Research Council (ARC) and the University of Queensland. SAA was the recipient of a postdoctoral fellowship (PDRF Phase II Batch-V) from Higher Education Commission (HEC Islamabad) Pakistan. EABU acknowledges funding from the University of Queensland (International Postgraduate Research Scholarship, UQ Centennial Scholarship, and UQ Advantage Top-Up Scholarship) and the Norwegian State Education Loans Fund.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fry, B.G.; Roelants, K.; Champagne, D.E.; Scheib, H.; Tyndall, J.D.A.; King, G.F.; Nevalainen, T.J.; Norman, J.A.; Lewis, R.J.; Norton, R.S.; et al. The toxicogenomic multiverse: Convergent recruitment of proteins into animal venoms. Annu. Rev. Genomics Hum. Genet. 2009, 10, 483–511. [Google Scholar] [CrossRef]
  2. Escoubas, P.; Sollod, B.; King, G.F. Venom landscapes: Mining the complexity of spider venoms via a combined cDNA and mass spectrometric approach. Toxicon 2006, 47, 650–663. [Google Scholar] [CrossRef]
  3. Fry, B.G. From genome to “venome”: Molecular origin and evolution of the snake venom proteome inferred from phylogenetic analysis of toxin sequences and related body proteins. Genome Res. 2005, 15, 403–420. [Google Scholar] [CrossRef]
  4. Inceoglu, B.; Lango, J.; Jing, J.; Chen, L.; Doymaz, F.; Pessah, I.N.; Hammock, B.D. One scorpion, two venoms: Prevenom of Parabuthus transvaalicus acts as an alternative type of venom with distinct mechanism of action. Proc. Natl. Acad. Sci. USA 2003, 100, 922–927. [Google Scholar]
  5. Olivera, B.M. CONUS VENOM PEPTIDES: Reflections from the biology of clades and species. Annu. Rev. Ecol. Syst. 2002, 33, 25–47. [Google Scholar] [CrossRef]
  6. Fry, B.G.; Wüster, W.; Kini, R.M.; Brusic, V.; Khan, A.; Venkataraman, D.; Rooney, A.P. Molecular evolution and phylogeny of elapid snake venom three-finger toxins. J. Mol. Evol. 2003, 57, 110–129. [Google Scholar]
  7. Fry, B.G.; Casewell, N.R.; Wüster, W.; Vidal, N.; Young, B.; Jackson, T.N.W. The structural and functional diversification of the Toxicofera reptile venom system. Toxicon 2012, 60, 434–448. [Google Scholar]
  8. Chang, D.; Duda, T.F. Extensive and continuous duplication facilitates rapid evolution and diversification of gene families. Mol. Biol. Evol. 2012, 29, 2019–2029. [Google Scholar]
  9. Fry, B.G.; Undheim, E.A.B.; Ali, S.A.; Jackson, T.N.W.; Debono, J.; Scheib, H.; Ruder, T.; Morgenstern, D.; Cadwallader, L.; Whitehead, D.; et al. Squeezers and leaf-cutters: differential diversification and degeneration of the venom system in toxicoferan reptiles. Mol. Cell. Proteomics 2013, 12, 1881–1899. [Google Scholar] [CrossRef]
  10. Kozminsky-Atias, A.; Zilberberg, N. Molding the business end of neurotoxins by diversifying evolution. FASEB J. 2012, 26, 576–586. [Google Scholar] [CrossRef]
  11. Casewell, N.R.; Wagstaff, S.C.; Harrison, R.A.; Renjifo, C.; Wuster, W. Domain loss facilitates accelerated evolution and neofunctionalization of duplicate snake venom metalloproteinase toxin genes. Mol. Biol. Evol. 2011, 28, 2637–2649. [Google Scholar] [CrossRef]
  12. Weinberger, H.; Moran, Y.; Gordon, D.; Turkov, M.; Kahn, R.; Gurevitz, M. Positions under positive selection--key for selectivity and potency of scorpion alpha-toxins. Mol. Biol. Evol. 2010, 27, 1025–1034. [Google Scholar]
  13. Zhu, S.; Peigneur, S.; Gao, B.; Lu, X.; Cao, C.; Tytgat, J. Evolutionary diversification of Mesobuthus α-scorpion toxins affecting sodium channels. Mol. Cell. Proteomics 2012, 11, M111.012054. [Google Scholar]
  14. Nei, M.; Gu, X.; Sitnikova, T. Evolution by the birth-and-death process in multigene families of the vertebrate immune system. Proc. Natl. Acad. Sci. USA 1997, 94, 7799–7806. [Google Scholar]
  15. Fry, B.G.; Lumsden, N.G.; Wüster, W.; Wickramaratna, J.C.; Hodgson, W.C.; Kini, R.M. Isolation of a neurotoxin (alpha-colubritoxin) from a nonvenomous colubrid: evidence for early origin of venom in snakes. J. Mol. Evol. 2003, 57, 446–452. [Google Scholar] [CrossRef]
  16. Doley, R.; Pahari, S.; Mackessy, S.P.; Kini, R.M. Accelerated exchange of exon segments in Viperid three-finger toxin genes (Sistrurus catenatus edwardsii; Desert Massasauga). BMC Evol. Biol. 2008, 8, 196. [Google Scholar]
  17. Fry, B.G.; Scheib, H.; van der Weerd, L.; Young, B.; McNaughtan, J.; Ramjan, S.F.R.; Vidal, N.; Poelmann, R.E.; Norman, J.A. Evolution of an arsenal: Structural and functional diversification of the venom system in the advanced snakes (Caenophidia). Mol. Cell. Proteomics 2008, 7, 215–246. [Google Scholar]
  18. Fry, B.G.; Wüster, W.; Ryan Ramjan, S.F.; Jackson, T.; Martelli, P.; Kini, R.M. Analysis of Colubroidea snake venoms by liquid chromatography with mass spectrometry: Evolutionary and toxinological implications. Rapid Commun. Mass Spectrom. 2003, 17, 2047–2062. [Google Scholar]
  19. Lumsden, N.G.; Fry, B.G.; Ventura, S.; Kini, R.M.; Hodgson, W.C. The in vitro and in vivo pharmacological activity of Boiga dendrophila (mangrove catsnake) venom. Auton. Autacoid Pharmacol. 2004, 24, 107–113. [Google Scholar]
  20. Pawlak, J.; Mackessy, S.P.; Fry, B.G.; Bhatia, M.; Mourier, G.; Fruchart-Gaillard, C.; Servent, D.; Ménez, R.; Stura, E.; Ménez, A.; et al. Denmotoxin, a three-finger toxin from the colubrid snake Boiga dendrophila (Mangrove Catsnake) with bird-specific activity. J. Biol. Chem. 2006, 281, 29030–29041. [Google Scholar] [CrossRef]
  21. Pawlak, J.; Mackessy, S.P.; Sixberry, N.M.; Stura, E.A.; le Du, M.H.; Ménez, R.; Foo, C.S.; Ménez, A.; Nirthanan, S.; Kini, R.M. Irditoxin, a novel covalently linked heterodimeric three-finger toxin with high taxon-specific neurotoxicity. FASEB J. 2009, 23, 534–545. [Google Scholar]
  22. Lumsden, N.G.; Fry, B.G.; Manjunatha Kini, R.; Hodgson, W.C. In vitro neuromuscular activity of “colubrid” venoms: Clinical and evolutionary implications. Toxicon 2004, 43, 819–827. [Google Scholar]
  23. Kini, R.M.; Doley, R. Structure, function and evolution of three-finger toxins: Mini proteins with multiple targets. Toxicon 2010, 56, 855–867. [Google Scholar] [CrossRef]
  24. Pahari, S.; Mackessy, S.P.; Kini, R.M. The venom gland transcriptome of the Desert Massasauga Rattlesnake (Sistrurus catenatus edwardsii): Towards an understanding of venom composition among advanced snakes (Superfamily Colubroidea). BMC Mol. Biol. 2007, 8, 115. [Google Scholar] [CrossRef]
  25. Utkin, Y.N.; Kukhtina, V.V.; Kryukova, E.V.; Chiodini, F.; Bertrand, D.; Methfessel, C.; Tsetlin, V.I. “Weak Toxin” from naja kaouthia is a nontoxic antagonist of α7 and muscle-type nicotinic acetylcholine receptors. J. Biol. Chem. 2001, 276, 15810–15815. [Google Scholar] [CrossRef]
  26. Rochelle, M.J.; Kardong, K.V. Constriction versus envenomation in prey capture by the brown tree snake, boiga irregularis (Squamata: Colubridae). Herpetologica 1993, 49, 301–304. [Google Scholar]
  27. Barber, C.M.; Isbister, G.K.; Hodgson, W.C. Alpha neurotoxins. Toxicon 2013, 66, 47–58. [Google Scholar]
  28. Gong, N.; Armugam, A.; Jeyaseelan, K. Postsynaptic short-chain neurotoxins from Pseudonaja textilis. cDNA cloning, expression and protein characterization. Eur. J. Biochem. 1999, 265, 982–989. [Google Scholar]
  29. Chiappinelli, V.A.; Weaver, W.R.; McLane, K.E.; Conti-Fine, B.M.; Fiordalisi, J.J.; Grant, G.A. Binding of native kappa-neurotoxins and site-directed mutants to nicotinic acetylcholine receptors. Toxicon 1996, 34, 1243–1256. [Google Scholar] [CrossRef]
  30. Carsi, J.M.; Potter, L.T. m1-toxin isotoxins from the green mamba (Dendroaspis angusticeps) that selectively block m1 muscarinic receptors. Toxicon 2000, 38, 187–198. [Google Scholar] [CrossRef]
  31. Fruchart-Gaillard, C.; Mourier, G.; Marquer, C.; Stura, E.; Birdsall, N.J.M.; Servent, D. Different interactions between MT7 toxin and the human muscarinic M1 receptor in its free and N-methylscopolamine-occupied states. Mol. Pharmacol. 2008, 74, 1554–1563. [Google Scholar] [CrossRef]
  32. Jianmongkol, S.; Marable, B.R.; Berkman, C.E.; Talley, T.T.; Thompson, C.M.; Richardson, R.J. Kinetic evidence for different mechanisms of acetylcholinesterase inhibition by (1R)- and (1S)-stereoisomers of isomalathion. Toxicol. Appl. Pharmacol. 1999, 155, 43–53. [Google Scholar] [CrossRef]
  33. Jolkkonen, M.; van Giersbergen, P.L.; Hellman, U.; Wernstedt, C.; Karlsson, E. A toxin from the green mamba Dendroaspis angusticeps: Amino acid sequence and selectivity for muscarinic m4 receptors. FEBS Lett. 1994, 352, 91–94. [Google Scholar] [CrossRef]
  34. Karlsson, E.; Jolkkonen, M.; Satyapan, N.; Adem, A.; Kumlin, E.; Hellman, U.; Wernstedt, C. Protein toxins that bind to muscarinic acetylcholine receptors. Ann. N. Y. Acad. Sci. 1994, 710, 153–161. [Google Scholar] [CrossRef]
  35. Koivula, K.; Rondinelli, S.; Nasman, J. The three-finger toxin MTalpha is a selective alpha(2B)-adrenoceptor antagonist. Toxicon 2010, 56, 440–447. [Google Scholar] [CrossRef]
  36. Kornisiuk, E.; Jerusalinsky, D.; Cerveñansky, C.; Harvey, A.L. Binding of muscarinic toxins MTx1 and MTx2 from the venom of the green mamba Dendroaspis angusticeps to cloned human muscarinic cholinoceptors. Toxicon 1995, 33, 11–18. [Google Scholar] [CrossRef]
  37. Krajewski, J.L.; Dickerson, I.M.; Potter, L.T. Site-directed mutagenesis of m1-toxin1: Two amino acids responsible for stable toxin binding to M(1) muscarinic receptors. Mol. Pharmacol. 2001, 60, 725–731. [Google Scholar]
  38. Kukhtina, V.V.; Weise, C.; Muranova, T.A.; Starkov, V.G.; Franke, P.; Hucho, F.; Wnendt, S.; Gillen, C.; Tsetlin, V.I.; Utkin, Y.N. Muscarinic toxin-like proteins from cobra venom. Eur. J. Biochem. 2000, 267, 6784–6789. [Google Scholar] [CrossRef]
  39. Quinton, L.; Girard, E.; Maiga, A.; Rekik, M.; Lluel, P.; Masuyer, G.; Larregola, M.; Marquer, C.; Ciolek, J.; Magnin, T.; et al. Isolation and pharmacological characterization of AdTx1, a natural peptide displaying specific insurmountable antagonism of the alpha1A-adrenoceptor. Br. J. Pharmacol. 2010, 159, 316–325. [Google Scholar] [CrossRef]
  40. Rouget, C.; Quinton, L.; Maïga, A.; Gales, C.; Masuyer, G.; Malosse, C.; Chamot-Rooke, J.; Thai, R.; Mourier, G.; de Pauw, E.; et al. Identification of a novel snake peptide toxin displaying high affinity and antagonist behaviour for the α2-adrenoceptors. Br. J. Pharmacol. 2010, 161, 1361–1374. [Google Scholar] [CrossRef]
  41. Toomela, T.; Jolkkonen, M.; Rinken, A.; Järv, J.; Karlsson, E. Two-step binding of green mamba toxin to muscarinic acetylcholine receptor. FEBS Lett. 1994, 352, 95–97. [Google Scholar] [CrossRef]
  42. Carsi, J.M.; Valentine, H.H.; Potter, L.T. m2-toxin: A selective ligand for M2 muscarinic receptors. Mol. Pharmacol. 1999, 56, 933–937. [Google Scholar]
  43. Diochot, S.; Baron, A.; Salinas, M.; Douguet, D.; Scarzello, S.; Dabert-Gay, A.S.; Debayle, D.; Friend, V.; Alloui, A.; Lazdunski, M.; et al. Black mamba venom peptides target acid-sensing ion channels to abolish pain. Nature 2012, 490, 552–555. [Google Scholar] [CrossRef]
  44. De Weille, J.R.; Schweitz, H.; Maes, P.; Tartar, A.; Lazdunski, M. Calciseptine, a peptide isolated from black mamba venom, is a specific blocker of the L-type calcium channel. Proc. Natl. Acad. Sci. USA 1991, 88, 2437–2440. [Google Scholar] [CrossRef]
  45. Karlsson, E.; Mbugua, P.M.; Rodriguez-Ithurralde, D. Fasciculins, anticholinesterase toxins from the venom of the green mamba Dendroaspis angusticeps. J. Physiol. (Paris) 1984, 79, 232–240. [Google Scholar]
  46. McDowell, R.S.; Dennis, M.S.; Louie, A.; Shuster, M.; Mulkerrin, M.G.; Lazarus, R.A. Mambin, a potent glycoprotein IIb-IIIa antagonist and platelet aggregation inhibitor structurally related to the short neurotoxins. Biochemistry 1992, 31, 4766–4772. [Google Scholar]
  47. Konshina, A.G.; Dubovskii, P.V.; Efremov, R.G. Structure and dynamics of cardiotoxins. Curr. Protein Pept. Sci. 2012, 13, 570–584. [Google Scholar] [CrossRef]
  48. Joubert, F.J.; Viljoen, C.C. Snake venom. The amino-acid sequence of the subunits of two reduced and S-carboxymethylated proteins (C8S2 and C9S3) from Dendroaspis angusticeps venom. Hoppe. Seylers. Z. Physiol. Chem. 1979, 360, 1075–1090. [Google Scholar] [CrossRef]
  49. Uniprot. Available online: www.uniprot.org (accessed on 14 November 2013).
  50. St Pierre, L.; Fischer, H.; Adams, D.J.; Schenning, M.; Lavidis, N.; de Jersey, J.; Masci, P.P.; Lavin, M.F. Distinct activities of novel neurotoxins from Australian venomous snakes for nicotinic acetylcholine receptors. Cell. Mol. Life Sci. 2007, 64, 2829–2840. [Google Scholar] [CrossRef]
  51. Bhaskaran, R.; Huang, C.C.; Chang, D.K.; Yu, C. Cardiotoxin III from the Taiwan cobra (Naja naja atra). Determination of structure in solution and comparison with short neurotoxins. J. Mol. Biol. 1994, 235, 1291–1301. [Google Scholar] [CrossRef]
  52. Harvey, A.L.; Marshall, R.J.; Karlsson, E. Effects of purified cardiotoxins from the Thailand cobra (Naja naja siamensis) on isolated skeletal and cardiac muscle preparations. Toxicon 1982, 20, 379–396. [Google Scholar] [CrossRef]
  53. Bougis, P.E.; Khélif, A.; Rochat, H. On the inhibition of [Na+,K+]-ATPases by the components of Naja mossambica mossambica venom: Evidence for two distinct rat brain [Na+,K+]-ATPase activities. Biochemistry 1989, 28, 3037–3043. [Google Scholar] [CrossRef]
  54. Gong, N.; Armugam, A.; Jeyaseelan, K. Molecular cloning, characterization and evolution of the gene encoding a new group of short-chain alpha-neurotoxins in an Australian elapid, Pseudonaja textilis. FEBS Lett. 2000, 473, 303–310. [Google Scholar] [CrossRef]
  55. Jiang, Y.; Li, Y.; Lee, W.; Xu, X.; Zhang, Y.Y.; Zhao, R.; Wang, W. Venom gland transcriptomes of two elapid snakes (Bungarus multicinctus and Naja atra) and evolution of toxin genes. BMC Genomics 2011, 12, 1. [Google Scholar] [CrossRef]
  56. Terrat, Y.; Sunagar, K.; Fry, B.; Jackson, T.; Scheib, H.; Fourmy, R.; Verdenaud, M.; Blanchet, G.; Antunes, A.; Ducancel, F. Atractaspis aterrima Toxins: The First Insight into the Molecular Evolution of Venom in Side-Stabbers. Toxins 2013, 5, 1948–1964. [Google Scholar] [CrossRef]
  57. Doley, R.; Mackessy, S.P.; Kini, R.M. Role of accelerated segment switch in exons to alter targeting (ASSET) in the molecular evolution of snake venom proteins. BMC Evol. Biol. 2009, 9, 146. [Google Scholar] [CrossRef]
  58. Sunagar, K.; Johnson, W.E.; O’Brien, S.J.; Vasconcelos, V.; Antunes, A. Evolution of CRISPs associated with toxicoferan-reptilian venom and mammalian reproduction. Mol. Biol. Evol. 2012, 29, 1807–1822. [Google Scholar] [CrossRef]
  59. Brust, A.; Sunagar, K.; Undheim, E.A.B.; Vetter, I.; Yang, D.C.; Casewell, N.R.; Jackson, T.N.W.; Koludarov, I.; Alewood, P.F.; Hodgson, W.C.; et al. Differential evolution and neofunctionalization of snake venom metalloprotease domains. Mol. Cell. Proteomics 2013, 12, 651–663. [Google Scholar] [CrossRef]
  60. Sunagar, K.; Fry, B.G.; Jackson, T.; Casewell, N.R.; Undheim, E.A.B.; Vidal, N.; Ali, S.; King, G.F.; Vasudevan, K.; Vasconcelos, V.; et al. Molecular evolution of vertebrate neurotrophins: Co-option of the highly conserved nerve growth factor gene into the advanced snake venom arsenal. PLoS One 2013, in press. [Google Scholar]
  61. Low, D.H.W.; Sunagar, K.; Undheim, E.A.B.; Ali, S.A.; Alagon, A.C.; Ruder, T.; Jackson, T.N.W.; Pineda Gonzalez, S.; King, G.F.; Jones, A.; Antunes, A.; Fry, B.G. Dracula’s children: Molecular evolution of vampire bat venom. J. Proteomics 2013, 89, 95–111. [Google Scholar] [CrossRef]
  62. Vidal, N.; Rage, J.; Couloux, A.; Hedges, S.B. Snakes (Serpentes). In The Timetree of Life; Hedges, S.B., Kumar, S., Eds.; Oxford University Press: New York, NY, USA, 2009; pp. 390–397. [Google Scholar]
  63. Kosakovsky Pond, S.L.; Murrell, B.; Fourment, M.; Frost, S.D.W.; Delport, W.; Scheffler, K. A random effects branch-site model for detecting episodic diversifying selection. Mol. Biol. Evol. 2011, 28, 3033–3043. [Google Scholar] [CrossRef]
  64. Fry, B.G.; Vidal, N.; Norman, J.A.; Vonk, F.J.; Scheib, H.; Ramjan, S.F.R.; Kuruppu, S.; Fung, K.; Hedges, S.B.; Richardson, M.K.; et al. Early evolution of the venom system in lizards and snakes. Nature 2006, 439, 584–588. [Google Scholar] [CrossRef]
  65. Fry, B.G.; Vidal, N.; van der Weerd, L.; Kochva, E.; Renjifo, C. Evolution and diversification of the Toxicofera reptile venom system. J. Proteomics 2009, 72, 127–136. [Google Scholar] [CrossRef]
  66. Fry, B.G.; Winter, K.; Norman, J.A.; Roelants, K.; Nabuurs, R.J.A.; van Osch, M.J.P.; Teeuwisse, W.M.; van der Weerd, L.; McNaughtan, J.E.; Kwok, H.F.; et al. Functional and structural diversification of the Anguimorpha lizard venom system. Mol. Cell. Proteomics 2010, 9, 2369–2390. [Google Scholar] [CrossRef]
  67. Fry, B.G.; Scheib, H.; de L M Junqueira de Azevedo, I.; Silva, D.A.; Casewell, N.R. Novel transcripts in the maxillary venom glands of advanced snakes. Toxicon 2012, 59, 696–708. [Google Scholar] [CrossRef]
  68. Li, M.; Fry, B.G.; Kini, R.M. Putting the brakes on snake venom evolution: The unique molecular evolutionary patterns of Aipysurus eydouxii (Marbled sea snake) phospholipase A2 toxins. Mol. Biol. Evol. 2005, 22, 934–941. [Google Scholar] [CrossRef]
  69. Junqueira de Azevedo, I.L.; Farsky, S.H.; Oliveira, M.L.; Ho, P.L. Molecular cloning and expression of a functional snake venom vascular endothelium growth factor (VEGF) from the Bothrops insularis pit viper. A new member of the VEGF family of proteins. J. Biol. Chem. 2001, 276, 39836–39842. [Google Scholar]
  70. Deufel, A.; Cundall, D. Prey transport in “palatine-erecting” elapid snakes. J. Morphol. 2003, 258, 358–375. [Google Scholar] [CrossRef]
  71. Deufel, A.; Cundall, D. Feeding in Atractaspis (Serpentes: Atractaspididae): A study in conflicting functional constraints. Zoology (Jena) 2003, 106, 43–61. [Google Scholar] [CrossRef]
  72. Gopalakrishnakone, P.; Kochva, E. Venom glands and some associated muscles in sea snakes. J. Morphol. 1990, 205, 85–96. [Google Scholar] [CrossRef]
  73. McCarthy, C.J. Adaptations of sea snakes that eat fish eggs; with a note on the throat musculature of Aipysurus eydouxi (Gray, 1849). J. Nat. Hist. 1987, 21, 1119–1128. [Google Scholar] [CrossRef]
  74. Salomão, M.G.; Laporta-Ferreira, I.L. The role of secretions from the supralabial, infralabial, and Duvernoy’s glands of the slug-eating snake Sibynomorphus mikani (Colubridae: Dipsadinae) in the immobilization of molluscan prey. J. Herpetol. 1994, 28, 369–371. [Google Scholar] [CrossRef]
  75. Voris, H.K. Fish eggs as the apparent sole food item for a genus of sea snake, emydocephalus (krefft). Ecology 1966, 47, 152. [Google Scholar] [CrossRef]
  76. VORIS, H.K.; VORIS, H.H. Feeding strategies in marine snakes: An analysis of evolutionary, morphological, behavioral and ecological relationships. Integr. Comp. Biol. 1983, 23, 411–425. [Google Scholar] [CrossRef]
  77. Limpus, C.J. The venom apparatus and venom yields of sub-tropical Queensland Hydrophiidae. Toxicon Suppl. 1978, 39–70. [Google Scholar]
  78. Cundall, D.; Deufel, A. Influence of the venom delivery system on intraoral prey transport in snakes. Zool. Anzeiger - A J. Comp. Zool. 2006, 245, 193–210. [Google Scholar] [CrossRef]
  79. Deufel, A.; Cundall, D. Functional plasticity of the venom delivery system in snakes with a focus on the poststrike prey release behavior. Zool. Anz. 2006, 245, 249–267. [Google Scholar] [CrossRef]
  80. Jackson, T.N.W.; Casewell, N.R.; Fry, B.G. Response to “Replies to Fry et al. (Toxicon 2012, 60/4, 434-448). Part A. Analyses of squamate reptile oral glands and their products: A call for caution in formal assignment of terminology designating biological function”. Toxicon 2013, 64, 106–112. [Google Scholar] [CrossRef]
  81. Jackson, T.N.W.; Casewell, N.R.; Fry, B.G. Response to “Replies to Fry et al. (Toxicon 2012, 60/4, 434–448). Part B. Analyses of squamate reptile oral glands and their products: A call for caution in formal assignment of terminology designating biological function”. Toxicon 2013, 64, 113–115. [Google Scholar] [CrossRef]
  82. Prado-Franceschi, J.; Hyslop, S.; Cogo, J.C.; Andrade, A.L.; Assakura, M.; Cruz-Hofling, M.A.; Rodrigues-Simioni, L. The effects of Duvernoy’s gland secretion from the xenodontine colubrid Philodryas olfersii on striated muscle and the neuromuscular junction: Partial characterization of a neuromuscular fraction. Toxicon 1996, 34, 459–466. [Google Scholar] [CrossRef]
  83. Shine, R. Comparative ecology of three australian snake species of the genus cacophis (Serpentes: Elapidae). Copeia 1980, 1980, 831–838. [Google Scholar] [CrossRef]
  84. Shine, R. Reproduction, feeding and growth in the australian burrowing snake vermicella annulata. J. Herpetol. 1980, 14, 71–77. [Google Scholar] [CrossRef]
  85. Shine, R. Ecology of Australian elapid snakes of the genera furina and glyphodon. J. Herpetol. 1981, 15, 219. [Google Scholar] [CrossRef]
  86. Shine, R. Food habits and reproductive biology of small australian snakes of the genera unechis and suta (Elapidae). J. Herpetol. 1988, 22, 307. [Google Scholar] [CrossRef]
  87. Shine, R. Constraints, allometry, and adaptation - food-habits and reproductive-biology of Australian brownsnakes (Pseudonaja, Elapidae). Herpetologica 1989, 45, 195–207. [Google Scholar]
  88. St Pierre, L.; Birrell, G.W.; Earl, S.T.; Wallis, T.P.; Gorman, J.J.; de Jersey, J.; Masci, P.P.; Lavin, M.F. Diversity of toxic components from the venom of the evolutionarily distinct black whip snake, Demansia vestigiata. J. Proteome Res. 2007, 6, 3093–3107. [Google Scholar] [CrossRef]
  89. Condrea, E. Hemolytic disorders associated with a primary red cell membrane defect. Experientia 1976, 32, 537–542. [Google Scholar] [CrossRef]
  90. Dufton, M.J.; Hider, R.C. Structure and pharmacology of elapid cytotoxins. Pharmacol. Ther. 1988, 36, 1–40. [Google Scholar] [CrossRef]
  91. Grant, G.A.; Chiappinelli, V.A. kappa-Bungarotoxin: Complete amino acid sequence of a neuronal nicotinic receptor probe. Biochemistry 1985, 24, 1532–1537. [Google Scholar] [CrossRef]
  92. Dewan, J.C.; Grant, G.A.; Sacchettini, J.C. Crystal structure of kappa-bungarotoxin at 2.3-A resolution. Biochemistry 1994, 33, 13147–13154. [Google Scholar] [CrossRef]
  93. Chang, L.-S.S.; Chung, C.; Lin, J.; Hong, E. Organization and phylogenetic analysis of kappa-bungarotoxin genes from Bungarus multicinctus (Taiwan banded krait). Genetica 2002, 115, 213–221. [Google Scholar] [CrossRef]
  94. Chang, L.; Lin, S.; Huang, H.; Hsiao, M. Genetic organization of alpha-bungarotoxins from Bungarus multicinctus (Taiwan banded krait): Evidence showing that the production of alpha-bungarotoxin isotoxins is not derived from edited mRNAs. Nucleic Acids Res. 1999, 27, 3970–3975. [Google Scholar] [CrossRef]
  95. Chang, L.; Lin, J.; Wu, P.; Chang, C.; Hong, E. cDNA sequence analysis and expression of kappa-bungarotoxin from Taiwan banded krait. Biochem. Biophys. Res. Commun. 1997, 230, 192–195. [Google Scholar] [CrossRef]
  96. Danse, J.M.; Garnier, J.M. cDNA deduced amino-acid sequences of two novel kappa-neurotoxins from Bungarus multicinctus. Nucleic Acids Res. 1990, 18, 1050. [Google Scholar] [CrossRef]
  97. Koludarov, I.; Sunagar, K.; Undheim, E.A.B.; Jackson, T.N.W.; Ruder, T.; Whitehead, D.; Saucedo, A.C.; Mora, G.R.; Alagon, A.C.; King, G.; et al. Structural and molecular diversification of the Anguimorpha lizard mandibular venom gland system in the arboreal species Abronia graminea. J. Mol. Evol. 2012, 75, 168–183. [Google Scholar] [CrossRef]
  98. Jansa, S.A.; Voss, R.S. Adaptive evolution of the venom-targeted vWF protein in opossums that eat pitvipers. PLoS One 2011, 6, e20997. [Google Scholar] [CrossRef]
  99. Voss, R.S. Opossums (Mammalia: Didelphidae) in the diets of Neotropical pitvipers (Serpentes: Crotalinae): Evidence for alternative coevolutionary outcomes? Toxicon 2013, 66, 1–6. [Google Scholar] [CrossRef]
  100. Nisani, Z.; Boskovic, D.S.; Dunbar, S.G.; Kelln, W.; Hayes, W.K. Investigating the chemical profile of regenerated scorpion (Parabuthus transvaalicus) venom in relation to metabolic cost and toxicity. Toxicon 2012, 60, 315–323. [Google Scholar] [CrossRef]
  101. Nisani, Z.; Dunbar, S.G.; Hayes, W.K. Cost of venom regeneration in Parabuthus transvaalicus (Arachnida: Buthidae). Comp. Biochem. Physiol. Part A Mol. Integr. Physiol. 2007, 147, 509–513. [Google Scholar] [CrossRef]
  102. Morgenstern, D.; King, G.F. The venom optimization hypothesis revisited. Toxicon 2013, 63, 120–128. [Google Scholar] [CrossRef]
  103. Kini, R.M.; Chan, Y.M. Accelerated evolution and molecular surface of venom phospholipase A2 enzymes. J. Mol. Evol. 1999, 48, 125–132. [Google Scholar] [CrossRef]
  104. Delpietro, H.A.; Russo, R.G. Acquired resistance to saliva anticoagulants by prey previously fed upon by vampire bats (Desmodus rotundus): Evidence for immune response. J. Mamm. 2009, 90, 1132–1138. [Google Scholar] [CrossRef]
  105. Pillet, L.; Trémeau, O.; Ducancel, F.; Drevet, P.; Zinn-Justin, S.; Pinkasfeld, S.; Boulain, J.C.; Ménez, A. Genetic engineering of snake toxins. Role of invariant residues in the structural and functional properties of a curaremimetic toxin, as probed by site-directed mutagenesis. J. Biol. Chem. 1993, 268, 909–916. [Google Scholar]
  106. Antil, S.; Servent, D.; Menez, A. Variability among the sites by which curaremimetic toxins bind to torpedo acetylcholine receptor, as revealed by identification of the functional residues of alpha-cobratoxin. J. Biol. Chem. 1999, 274, 34851–34858. [Google Scholar] [CrossRef]
  107. Antil-Delbeke, S.; Gaillard, C.; Tamiya, T.; Corringer, P.J.; Changeux, J.P.; Servent, D. Ménez, a Molecular determinants by which a long chain toxin from snake venom interacts with the neuronal alpha 7-nicotinic acetylcholine receptor. J. Biol. Chem. 2000, 275, 29594–29601. [Google Scholar]
  108. Trémeau, O.; Lemaire, C.; Drevet, P.; Pinkasfeld, S.; Ducancel, F.; Boulain, J.C.; Ménez, A. Genetic engineering of snake toxins. The functional site of Erabutoxin a, as delineated by site-directed mutagenesis, includes variant residues. J. Biol. Chem. 1995, 270, 9362–9369. [Google Scholar]
  109. Heyborne, W.H.; Mackessy, S.P. Identification and characterization of a taxon-specific three-finger toxin from the venom of the Green Vinesnake (Oxybelis fulgidus; family Colubridae). Biochimie 2013, 95, 1923–1932. [Google Scholar] [CrossRef]
  110. Gotz, S.; Garcia-Gomez, J.M.; Terol, J.; Williams, T.D.; Nagaraj, S.H.; Nueda, M.J.; Robles, M.; Talon, M.; Dopazo, J.; Conesa, A. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008, 36, 3420–3435. [Google Scholar] [CrossRef]
  111. NCBI. Available online: www.ncbi.nlm.nih.gov (accessed on 14 November 2013).
  112. Altschul, S.F.; Madden, T.L.; Schaffer, A.A.; Zhang, J.; Zhang, Z.; Miller, W.; Lipman, D.J. Gapped BLAST and PSI-BLAST: A new generation of protein database search programs. Nucleic Acids Res. 1997, 25, 3389–3402. [Google Scholar] [CrossRef]
  113. Edgar, R.C. MUSCLE: Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004, 32, 1792–1797. [Google Scholar] [CrossRef]
  114. Ronquist, F.; Teslenko, M.; van der Mark, P.; Ayres, D.L.; Darling, A.; Hohna, S.; Larget, B.; Liu, L.; Suchard, M.A.; Huelsenbeck, J.P. MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 2012, 61, 539–542. [Google Scholar] [CrossRef]
  115. Whelan, S.; Goldman, N. A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach. Mol. Biol. Evol. 2001, 18, 691–699. [Google Scholar] [CrossRef]
  116. Guindon, S.; Dufayard, J.F.; Lefort, V.; Anisimova, M.; Hordijk, W.; Gascuel, O. New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0. Syst. Biol. 2010, 59, 307–321. [Google Scholar] [CrossRef]
  117. Posada, D.; Crandall, K.A. The effect of recombination on the accuracy of phylogeny estimation. J. Mol. Evol. 2002, 54, 396–402. [Google Scholar]
  118. Delport, W.; Poon, A.F.; Frost, S.D.; Kosakovsky Pond, S.L. Datamonkey 2010: A suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics 2010, 26, 2455–2457. [Google Scholar] [CrossRef]
  119. Kosakovsky Pond, S.L.; Posada, D.; Gravenor, M.B.; Woelk, C.H.; Frost, S.D.W.; Pond, S.L.K. Automated phylogenetic detection of recombination using a genetic algorithm. Mol. Biol. Evol. 2006, 23, 1891–1901. [Google Scholar] [CrossRef]
  120. Goldman, N.; Yang, Z. A codon-based model of nucleotide substitution for protein-coding DNA sequences. Mol. Biol. Evol. 1994, 11, 725–736. [Google Scholar]
  121. Yang, Z. Likelihood ratio tests for detecting positive selection and application to primate lysozyme evolution. Mol. Biol. Evol. 1998, 15, 568–573. [Google Scholar] [CrossRef]
  122. Yang, Z. PAML 4: Phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 2007, 24, 1586–1591. [Google Scholar] [CrossRef]
  123. Zhang, J.; Nielsen, R.; Yang, Z. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol. Biol. Evol. 2005, 22, 2472–2479. [Google Scholar] [CrossRef]
  124. Yang, Z.; Nielsen, R. Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol. Biol. Evol. 2002, 19, 908–917. [Google Scholar] [CrossRef]
  125. Nielsen, R.; Yang, Z. Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene. Genetics 1998, 148, 929–936. [Google Scholar]
  126. Yang, Z.; Wong, W.S.W.; Nielsen, R. Bayes empirical Bayes inference of amino acid sites under positive selection. Mol. Biol. Evol. 2005, 22, 1107–1118. [Google Scholar] [CrossRef]
  127. Murrell, B.; Moola, S.; Mabona, A.; Weighill, T.; Sheward, D.; Kosakovsky Pond, S.L.; Scheffler, K. FUBAR: A fast, unconstrained bayesian approximation for inferring selection. Mol. Biol. Evol. 2013, 30, 1196–1205. [Google Scholar] [CrossRef]
  128. Pond, S.L.K.; Frost, S.D.W.; Muse, S.V. HyPhy: Hypothesis testing using phylogenies. Bioinformatics 2005, 21, 676–679. [Google Scholar] [CrossRef]
  129. Murrell, B.; Wertheim, J.O.; Moola, S.; Weighill, T.; Scheffler, K.; Kosakovsky Pond, S.L. Detecting individual sites subject to episodic diversifying selection. PLoS Genet. 2012, 8, e1002764. [Google Scholar] [CrossRef]
  130. Woolley, S.; Johnson, J.; Smith, M.J.; Crandall, K.A.; McClellan, D.A. TreeSAAP: Selection on amino acid properties using phylogenetic trees. Bioinformatics 2003, 19, 671–672. [Google Scholar] [CrossRef]
  131. Bielawski, J.P.; Yang, Z. A maximum likelihood method for detecting functional divergence at individual codon sites, with application to gene family evolution. J. Mol. Evol. 2004, 59, 121–132. [Google Scholar]
  132. Pond, S.L.K.; Scheffler, K.; Gravenor, M.B.; Poon, A.F.Y.; Frost, S.D.W. Evolutionary fingerprinting of genes. Mol. Biol. Evol. 2010, 27, 520–536. [Google Scholar] [CrossRef]
  133. Kelley, L.A.; Sternberg, M.J.E. Protein structure prediction on the Web: A case study using the Phyre server. Nat. Protoc. 2009, 4, 363–371. [Google Scholar] [CrossRef]
  134. DeLano, W.L. The PyMOL Molecular Graphics System; DeLano Scientific: San Carlos, CA, USA, 2002. [Google Scholar]
  135. Armon, A.; Graur, D.; Ben-Tal, N. ConSurf: An algorithmic tool for the identification of functional regions in proteins by surface mapping of phylogenetic information. J. Mol. Biol. 2001, 307, 447–463. [Google Scholar] [CrossRef]
  136. Fraczkiewicz, R.; Braun, W. Exact and efficient analytical calculation of the accessible surface areas and their gradients for macromolecules. J. Comput. Chem. 1998, 19, 319–333. [Google Scholar] [CrossRef]

Supplementary Files

  • Supplementary File 1:

    Supplementary (ZIP, 3550 KB)

  • Share and Cite

    MDPI and ACS Style

    Sunagar, K.; Jackson, T.N.W.; Undheim, E.A.B.; Ali, S.A.; Antunes, A.; Fry, B.G. Three-Fingered RAVERs: Rapid Accumulation of Variations in Exposed Residues of Snake Venom Toxins. Toxins 2013, 5, 2172-2208. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins5112172

    AMA Style

    Sunagar K, Jackson TNW, Undheim EAB, Ali SA, Antunes A, Fry BG. Three-Fingered RAVERs: Rapid Accumulation of Variations in Exposed Residues of Snake Venom Toxins. Toxins. 2013; 5(11):2172-2208. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins5112172

    Chicago/Turabian Style

    Sunagar, Kartik, Timothy N. W. Jackson, Eivind A. B. Undheim, Syed. A. Ali, Agostinho Antunes, and Bryan G. Fry. 2013. "Three-Fingered RAVERs: Rapid Accumulation of Variations in Exposed Residues of Snake Venom Toxins" Toxins 5, no. 11: 2172-2208. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins5112172

    Article Metrics

    Back to TopTop