Abstract

Toads are chemically defended by bufadienolides, a class of cardiotonic steroids that exert toxic effects by binding to and disabling the Na+/K+-ATPases of cell membranes. Some predators, including a number of snakes, have evolved resistance to the toxic effects of bufadienolides and prey regularly on toads. Resistance in snakes to the acute effects of these toxins is conferred by at least two amino acid substitutions in the cardiotonic steroid binding pocket of the Na+/K+-ATPase. We surveyed 100 species of snakes from a broad phylogenetic range for the presence or absence of resistance-conferring mutations. We found that such mutations occur in a much wider range of taxa than previously believed. Although all sequenced species known to consume toads exhibited the resistance mutations, many of the species possessing the mutations do not feed on toads, much less specialize on that food source. This suggests that either there is little performance cost associated with these mutations or they provide an unknown benefit. Furthermore, the distribution of the mutation among major clades of advanced snakes suggests that the origin of the mutation reflects evolutionary retention more than dietary constraint.

1. Introduction

Many plants and animals are defended by toxic compounds, and circumvention of those defences has often involved the evolution of elaborate mechanisms for resistance to the toxins [1–5]. Among animals, some defensive compounds are synthesized from non-toxic precursors, whereas others are acquired and stored (sequestered) from environmental sources, primarily diet, and redeployed defensively [3,6–8]. Toads (Bufonidae) synthesize potent cardiotonic steroids known as bufadienolides (BDs) from cholesterol and store those toxins in high concentrations in their cutaneous glands [6]. Those toxins protect toads from the majority of predators, including many snakes that readily consume other species of frogs. Nonetheless, specialized bufophagy (the consumption of toads as a major (more than 80%) dietary component) or the sequestration of dietary BDs for defence occur in at least three genera of snakes on different continents (electronic supplementary material, table S1) [9–11]. Such systems of chemical defence and resistance between prey and predators exemplify the coevolutionary arms races that can drive molecular and physiological changes in the participating organisms [12].

Studies involving independently evolved systems of chemical defence and their circumvention by predators have shown that resistance to toxic dietary items can be accomplished by similar molecular changes that confer target-site insensitivity [11–15]. Furthermore, there is evidence that tight coupling between the resistance of predators and toxicity of prey in a coevolutionary arms race can result when the predator incurs a performance cost from resistance [12]. When a toxin targets a critical cellular function, phenotypes resistant to that toxin sometimes evolve through parallel or convergent mutations of the same gene across multiple species or populations [11–13,16–18].

BDs and other cardiotonic steroids, including the plant-derived cardenolides (CDs) [19,20], target the ubiquitous cell-membrane enzyme Na+/K+-ATPase, inhibiting the transport of sodium ions out of the cell and potassium ions in. The Na+/K+-ATPase molecule consists of a catalytic α subunit, with 10 transmembrane segments, which carries out the functional activities of the enzyme, and a glycosylated β subunit, with one transmembrane segment, which provides structural stability (figure 1a) [22–25]. The H1–H2 extracellular loop of the α subunit comprises a significant portion of the binding pocket for cardiotonic steroids (figure 1b). However, some sites on the H3–H4 and the H5–H6 extracellular loops also contribute to the binding pocket to a lesser extent [13], and recent evidence suggests that mutations near the C-terminal end also confer resistance [26].

(a) Schematic diagram of an unfolded Na+/K+-ATPase molecule, showing the location of the Q111 L and G120R substitutions on the H1–H2 extracellular loop. The box shows the area enlarged in b. (b) Enlargement of the H1–H2 extracellular loop and the first two transmembrane domains. (c) The amino acid sequence of the first extracellular loop, indicating the substitutions at positions 111 and 120 that confer resistance to bufadienolides in snakes. Redrawn from Köksoy [21]. For a more detailed structural reconstruction of the interactions of the Na+/K+-ATPase and a BD, see Laursen et al. [22]. (Online version in colour.)

Several studies have found that two amino acid substitutions, one at position 111 and the other at or near position 122 of the H1–H2 loop, confer lowered affinity for cardiotonic steroids in resistant insects, lizards, and snakes [11,13–15,27–30]. Dobler et al. [13] identified various amino acid replacements encoded by the gene for Na+/K+-ATPase, including replacements at two positions in the H1–H2 extracellular loop of several independent lineages of CD-resistant insects. Ujvari et al. [15] identified two similar amino acid substitutions in the H1–H2 extracellular loop of BD-resistant species of varanid lizards: a leucine (L) replaces a glutamine (Q) at position 111 (Q111 L) and an arginine (R) replaces a glycine (G) at position 120 (G120R). Together, the Q111 L and G120R substitutions were shown to confer significantly reduced affinity of Na+/K+-ATPase for BDs when transfected into human kidney cells [15]. These two substitutions were subsequently found to occur in several species of snakes that are resistant to BDs [11]. Mohammadi et al. [29] investigated the extent of resistance conferred at the whole-body level in a species of snake known to possess the Q111 L and G120R mutations but feeds only rarely on toads. Those snakes were able to withstand concentrations of intraperitoneally injected BDs 200-fold higher than the mouse LD50.

Site-directed mutagenesis of the gene that encodes for Na+/K+-ATPase in Drosophila melanogaster has revealed that the Q111 L substitution alone, in CD-resistant insects, contributes slightly to resistance, whereas the mutation at position 122 (N122H, with histidine replacing asparagine) contributes significantly [13]. The N122H substitution results in a substantial biochemical change, the replacement of a neutral amino acid by a positive one [13], so it is no surprise that this mutation significantly alters the sensitivity of insect Na+/K+-ATPases to CDs. Furthermore, because the G120R mutation in BD-resistant reptiles results in the same biochemical change two amino acid positions away from the N122H mutation of insects [11,15,29], it is safe to assume that the G120R mutation is likely to be the major contributor to lowered BD affinity in mutant reptile Na+/K+-ATPases [15].

Adaptation to similar dietary toxins through the evolution of virtually identical amino acid substitutions across a range of taxa, spanning invertebrates and vertebrates, represents a striking example of convergent evolution. Previous studies have typically concluded that the resistance-conferring mutations co-occur with dietary specialization [11–14], suggesting that resistance to dietary toxins occurs only under strong selective pressure.

A recent survey of genetic resistance to cardiotonic steroids across the Metazoa [11] noted a strong correspondence in snakes between Q111 L and G120R mutations and predation on toads. Ujvari et al. [11] recovered three origins of resistance in their sample of 27 snakes and associated each origin with the assumption of toad-eating habits. Here, we expand the sample to 100 species of snakes, finding many more species that possess the resistance-conferring Q111 L and G120R substitutions. We reconstruct the ancestral states for these taxa and find, in contrast to the more limited taxon sampling of Ujvari et al. [11], that the mutations for resistance occur more widely among snakes than does either specialized or opportunistic bufophagy. Indeed, the mutations occur in some species with diets restricted to prey that lack cardiotonic steroids. This lack of close correspondence between the mutation and toad consumption suggests that, in contrast to similar studies of predators on amphibians defended by tetrodotoxin (which blocks sodium channels) [12], resistance to BDs does not impose a significant performance cost or there is some other, unknown benefit. Furthermore, our results suggest that target-site insensitivity is necessary, but may not be sufficient, to permit dietary specialization on toads or the sequestration of dietary BDs.

2. Material and methods

(a) Selection and classification of study taxa

Because our goal was to test whether a mutated Na+/K+-ATPase alone could explain the evolution of specialized bufophagous habits, we began by sequencing part of the ATP1a3 gene in previously unsequenced toad-eating species. Our criteria for toad-eating habits are restrictive and do not include diet data from general and secondary sources, which classify some species based upon anecdotal dietary records or assumptions of diet based on the predator's habitat. We based our dietary classification (electronic supplementary material, table S1) almost exclusively on dietary studies and other published reports in the primary literature. However, we also scored as (potentially) bufophagous those species that possess the resistance mutations and have a broad diet of amphibians, even in the absence of specific reports of toads in their diet. This conservative assumption maximizes the perceived co-occurrence of the mutations and bufophagy.

We also broadened our genetic sampling to include non-toad-eating sister taxa of the specialized and facultative toad-eaters, to determine whether the distribution of mutations for BD resistance closely correspond to toad-eating habits. Combined with other sequence data available in GenBank, we included data for 100 species of snakes (see electronic supplementary material). Although not shown in the figures, our phylogenetic analysis included 21 lizards (18 of them in the genus Varanus, plus one dactyloid, one agamid, and one scincid), as well as Sphenodon, Alligator, Gallus, and Xenopus, as outgroups.

(b) Sequence data

Ethanol-preserved tissue samples were obtained from the Louisiana State University Museum of Zoology. Other tissue samples were acquired from snakes collected in the field in Louisiana (permit no. LNHP-12-064), Utah (permit no. 1COLL9134), and Nebraska (permit no. 477) or from previously frozen and ethanol-preserved specimens in our laboratories at Utah State University and Kyoto University. DNA was isolated from tissues, using either the DNeasy Blood and Tissue Kit (Qiagen™) or the Phire direct polymerase chain reaction (PCR) kit (Thermo Fisher™). The H1–H2 extracellular loop of the ATP1a3 gene was amplified, using the primers (Integrated DNA Technologies®) ATP1a3Fwd (CGAGATGGCCCCAATGCTCTCA) and ATP1a3Rvs (TGGTAGTAGGAGAAGCAGCCGGT), which were modified from Moore et al. [31]. Amplified DNA product was purified, using the Wizard SV Gel and PCR clean-up system (Promega™). DNA was sequenced on an ABI 3730xo Genetic Analyzer, using BigDye Terminator kit v. 3.1. Sequences were aligned and analysed, using Geneious v. 6.1.8 (Biomatters Ltd.). Sequences were translated to amino acids, and the hydrophobicity and isoelectric points of each residue were calculated, using Geneious v. 6.1.8 (Biomatters Ltd.). Some sequences were obtained from GenBank (accession numbers KP238131-KP238176, XM_007437572, and AZIM01015831). DNA sequence data generated in this study have been submitted to GenBank with the accession numbers KU738063-KU738139, KU933521, and KU950321.

(c) Ancestral state reconstruction

Stochastic character mapping [32] was used to obtain probabilistic reconstructions of ancestral states for non-synonymous substitutions at amino acid positions 111 (second codon position, A versus T) and 120 (first codon position, A versus G or C; see Results). Character mapping was based on the time-calibrated squamate reptile phylogeny from Zheng & Weins [33], which we pruned to retain only the 121 species for which we had ATP1a3 sequences. Sequences were available for seven species that were not included in Zheng and Wiens' phylogeny, although congeners were included in that tree. Therefore, we considered the sequenced species to be equivalent to their congeners for the purposes of the analysis, placing them at the same point on the phylogeny. The included species were Conophis vittatus (replacing C. lineatus), Rhabdophis leonardi (replacing R. nuchalis), Cyclophiops semicarinatus (replacing C. major), Dendrelaphis subocularis (replacing D. tristis), D. punctulatus (replacing D. bifrenalis), Causus maculatus (replacing C. rhombeatus), and Vipera aspis (replacing V. ammodytes).

We first obtained the joint posterior distribution of ancestral character states and the transition matrix for each nucleotide position. We allowed all transition rates to differ and assumed a flat prior on the root states. Posteriors were inferred using Markov chain Monte Carlo (MCMC) with a 10 000 iteration burn-in, 100 000 sampling iterations, and a thinning interval of 50 iterations. Two chains were run to evaluate convergence of the MCMC algorithm to the stationary distribution, and results were combined for inference. Stochastic character maps were simulated over the posterior distribution of the transition matrix and thus account for uncertainty in transition rates. This analysis was conducted using the make.simmap function in the R package phytools [34].

In contrast to CD-resistant insects, in which position 122 of Na+/K+-ATPase has a positive amino acid, we always observed a shift from a neutral to positive amino acid at position 120 in BD-resistant snakes. Therefore, we asked whether such a pattern could be explained by chance under the null model that both changes would have an equivalent effect and are equally likely. Our calculations were made over the posterior distribution of the character state reconstructions, thus accounting for uncertainty in the number of changes.

We next compared models of correlated versus independent evolution for the two amino acids. We fitted continuous-time Markov models for trait evolution, either with separate transition matrixes for each amino acid (independent evolution model) or with a single transition matrix for the pair of amino acids (dependent or correlated, evolution model) [34]. We obtained maximum-likelihood parameter estimates, using BayesTraitsV2 (25 starts were used for parameter estimation). Models were then compared, using a likelihood ratio test as described by Pagel [35].

3. Results

(a) Amino acid substitutions

Our analysis of Na+/K+-ATPase α3 subunit sequences for the H1–H2 extracellular loop revealed substitutions at several residues. In our survey of 100 species of snakes, we found two consistent amino acid substitutions in all taxa known to consume toads. At position 111, we found the Q111 L substitution and at position 120, the G120R substitution, a result that is consistent with previously published findings (figure 1c) [11]. These two substitutions always occurred together and usually involved the same nucleotide changes. Notably, in addition to the occurrence of the Q111 L and G120R substitutions in all bufophagous species included in our analysis, we found many instances of these same two substitutions (hereinafter, ‘resistance mutations’) across the phylogeny of snakes, including many species that rarely or never feed on toads. In addition to the Q111 L and G120R substitutions, several other substitutions occurred at various positions among different groups (electronic supplementary material, figure S1). However, these changes are not known to alter the biochemical properties of the binding pocket to any significant extent, and thus are believed to be of no major functional consequence [11].

(b) Ancestral state reconstruction

Stochastic character mapping supported multiple gains and losses of the resistance mutations across the tree (figure 2). We found evidence of 11 (95% credible intervals [CIs] 2–18) transitions from A to T (with eight reversals, 95% CIs 2–19) at position 2 of amino acid 111 and 11 (95% CIs 1–33) transitions from A to G (with 10 reversals, 95% CIs 1–32) at position 1 of amino acid 120. Given the total number of transitions and the paucity of taxa with only one of the resistance mutations, we can confidently reject the null model that changes at amino acid 111 and 120 are equivalent (binomial probability: median across MCMC replicates, p = 0.0005; mean across MCMC replicates = 0.0511). Moreover, a model of correlated or dependent evolution for amino acid positions 111 and 120 was preferred relative to a null model of independent evolution (likelihood ratio = 88.85, χ2 with 4 d.f., p < 0.001).

Ancestral state reconstruction of the mutations to Na+/K+-ATPase that confer resistance to bufadienolides in snakes. Coloured circles beside terminal taxa indicate the amino acid composition of Na+/K+-ATPase at the H1–H2 extracellular loop (red = resistant L111 and R120, and black = non-resistant Q111 and G120). The names of bufophagous species are shown in red, and non-bufophagous species are shown in black (electronic supplementary material table S1) [10,35,36]. Inferred ancestral states are denoted by pie charts that indicate the posterior probability of having (red) versus not having (black) the resistance mutations. Selected higher taxa are identified, and circled numbers identify nodes that are discussed in the text. Note that the resistance mutations to Na+/K+-ATPase are not limited to taxa that prey on toads (as shown by the number of black species names beside red circles). Two highly bufophagous species (Heterodon platirhinos and Xenodon merremii) and one that sequesters BDs (Rhabdophis tigrinus), indicated by green arrows, possess enlarged adrenal glands that may function in countering the effects of chronic exposure to bufadienolides. Inset: Heterodon platirhinos feeding on the toad Anaxyrus terrestris. Note the viscous skin secretion of the toad, which is rich in bufadienolides.

4. Discussion

Our results demonstrate that the mutations at positions 111 and 120 in the H1–H2 extracellular loop of the Na+/K+-ATPase α3 subunit, which confer resistance to toad toxins in snakes [11,15,29], occur in a much wider range of species than previously believed, including many species that do not regularly, if ever, prey on toads (electronic supplementary material, table S1). The first extracellular loop (amino acid positions 111–122) comprises a major part of the cardiotonic steroid binding pocket, and the non-resistant gene codes for the amino acids glutamine at position 111 and glycine at position 120. The observed substitution of leucine at position 111 results in a biochemical residue change from polar-neutral to non-polar-neutral, which increases the hydrophobicity of the residue but does not change the isoelectric point (electronic supplementary material, figure S1). The substitution of arginine at position 120 results in a biochemical residue change from non-polar-neutral to polar-positive, which increases the isoelectric point and decreases the hydrophobicity of the residue (electronic supplementary material, figure S1). Although the physiological effect of the G120R substitution has not been specifically tested, it is biochemically similar to the residue change caused by the N122H substitution in milkweed insects (change from a neutral to positively charged residue), which has been shown to confer significantly lower binding affinity to CDs [13]. Furthermore, the combined substitutions of Q111 L and G120R have been shown to confer a much lower affinity to BDs than the wild-type gene in vertebrate cells [15]. Although the Q111 L substitution does contribute to resistance in cardenolide-resistant invertebrates [13], evidence from insect phylogenies suggests that it is a necessary intermediate step that leads to the later appearance of other substitutions in the H1–H2 extracellular loop [4]. A similar relationship presumably explains the consistent co-occurrence of substitution Q111 L with G120R in BD-resistant snakes.

The distribution of the resistance mutations Q111 L and G120R among colubroid snakes, including their widespread appearance in many species that do not prey on toads, raises a number of important issues that were obscured by the limited taxon sampling of previous studies. The occurrence of such mutations among non-toad-eating species suggests that possession of these mutations does not carry a high performance cost. Paradoxically, however, investigations of the ion transport activities of some resistant Na+/K+-ATPases have shown that the mutations do in fact lower the enzyme's function [36,37]. It appears that in insects, the deleterious effects of the lowered enzyme function are circumvented by the presence of a secondary Na+/K+-ATPase isoform [37]. It is thus likely that in snakes Na+/K+-ATPase isoforms α1 and α2 might be contributing similarly. This contrasts sharply with the limited distribution of mutations for resistance in snakes that consume amphibians defended by tetrodotoxin (TTX) [12]. TTX targets voltage-gated sodium channels, and the mutations that confer resistance to TTX also impart a substantial countervailing cost, as measured by locomotor performance [38,39]. Several point mutations have been documented among TTX-resistant snakes [12], and their occurrence is tightly linked to the consumption of TTX-defended amphibians. In the best-documented case, involving the natricine snake Thamnophis sirtalis and the newt Taricha granulosa, geographical differences in the level of TTX in the newt are mirrored by differences in the level of innate resistance of the sodium channels in the snake [40]. Populations of snakes that consume low levels of TTX exhibit low levels of resistance.

In contrast to such close correspondence between prey toxicity and sodium channel resistance to TTX, our data document the occurrence of BD-resistance mutations in Na+/K+-ATPase among many species that are not known to consume toads at all, including specialized predators on molluscs (Dipsas catesbyi, Storeria occipitomaculata), crayfishes (Regina spp.), fishes (Nerodia rhombifer, Helicops angulatus), and mammals (Crotalus atrox, Agkistrodon contortrix; electronic supplementary material, table S1). It is possible that the same mutation protects a few of these taxa against other toxins that exert yet unknown pharmacological effects (as in some molluscs [41] and other invertebrates [42–46]). However, most non-toad-eating species that possess the resistance mutations appear not to derive any benefit nor incur a substantial cost from the mutated form of Na+/K+-ATPase. Future studies should address that inference directly.

Our study includes relatively broad taxon sampling of the three largest subfamilies of colubrid snakes. The sister lineages Natricinae and Dipsadinae (figure 2) both include large numbers of taxa that possess the resistance mutations. The species in both lineages consume a wide variety of prey, and both subfamilies include species that are highly bufophagous or sequester BDs. Our results indicate a 0.61 probability that resistance to BDs originated in the common ancestor of both of these groups (figure 2, node 3) and was independently lost in two extant lineages. Interestingly, both these subfamilies include early branches that are highly bufophagous or sequester BDs, suggesting the possibility that consumption of toads characterized early members of both radiations, and possibly their common ancestor. Arguably the most highly bufophagous genus of snakes, Heterodon of North America, lies near the base of the Dipsadinae, although members of the sister lineage of Heterodon (which includes Farancia, Diadophis, Carphophis, and Contia) lack the resistance mutations, most likely through secondary loss of those residues (figure 2, node 2). The other highly bufophagous lineage of dipsadines, the genus Xenodon (sensu lato, including the nominal genera Waglerophis and Lystrophis [47]), is nested more deeply within the Dipsadinae, as a member of a diverse South American radiation (figure 2, node 1).

Among the Natricinae, the overwhelming majority of species possesses the resistance mutations. A basal split (figure 2, node 5) separates an exclusively Asian lineage from one that includes a few Asian species, as well as a number in Europe and the large North American tribe Thamnophiini. The exclusively Asian lineage includes the genus Rhabdophis, most members of which possess defensive structures known as nuchal or nucho-dorsal glands [48]. Some species accumulate high concentrations of BDs from toads consumed as prey and release the toxins to deter the snakes' predators [3,49]. Furthermore, female Rhabdophis tigrinus can transfer prey toxins to their embryos before oviposition [50], and a recent field study demonstrated that females preferentially forage for toads during the reproductive season, presumably increasing the likelihood of provisioning their embryos with BDs [51]. Therefore, although toads comprise a numerically small element of the diet in this species [52], the snakes often retain high levels of BDs in their tissues. The sampled species of Natricinae exhibit only a single example of likely secondary loss of the resistance mutations, the East Asian genus Sinonatrix (figure 2, node 4). That genus of highly aquatic fish-eaters (electronic supplementary material, table S1), is sister to the West Asian, European, and North American natricines, all of which possess the resistance mutations and one of which, Natrix natrix of Europe, frequently consumes toads.

The third broadly sampled colubrid subfamily, the Colubrinae, displays a different evolutionary picture. Of 15 species sampled, only two distantly related species possess the resistance mutation (and are known to feed on toads). Our analysis suggests a 91% probability that the ancestor of the Colubrinae (figure 2, node 6) lacked the mutation. Thus, although all species in our sample that are known to prey on toads possess the resistance mutations, the recovered pattern is one in which only a few members of the ancestrally non-resistant Colubrinae evolved resistance, whereas many members of the ancestrally resistant Natricinae and Dipsadinae retained the mutations despite consuming prey other than toads. These results argue that historical contingency plays a strong role in the origin of genetic resistance to BDs and suggest that the performance cost of retaining the mutations is low, even when it confers no evident advantage.

Two other large and widespread families, Viperidae and Elapidae, possess bufophagous genera, both of which occur in Africa. Viperids exhibit greater variation in the distribution of the resistance mutations than other major lineages. However, the bufophagous viper Causus is basal to a radiation of African Viperinae (figure 2, node 8), of which all sequenced taxa retain those mutations. The Eurasian genus Vipera diverged before the African radiation and lacks the mutations. The other major lineage of vipers, the Crotalinae, exhibits a diversity of Na+/K+-ATPase sequences. However, our sampling of Viperidae is limited, and further sequencing is needed to resolve the origins of resistance in this family with confidence. Among the Elapidae, the bufophagous African genus Hemachatus is sister to the diverse radiation of African and Asian Naja (figure 2, node 7), which includes a number of facultative toad-eaters. Unfortunately, our sampling includes only a single species of Lamprophiidae, a large and ecologically diverse African colubroid radiation that is sister to the Elapidae. The only sequenced lamprophiid, Boaedon fuliginosus, lacks the resistance mutations, but further sampling of this family is desirable, especially considering that one species of the lamprophiid genus Crotaphopeltis feeds largely on toads [9].

It is possible that some generally non-toad-eating taxa may rarely encounter and consume bufonids, either as adults or as juveniles. (Information on the diets of young snakes is especially limited.) However, it is doubtful that rare instances of toad consumption would have driven the evolution of the resistance mutations, especially in the varied dietary specialists discussed above. It is more likely that species belonging to ancestrally resistant lineages may occasionally consume bufonids when encountered, despite a diet consisting primarily of other prey. Such a scenario is not inconsistent with our conclusions. Conversely, however, we predict that non-resistant snakes would rarely, if ever, survive such a rare predatory event.

Recent time-calibrated phylogenetic studies of toads and snakes support the hypothesis that toads arrived on all occupied continents prior to, or coincident with, the evolution of specialized bufophagous snakes. The Bufonidae are believed to have arisen in South America, after the breakup of Gondwana approximately 88 Ma, and dispersed north and then west across the Bering Strait, establishing lineages in Eurasia and Africa before returning to the Americas across a North Atlantic land connection about 43 Ma [53]. The presence of BDs in the most basal genus of bufonids, Melanophryniscus, is uncertain, although its species are chemically defended by the indolealkylamine bufotenine and by alkaloids sequestered from arthropod prey [54]. Skin extracts from Melanophryniscus appear to inhibit Na+/K+-ATPase [55]. BDs are present in Atelopus [55], which is sister to all remaining bufonids in both the New and Old Worlds. Furthermore, a more derived clade of Neotropical bufonids (Rhinella) diverged from a North American ancestor about 41 Ma, so the presence of BD-defended toads in South America clearly predates the origin of the bufophagous South American genus Xenodon about 18 Ma. Similarly, the origin of the notably bufophagous taxa Heterodon (North America, 30 Ma), Rhabdophis (East Asia, 29 Ma), Natrix natrix (Europe, 13 Ma), Hemachatus (Africa, 23 Ma), and Causus (Africa, 30 Ma) all postdate the completion of the bufonid diaspora [33].

Importantly, dates vary for the divergence of the related colubrid subfamilies Natricinae and Dipsadinae, both of which include relatively basal bufophagous genera and large numbers of genera that possess the resistance mutations. Indeed, our analysis indicated a fairly high probability that the ancestor of Natricinae + Dipsadinae (figure 2, node 3) possessed those mutations. A detailed analysis of the Natricinae suggested that they originated in eastern Asia approximately 35 Ma (figure 2, node 5), later establishing lineages in Africa, Europe, and North America [56]. Consistent with that scenario, the most basal dipsadines, Thermophis and Stichophanes, occur in Asia, whereas all other members of the subfamily are limited to the New World, including the highly bufophagous Heterodon and Xenodon. More recently, Zheng & Wiens [33] dated the divergence of Natricinae and Dipsadinae to about 48 Ma. These two divergence estimates bracket the date for the re-entry of toads into the New World from Eurasia (about 43 Ma; [38]) and are reasonably consistent with an origin of the mutations in the ancestor of Natricinae + Dipsadinae in Asia while in sympatry with bufonids.

The earliest lineage of snakes in which we found the resistance mutations is Epictia albipuncta (Leptotyphlopidae; figure 2, node 9), our only sample from the large but poorly understood radiation of Scolecophidia. That possibly paraphyletic group of five families and about 400 species [33,57] underwent an initial split into three major lineages about 123–125 Ma [57], long before the origin of toads. Given the presence of the resistance mutations in some varanid lizards [11,15], it is possible that the mutations may have been present in the common ancestor of all snakes. Alternatively, resistance may have evolved independently in scolecophidians, nearly all of which specialize on a diet of social insects, either ants or termites [8]. Ants are defended by a diverse array of alkaloids [42,43], and some termites possess terpene toxins [43–46], although the diversity of toxins in the early life-history stages on which scolecophidians largely feed is not well known. The pharmacological effects of most of ant and termite compounds are also poorly understood, and thus it is uncertain whether the same mutations to Na+/K+-ATPase that confer resistance to BDs would protect against those prey. However, evidence of resistance to both BDs and alkaloids has been demonstrated in Heterodon platirhinos, a highly bufophagous species that is also resistant to TTX in newts, although the mechanism for its resistance to TTX remains unknown [58]. Clearly, additional sequencing of diverse scolecophidian species is needed, together with a better understanding of the physiological effects of their prey's defensive toxins.

Although our sampling efforts were biased toward certain lineages, our data are sufficient to demonstrate that the presence of BD resistance mutations across the phylogeny of snakes does not correspond closely to the consumption of toads. Nonetheless, it is true that all bufophagous snakes in our sample possess the resistance mutations. Comparison of the physiological effects of a controlled dose of the BD gamabufotalin on a species with the resistance mutations, Nerodia rhombifer, and one that lacks the mutations, Pituophis catenifer, revealed sharply different responses [29]. N. rhombifer, a specialized predator on fishes, exhibited no signs of distress when subjected to mass-specific doses 20 times greater than those that proved lethal to P. catenifer. This suggests that the resistance mutations are indeed protective against the acute effects of BD consumption. Furthermore, the widespread occurrence of the resistance mutations among non-toad-eating snakes suggests that possession of the resistance mutations does not impose a large performance cost. However, the inclusion of toads as an occasional prey item among the genetically resistant species, together with the relative rarity of species that either specialize on a diet of toads or sequester dietary BDs, suggests that the cost of chronic exposure to those toxins poses a greater physiological challenge and may require additional mechanisms of resistance.

In parallel invertebrate systems, a variety of mechanisms confer resistance to CDs, including detoxification by enzymes, excretion and exclusion from sensitive tissues, in addition to target-site insensitivity [1,2,4,5]. As such, target-site insensitivity may not always be required as a precursor for tolerance to toxic prey. Furthermore, several resistance mechanisms against environmental toxins can exist within a single species. For example, dichlorodiphenyltrichloroethane resistance in insects can be achieved by amino acid substitutions in the voltage-gated sodium channel [17,59] and also through increased expression of detoxifying enzymes [17,60]. This further suggests that, although the resistance mutations of Na+/K+-ATPase may be sufficient to protect a snake from the acute effects of BD ingestion, additional physiological adaptations may be required to protect highly bufophagous taxa and those that sequester BDs against the chronic effects of toxin exposure. Some CD-resistant insects show tissue-specific increases of mutant Na+/K+-ATPase expression that reflect the relative exposure of different tissues to CDs [14]. A similar mechanism in bufophagous snakes may involve a specific endocrine response to BDs. Mohammadi et al. [10] demonstrated that three independently evolved snakes that specialize on toads possess highly enlarged adrenal glands in comparison with their non-toad-eating sister taxa (figure 2). Adrenal corticosteroid hormones exert diverse physiological effects, and aldosterone in particular has been linked to increased expression of Na+/K+-ATPase in cardiac cells [61].

5. Conclusion

Although target-site insensitivity clearly plays an essential role in the survival of toad-eating snakes, it is not tightly coupled to bufophagy. The occurrence of the resistance mutations among species with varied diets suggests that the mutations impose little performance cost. To understand the evolution of highly specialized bufophagous snakes, it will be necessary to move beyond target-site insensitivity and investigate the complex, interconnected physiological pathways that confer resistance to chronic BD exposure.

Authors' contributions

S.M. and A.H.S. designed the research; S.M., J.G., and H.T. performed the research and collected data; Z.G., S.M., and A.H.S. analysed the data; A.M. contributed to the intellectual development of the research; S.M., A.H.S., Z.G., and J.G. wrote the paper.

Competing interests

We declare we have no competing interests.

Funding

This study was supported by a SICB Grant-in-Aid of Research (to S.M.), a USU Ecology Center Research Grant (to S.M.), a USU Center for Women and Gender Research Grant (to S.M.), the National Science Foundation (EAPSI SP15039 to S.M.), and a grant from the Japan Society for the Promotion of Science (Scientific Research C, no. 26440213 to A.M.).

Acknowledgements

We thank Caleb D. McMahan, Bryan M. Kluever, Andrew M. Durso, and Lorin A. Neuman-Lee for assistance in collecting snakes in the field. We especially thank the Louisiana State University Museum of Natural Science for providing the tissues that made this study possible.

. 2016Mutations to the cardiotonic steroid binding site of Na+/K+-ATPase are associated with high level of resistance to gamabufotalin in a natricine snake. Toxicon114, 13–15. (doi:10.1016/j.toxicon.2016.02.019)

. 2002The evolutionary response of predators to dangerous prey: hotspots and coldspots in the geographic mosaic of coevolution between newts and snakes. Evolution56, 2067–2082. (doi:10.1111/j.0014-3820.2002.tb00132.x)

. 2009Molecular phylogeny of advanced snakes (Serpentes: Caenophidia) with an emphasis on South American Xenodontines: a revised classification and descriptions of new taxa. Pap. Avulsos. Zool. (São Paulo)49, 115–153. (doi:10.1590/S0031-10492009001100001)