Royal Society Publishing

Genetic architecture of a feeding adaptation: garter snake (Thamnophis) resistance to tetrodotoxin bearing prey

Chris R. Feldman, Edmund D. Brodie, Edmund D. Brodie, Michael E. Pfrender


Detailing the genetic basis of adaptive variation in natural populations is a first step towards understanding the process of adaptive evolution, yet few ecologically relevant traits have been characterized at the genetic level in wild populations. Traits that mediate coevolutionary interactions between species are ideal for studying adaptation because of the intensity of selection and the well-characterized ecological context. We have previously described the ecological context, evolutionary history and partial genetic basis of tetrodotoxin (TTX) resistance in garter snakes (Thamnophis). Derived mutations in a voltage-gated sodium channel gene (Nav1.4) in three garter snake species are associated with resistance to TTX, the lethal neurotoxin found in their newt prey (Taricha). Here we evaluate the contribution of Nav1.4 alleles to TTX resistance in two of those species from central coastal California. We measured the phenotypes (TTX resistance) and genotypes (Nav1.4 and microsatellites) in a local sample of Thamnophis atratus and Thamnophis sirtalis. Allelic variation in Nav1.4 explains 23 per cent of the variation in TTX resistance in T. atratus while variation in a haphazard sample of the genome (neutral microsatellite markers) shows no association with the phenotype. Similarly, allelic variation in Nav1.4 correlates almost perfectly with TTX resistance in T. sirtalis, but neutral variation does not. These strong correlations suggest that Nav1.4 is a major effect locus. The simple genetic architecture of TTX resistance in garter snakes may significantly impact the dynamics of phenotypic coevolution. Fixation of a few alleles of major effect in some garter snake populations may have led to the evolution of extreme phenotypes and an ‘escape’ from the arms race with newts.

1. Introduction

Since the early days of the modern synthesis, genetic studies combined with the analysis of phenotypic variation have been fundamental to our understanding of the process of adaptive evolution (Fisher 1958; Futuyma 1998; Orr 2005). In recent years, the study of adaptive variation has benefited from major gains in molecular genetic techniques and a growing body of work is now documenting the molecular basis of adaptation (e.g. Abzhanov et al. 2004; Albertson et al. 2005; Geffeney et al. 2005; Joron et al. 2006; Storz et al. 2007; Feldman et al. 2009). The results of these studies are rapidly increasing our understanding of the process and constraints of adaptive evolution. Nevertheless, a number of unresolved issues remain regarding the generality of the mechanisms suggested by empirical work (Orr 2005). Chief among these is whether adaptive evolution is typically driven by changes in many genes of minor effect (micromutational or polygenic model), or mutations in a few genes of major effect (macromutational or oligogenic model; Orr & Coyne 1992; Barton & Keightley 2002; Orr 2005).

Historically, it has been thought that adaptive evolution is primarily fuelled by changes in many genes with minor effects (Fisher 1958; Lande 1983; Falconer & Mackay 1996). However, it has become increasingly clear that much of the variation in complex phenotypic traits is controlled by a small number of loci of major effect (e.g. Doebley & Stec 1991; Cohn & Tickle 1999; Krieger & Ross 2002; Abzhanov et al. 2004; Albertson et al. 2005; Hoekstra et al. 2006; Sutter et al. 2007; Storz et al. 2009), challenging the polygenic model of adaptive evolution (Orr & Coyne 1992; Orr 2005). While our view of the genetic architecture of adaptations is rapidly changing, both the polygenic and oligogenic viewpoints are certainly valid (Orr & Coyne 1992; Barton & Keightley 2002; Orr 2005). Each genetic model makes specific predictions about the tempo and mode of evolution, as well as the nature of selection (Fisher 1958; Lande 1983; Barton & Keightley 2002; Orr 2005). When the genetic architecture of a trait under selection is relatively simple, dramatic adaptation can occur rapidly through the fixation of a few beneficial mutations of large effect (Orr 2005; Barrett & Schluter 2008; e.g. Hawthorne & Via 2001; ffrench-Constant et al. 2004; Colosimo et al. 2005; Hartley et al. 2006). Yet phenotypic evolution will proceed more slowly and through incremental steps if adaptive traits are controlled by large numbers of loci (Orr 2005). Polygenic adaptations will require many mutations, most of which are expected to be recessive and take longer to spread to fixation than dominant alleles (Kimura 1983). Furthermore, some loci might back-mutate and accumulate disadvantageous recombinants while waiting on adaptive mutations, whereas other loci may be constrained by epistatic interactions or pleiotropy. Thus, adaptive evolution should proceed gradually unless multiple loci can change in concert.

The genetic architecture of adaptations is also an important component of the dynamics of coevolution (Bohannan & Lenski 2000; Thompson 2005; Wade 2007; e.g. Yoshida et al. 2007). Cyclical coevolution wherein predator (exploiter) phenotypes continually ‘chase’ prey (victim) phenotypes are expected regardless of whether the genetic architecture underlying the traits is simple or complex (Sasaki 2000; Agrawal & Lively 2002, 2003; Kopp & Gavrilets 2006; Nuismer et al. 2007). However, when the traits that mediate the phenotypic interface of coevolution have a simple genetic basis, reciprocal selection can more easily fix adaptive alleles, leading to the evolution of extreme phenotypes in which the one species ‘loses’ the arm-race (Sasaki 2000; Kopp & Gavrilets 2006; Nuismer et al. 2007). Conversely, coevolutionary cycles are predicted to be more stable over time and more dynamic under a multilocus model because adaptive alleles in the predator rarely change in harmony, allowing the one species to stay ‘ahead’ in the arms race (Kopp & Gavrilets 2006; Nuismer et al. 2007).

The interaction between toxic newts (Taricha) and several resistant predatory garter snakes (Thamnophis) provides a model system for the study of adaptive variation and predator–prey coevolution (e.g. Brodie & Brodie 1999; Brodie et al. 2002a; Hanifin et al. 2008). This system is ideal because the traits that mediate coevolution are identified, geographically variable, and at least partly controlled by a well-studied gene family. Newts of the genus Taricha possess the neurotoxin tetrodotoxin (TTX; Mosher et al. 1964; Wakely et al. 1966; Brodie et al. 1974; Yotsu et al. 1990), which acts as a powerful chemical defence against vertebrate predators (Brodie 1968; Brodie et al. 1974). Tetrodotoxin binds selectively to the outer pore of voltage-gated sodium channels in nerves and muscles (Lipkind & Fozzard 2000; Hille 2001), blocking the movement of sodium ions (Na+) across the cell membrane and halting the propagation of action potentials (Kao & Levinson 1986; Hille 2001). By arresting nerve impulses in muscle and nervous tissue, TTX causes immobilization, respiratory failure and often death (Brodie 1968; How et al. 2003; Isbister & Kiernan 2005). Despite the fact that TTX is one of the most powerful neurotoxins known (Medinsky & Klaassen 1996), three species of Thamnophis prey on sympatric newts (Brodie & Brodie 1990; Brodie et al. 2005) and have independently evolved high tolerance of TTX (Feldman et al. 2009).

The physiological and genetic mechanisms at least partially responsible for elevated TTX resistance involve slight alterations in the outer pore (P-loop) of the skeletal muscle sodium channel (Nav1.4) that dramatically reduce the affinity of TTX to this protein (Geffeney et al. 2002, 2005; Feldman et al. 2009). Certain P-loop replacements in Nav genes alter TTX ligation to sodium channels and these changes are thought to be at the root of physiological adaptation to TTX (Geffeney et al. 2005; Venkatesh et al. 2005; Jost et al. 2008; Feldman et al. 2009). Thus, Nav loci probably constitute genes of major effect, yet we still lack a firm grasp on the contribution of individual alleles to TTX resistance. Although all populations of resistant garter snakes display some phenotypic variation, only a few exhibit the high degree of variance suggestive of multiple underlying NaV genotypes. Here, we take advantage of striking variability in TTX resistance seen in T. atratus and T. sirtalis along the central coast of California. We focus on these two species at a local scale so that intraspecific comparisons occur between individuals with similar genetic backgrounds. We examine the relationship between allelic variation in Nav1.4 and TTX resistance in these populations, and how this relationship bears on our understanding of adaptive evolution in garter snakes and the coevolutionary dynamics between newts and snakes.

2. Material and methods

(a) Bioassays

We collected TTX resistance data from 76 T. atratus from Santa Lucia Preserve, Monterey Co. and 47 T. sirtalis from Gilroy, Santa Clara Co. and Santa Lucia Preserve, Monterey Co.

We measured TTX resistance using a highly repeatable and previously validated bioassay of whole organism performance (Brodie & Brodie 1990; Brodie et al. 2002b). We first established an individual's ‘baseline speed’ by racing it down a 4 m racetrack equipped with infrared sensors. We averaged the speed of two time trials to obtain an individual's baseline crawl speed. Following a day of rest (48 h between injections), we gave each snake an intra-peritoneal injection of a known, mass-adjusted dose of TTX (Sigma). Thirty minutes later we raced snakes again to determine ‘post-injection speed.’ We repeated this process, again resting snakes for a day (48 h between injections) then increasing the dose of TTX (0.5, 1, 2, 5 and 10 µg) and running snakes, up to five total sequential TTX tests per snake. We scored ‘resistance’ as the reduction of an individual's baseline sprint speed following an injection of TTX (post-injection speed/baseline speed). We calculated an individual's dose–response curve from the serial TTX injections using a simple linear regression (Ridenhour et al. 2004). From this regression model we estimated the ‘50 per cent dose,’ the amount of TTX required to reduce the snake to 50 per cent of its baseline speed. Because TTX resistance is related to body size (Brodie & Brodie 1990; Ridenhour et al. 2004), we transformed doses into mass-adjusted mouse units (MAMU), the amount of TTX (mg) required to kill a 20 g mouse in 10 min (see Ridenhour et al. 2004).

Some of the T. sirtalis were so resistant to TTX that they were essentially unaffected by our standard injections. We administered additional TTX doses to a subset of these snakes and found they could still run above 50 per cent of their normal ability at over 200, 500 and 1000 MAMUs. We thus assigned a minimum 50 per cent estimate of 100 MAMUs to all unaffected T. sirtalis (based on the highest common dose given), recognizing that the actual measures of TTX resistance must be much higher. Because we assigned a threshold value of resistance, we could not estimate variance for samples including these individuals (i.e. T. sirtalis genotypes).

(b) Sequence data

We collected DNA sequence data from Nav1.4 from each garter snake assayed for TTX resistance (T. atratus n = 76; T. sirtalis n = 47). The entire a-subunit of Nav1.4 in garter snakes encodes 1875 residues (5658 bp) and shows high structural and amino acid homology, and conservation of intron/exon boundaries with mammalian Nav1.4 (Feldman et al. 2009). However, we focused on variation in portions of the four domains (DI-DIV) that code for the outer pore (P-loops) because TTX interacts with residues of the outer pore (Lipkind & Fozzard 2000; Hille 2001) and changes at some of these sites in Nav1.4 are thought to contribute to TTX resistance in snakes (Geffeney et al. 2005; Feldman et al. 2009) and pufferfish (Venkatesh et al. 2005; Jost et al. 2008; Maruta et al. 2008).

We isolated and purified genomic DNA from muscle or liver tissue with the DNeasy Tissue Kit (Qiagen, Inc.). We amplified the four P-loops of Nav1.4 (figure 1) using primers we designed specifically to capture the regions between the S5 and S6 transmembrane segments (Feldman et al. 2009) that form the outer pore, yielding five amplicons of portions of the four domains (two from DIII) as follows: DI approximately 900bp (partial exon 8, all intron 8, partial exon 9); DII approximately 210bp (partial exon 13); DIIIa approximately 450bp (partial exon 19, partial intron 19) DIIIb approximately 800–1600bp (partial intron 19, all exon 20, all intron 20, partial exon 21); DIV approximately 450bp (partial exon 24). We cleaned amplified products using the ExcelaPure PCR Purification Kit (Edge Biosystems) and used purified template in cycle-sequencing reactions with Big Dye 3.1 (Applied Biosystems, Inc.). Following an isopropanol/ethanol precipitation, we ran cycle-sequenced products on an ABI 3130 or ABI 3730 DNA Analyzer (Applied Biosystems, Inc.). We sequenced all samples in both directions.

Figure 1.

Structure of the a-subunit of the skeletal muscle sodium channel (Nav1.4) and the functional variation that characterizes Nav1.4 alleles in the garter snakes Thamnophis atratus and T. sirtalis. (a) The membrane bound Nav1.4 protein is composed of four domains (DI–DIV) with six transmembrane segments (S1–S6) and linkers that connect segments. The four, polypeptide chains linking S5 to S6 (bold) are referred to as the pore loops (P-loops) because they form the cone-shaped outer pore of the channel, allowing the selective permeation of Na+ ions vital for action potential generation and propagation. However, a number of P-loop residues bind strongly to the neurotoxin TTX, which occludes the pore and halts Na+ movement. Derived P-loop replacements (grey circles) occur at critical residues in DIII and DIV that change the structure and electrostatic environment of the pore and alter TTX binding affinity (see text). (b) Nav1.4 alleles identified in central California populations of T. atratus and T. sirtalis alongside the DIII and DIV P-loop amino acid sequences that characterize the alleles, and the frequencies at which these alleles were found. A TTX-sensitive T. elegans displaying the ancestral garter snake Nav1.4 sequence (1.4+) is provided as a reference, and a human sequence is given for comparison (GenBank M81 758). Structures of the pore (*, selectivity filter; α, α-helix; β, β-strand) follow Lipkind & Fozzard (2000); invariant DI and DII P-loop sequences not shown.

We edited and aligned sequences in Sequencher 4.2 (Gene Codes Corp.) and translated coding regions into amino acid sequences using MacClade 4.08 (Maddison & Maddison 2005). We deposited all sequences in GenBank (HM 365341–HM 365925).

We scored a Nav1.4 sequence as a unique allele if it possessed amino acid replacements in or near the pore forming structures of the P-loops (pore a-helix, selectivity filter, β-strand; Lipkind & Fozzard 2000) that interact with TTX. Our notation for alleles includes the numerical designation of the sodium channel family member (i.e. 1.4) followed by a superscript of one letter amino acid abbreviations given in the order those derived allelic substitutions occur in the locus. This nomenclature reflects the molecular differences between the ancestral garter snake gene sequence (here termed 1.4+) and derived Nav1.4 variants, rather than putative phenotypic effects or dominance attributes of alleles.

(c) Microsatellite data

To provide a rough estimate of background levels of variation in the garter snake genome to contrast against patterns in our candidate locus, we collected fragment data from putatively neutral loci for half our T. atratus sample (n = 38) and all our T. sirtalis sample (n = 47). We amplified 16 microsatellite markers developed for natricines (Nerodia and Thamnophis): NSU2, NSU3, NSU7, NSU8, NSU10 (Prosser et al. 1999), Ts2, Ts3, Ts4 (McCracken et al. 1999), 2Ts (Garner et al. 2002), Te1Ca2, Te1Ca18, Te1Ca29, Ts1Ca4 (Garner et al. 2004), TEO51B, TSO10, TSO42 (Manier & Arnold 2005). Microsatellites were end-labelled (forward primer) with one of four fluorescent tags (HEX, NED, 6FAM, PET), allowing us to multiplex amplified products (2–5 markers per panel) and resolve genotypes on an ABI 3730 DNA Analyzer (using LIZ size standard). We re-genotyped samples with rare alleles, faint signal, or ambiguous peaks, and checked for error by re-genotyping a subset of all samples. We scored alleles with GeneMarker v1.85 (SoftGenetics).

(d) Genotypic analyses

We calculated numbers of alleles (Na), frequencies, and observed (Ho) and expected (He) heterozygosities for each locus in GenAlEx v6.2 (Peakall & Smouse 2006). We tested for deviations from Hardy–Weinberg equilibrium (HWE) in each locus using the exact test in GenePop v4 (Raymond & Rousset 1995; Rousset 2008). We also assessed linkage disequilibrium (LD) among loci using the likelihood ratio test in GenePop. We estimated p-values for HWE and LD tests in GenePop with the Markov chain method, runing100 batches with 1000 iterations per batch, then adjusted p-values using the sequential Bonferroni correction (Rice 1989). Finally, we assessed the degree of population structure within our T. atratus and T. sirtalis data by calculating the inbreeding coefficient (FIS) in GenAlEx.

(e) Phenotype–genotype matching

We assigned each phenotype to its respective genotype for each locus and then determined the statistical association between genotypes (allelic variation in Nav1.4 or in microsatellites) and phenotypes (variation in TTX resistance) with a single factor analysis of variance (ANOVA) using restricted maximum likelihood in SAS v9.2 (SAS Institute). This method also yields an estimate of the percentage of variance explained (PVE) by the genotype, that is, the proportion of the total trait variance attributed to the genotype (e.g. Hoekstra et al. 2006). We then tested for overall differences in TTX resistance between Nav1.4 genotypes with post hoc tests to assess the significance of pairwise differences between Nav1.4 genotypes; we used t-tests of the least squares means from the single factor ANOVA to assess differences among T. atratus genotypes, and Chi-square tests of the Wilcoxon scores from a Kruskal–Wallis one-way ANOVA to assess differences among T. sirtalis genotypes. We used the non-parametric Kruskal–Wallis test for T. sirtalis to accommodate the truncated values of resistant snakes.

We qualitatively note apparent patterns of allelic dominance in Nav1.4 but did not formally evaluate dominance here because we could not estimate phenotypic variances around some T. sirtalis genotypes (see above) and because we did not identify any fully heterozygous T. atratus (see results).

3. Results and discussion

(a) Variation in TTX resistance

Thamnophis atratus and T. sirtalis from central coastal California possess extensive variation in TTX resistance. For example, the amount of TTX required to slow T. atratus to 50 per cent of their normal crawl speed ranged from 2 to 100 MAMUs, likewise, 50 per cent TTX dosages in T. sirtalis ranged from 6 to more than 100 MAMUs. In fact, the range of phenotypes observed here is equivalent to the range of TTX resistance seen across the entire Thamnophis clade, with some T. atratus and T. sirtalis possessing ancestral, low levels of TTX resistance (e.g. 2–5 MAMUs; Motychak et al. 1999; Feldman et al. 2009) and others exhibiting some of the most dramatic resistance ever recorded (e.g. more than 100 MAMUs; Brodie et al. 2002a; Feldman et al. 2009).

(b) Variation in Nav1.4 and neutral markers

We found four Nav1.4 alleles in our central coastal population of T. atratus (1.4+, 1.4EPN, 1.4EP, 1.4N; figure 1). The first, 1.4+, is the same TTX-sensitive Nav1.4 allele (‘wild-type’) found across the garter snake phylogeny (Feldman et al. 2009). The three derived T. atratus alleles can be distinguished from 1.4+ as follows: 1.4EPN contains a D1277E in the DIII β-strand, an A1281P C-terminal to the DIII β-strand, and a D1568N in the DIV β-strand; 1.4EP holds the two DIII changes, D1277E and A1281P; 1.4N possesses only the DIV replacement, D1568N. The four alleles are distributed among only four genotypes in this T. atratus population (figure 2).

Figure 2.

TTX resistance in T. atratus (dark grey) and T. sirtalis (light grey) allocated by Nav1.4 genotype. For each genotype, mean TTX resistance, standard deviation (error bars), and observed versus expected genotypes are given. Note that some T. atratus genotypes were not detected in our sample (expected proportions do not total 100%). Due to the heightened resistance of T. sirtalis with derived Nav1.4 alleles, we were forced to assign a minimum threshold value of TTX resistance and could not estimate phenotypic variance for these individuals (see text).

Though Nav1.4 in this population of T. atratus displays similar allelic richness to the neutral markers, actual heterozygosity in Nav1.4 is lower than all but one variable microsatellite (table 2). Furthermore, Nav1.4 is severely out of HWE, consistent with the premise that this gene is under strong natural selection. Only one microsatellite is out of HWE despite the fact that some non-random mating is evident (FIS = 0.186), but this marker is not associated with Nav1.4 (all markers statistically unlinked; data not shown) or the phenotype (table 2).

We found only two Nav1.4 alleles in central coastal T. sirtalis (1.4+, 1.4LVNV; figure 1). As with T. atratus, some T. sirtalis possess the same functional allele that appears to represent the ancestral garter snake Nav1.4 sequence (Feldman et al. 2009). The derived T. sirtalis allele differs from this putative ancestor through four mutations in the P-loop of DIV; 1.4LVNV is characterized by an I1555L and I1561V in the a-helix and a D1568N and G1569V in the β-strand. We found the two T. sirtalis Nav1.4 alleles in all three possible genotypic combinations (figure 2).

As in T. atratus, Nav1.4 in sympatric T. sirtalis also exhibits a similar number of alleles as neutral loci, but generally displays lower overall heterozygosity (table 2). The locus is also out of HWE, however, so are four microsatellites (table 2). It is not clear if such widespread disequilibrium is the hallmark of selection in this genome (all markers statistically unlinked; data not shown), or simply the fact that our sample of T. sirtalis does not represent a single panmictic unit, though FIS is negligible (0.086), suggesting little structure attributable to non-random breeding.

(c) Relationship between variation in Nav1.4 and TTX resistance

We identified five Nav1.4 alleles (four derived) in these garter snakes. Three of the Nav1.4 alleles have been identified and detailed previously (Geffeney et al. 2005; Feldman et al. 2009), and the two T. sirtalis alleles have been functionally expressed (Geffeney et al. 2005). Most of the P-loop replacements that characterize the derived Nav1.4 alleles have been documented to reduce the binding affinity of TTX to sodium channels at the molecular level (Terlau et al. 1991; Pérez-García et al. 1996; Chen et al. 1997; Penzotti et al. 1998; Choudhary et al. 2003; Geffeney et al. 2005; Jost et al. 2008; Maruta et al. 2008), and a number are coincident with those seen in TTX-bearing pufferfish (Venkatesh et al. 2005; Jost et al. 2008; Maruta et al. 2008). Individual P-loop substitutions modulate TTX ligation by changing the structure and molecular environment of the channel's outer pore, reducing the number and strength of chemical bonds as well as dampening steric attraction between TTX and the outer pore (e.g. Tikhonov & Zhorov 2005; Scheib et al. 2006). For example, the D1568N replacement observed in both snake species occurs at a site known to play a major role in TTX ligation (Terlau et al. 1991; Pérez-García et al. 1996; Chen et al. 1997; Penzotti et al. 1998; Choudhary et al. 2003). This specific amino acid substitution has been shown to increase TTX resistance 30–40-fold in rat Nav1.4 (Penzotti et al. 1998; Choudhary et al. 2003), probably because the hydrogen bond normally formed between TTX and D1568 (Chen et al. 1997; Choudhary et al. 2003; Scheib et al. 2006) is neutralized by uncharged N, and because the larger functional group of N alters the volume of the outer pore. Nevertheless, a direct relationship between allelic variation in Nav1.4 and TTX resistance at the organismal level has yet to be established.

Partitioning T. atratus phenotypes into their respective Nav1.4 genotypes reveals obvious phenotypic differences between the groups (table 1; figure 2). Snakes homozygous for the ancestral Nav1.4 allele exhibit the lowest TTX resistance (=7.41, s.d. = 3.99) while the next lowest group is composed of snakes heterozygous for 1.4N and 1.4+ (=16.47, s.d. 18.08). These two genotypes display considerably lower TTX resistance (t-value = −5.06, p < 0.0001; t-value = −3.77, p = 0.0003) than snakes with the 1.4EPN allele (=47.84, SD 39.80). Whether the effects of the individual amino acid substitutions is additive remains to be determined. Nevertheless, the association between Nav1.4 polymorphism and TTX resistance is significant in T. atratus (F3,78 = 6.764, p < 0.0004), and 23.4 per cent of the phenotype is explained by the Nav1.4 genotype (table 2). This pattern of association contrasts that seen among neutral loci, which show no relationship with TTX resistance (table 2).

View this table:
Table 1.

Differences in TTX resistance between Nav1.4 genotypes in central California populations of T. atratus and T. sirtalis. Statistically significant differences between genotypes in italic.

Matching T. sirtalis phenotypes to their respective Nav1.4 genotypes also reveals dramatic differences among the genotypes (table 1; figure 2). Snakes homozygous for 1.4+ show moderate levels of TTX resistance (=18.97, s.d. = 8.58) while those heterozygous for 1.4LVNV and 1.4+ or homozygous for 1.4LVNV display exceedingly high resistance (>100.00; >100.00). Snakes with these latter two genotypes are significantly less affected by TTX than those homozygous for 1.4+ (χ2-value = 21.98, p < 0.0001; χ2-value = 29.53, p < 0.0001), but indistinguishable from one another, at least as quantified here (χ2-value = 0, p = 1). Thus, the derived 1.4LVNV allele appears to exhibit complete dominance. However, our capacity to detect dominance is limited by our inability to estimate the phenotype variance around these two genotypes (see §2). At the molecular genetic level, it seems reasonable to expect that heterozygotes would be less resistant to TTX because presumably half of the skeletal muscle sodium channels in these snakes should be translated from the 1.4+ allele. Yet at the organismal level, T. sirtalis heterozygous or homozygous for 1.4LVNV appear so immune to TTX that these genotypes may be functionally equivalent. The relationship between variation in Nav1.4 and TTX resistance in these T. sirtalis populations is strong (F2,44 = 1330.820, p < 0.0001; we cannot consider the PVE of 98% from our ANOVA reliable given our truncated dataset). Conversely, our microsatellite markers, which provide a rough measure of background variation in the garter snake genome, show no association with the phenotype (table 2), supporting the conclusion that the robust link between Nav1.4 genotype and phenotype is not an artefact of hidden population structure.

View this table:
Table 2.

Genetic variation in Nav1.4 and 16 microsatellite markers in central California T. atratus and T. sirtalis as measured in number of alleles (Na) and observed heterozygosities (Ho), along with tests for Hardy–Weinberg equilibrium (HWE) and relationships between genotype and phenotype (using a single factor ANOVA). Loci that did not produce reliable amplicons denoted as not applicable (N/A); loci that were monomorphic denoted with dash.

(d) Adaptation and gene family evolution

We show that Nav1.4 is clearly a large effect locus, explaining 23 per cent of TTX resistance in T. atratus and correlating completely with qualitative differences in TTX resistance in T. sirtalis. Nevertheless, allelic variation in Nav1.4 cannot solely account for whole animal resistance in garter snakes. First, the proportion of the phenotype explained by the genotype in T. atratus suggests the contribution of additional loci or molecular mechanisms. Second, the levels of TTX resistance displayed by T. sirtalis that are homozygous for the ancestral Nav1.4 allele (1.4+) are higher than T. atratus and other garter snakes with this same non-resistant genotype (Feldman et al. 2009). Third, Nav1.4 expression is confined to skeletal muscle tissue (Trimmer et al. 1989, 1990; Goldin 2001), yet the central nervous system and some peripheral nerves are natively sensitive to TTX (Goldin 2002). It seems evident that these tissues must also be defended against TTX to produce a fully resistant organism.

In mammals, nine different Nav genes are functionally expressed in specific tissues with excitable cells (Plummer & Meisler 1999; Goldin 2001; Hille 2001). A third of these (Nav1.5, Nav1.8, Nav1.9) are natively insensitive to TTX (Goldin 2002), while the remaining sodium channels (Nav1.1, Nav1.2, Nav1.3, Nav1.4, Nav1.6, Nav1.7) are highly sensitive to the toxin (Goldin 2002). Modification of other TTX-sensitive Nav loci may have encountered adaptive mutations similar to those in the derived Nav1.4 alleles. Because even single mutations can dramatically alter TTX binding affinity (e.g. Noda et al. 1989; Terlau et al. 1991; Geffeney et al. 2005; Tikhonov & Zhorov 2005; Venkatesh et al. 2005; Du et al. 2009), only slight modification of central and peripheral nerve Nav genes may be required. Indeed, parallel changes have occurred across Nav paralogues in TTX-bearing pufferfish (Jost et al. 2008). Interestingly, in vertebrates most Nav paralogues are physically close to their nearest relatives (Plummer & Meisler 1999; Goldin 2002; Novak et al. 2006), and in mammals, seven of the nine genes are regionally clustered within two chromosomes (Plummer & Meisler 1999; Goldin 2002). Thus, coincident P-loop changes in separate sodium channel genes might immediately become linked. If recombination rates are low, selection could easily fix such linkage groups in populations before recombination breaks apart adaptive allelic combinations.

We suspect that convergent gene family evolution (e.g. Jost et al. 2008; Storz et al. 2009) explains whole animal TTX resistance in garter snakes and suggest that both T. atratus and T. sirtalis demonstrate the extremes of this scenario. In the case of T. sirtalis, derived Nav variants may have become (or remain) linked, forming a coadapted gene complex that produces an extreme phenotype. In T. atratus, on the other hand, recombination may have disassociated the linkages between TTX resistant alleles at separate Nav loci. Some T. atratus would still inherit this suite of adaptive variation and possess the ‘Nav supergene’ that produces acute resistance. Others would possess varying degrees of mismatch between resistant Nav genes and express corresponding levels of TTX resistance. Alternatively, species differences in resistance variation might exist because of differences in the age of the ecological interaction with newts or of the mutations themselves. Identification and characterization of the P-loops in all TTX-sensitive sodium channels in Thamnophis should eventually provide a complete picture of the molecular mechanisms responsible for elevated TTX resistance.

(e) Simple genetic architecture and coevolution

The molecular basis of adaptive TTX resistance in Thamnophis appears relatively simple, probably involving a few replacements in a handful of genes. Simple genetic architecture underlying the traits that mediate coevolution should allow for a rapid back-and-forth evolutionary arms race, but also potential ‘escape’ from the arms race by the predator (Sasaki 2000; Kopp & Gavrilets 2006; Nuismer et al. 2007). Though we only have a handle on the genetics underlying the predatory counter adaptation, both of these predictions appear to have been met. Although coevolution between toxic newts and predatory garter snakes has deep roots, and probably extends beyond the diversification of many contemporary lineages (Feldman et al. 2009), a number of T. sirtalis and T. granulosa populations in the Pacific Northwest seem to have initiated an arms race only recently, probably during the Late Pleistocene (Ridenhour et al. 2007). A range-wide survey of trait matching between T. sirtalis and Taricha showed that one-third of western T. sirtalis populations possess such acute TTX resistance that they are essentially ‘winning’ the arms race (Hanifin et al. 2008). Characterization of additional Nav loci in Thamnophis and continued phenotype–genotype matching at the population level will help determine whether the evolution of a Nav supergene in garter snakes frequently creates an insurmountable counter measure to which newts cannot respond, or whether recombination and gene flow consistently decouple adaptive alleles and keep newts competitive in this arms race.


All research was conducted under approval of Utah State University IACUC (Protocol no. 1008). Snakes injected with small doses of TTX fully recover after several hours (Brodie et al. 2002b), and as a sodium channel blocker, TTX impedes pain reception (Lai et al. 2003). At the completion of the study, all specimens were humanely euthanized and deposited as voucher specimens in the California Academy of Sciences (CAS) and University of Texas at Arlington (UTA) herpetology collections.

We thank M. Edgehouse, A. Mortensen and A. Wilkinson for assistance in the field and R. Lawson for field advice. We are grateful for access to land by the Santa Lucia Preserve. For aid with captive specimens and bioassays we thank A. Mortensen, J. Scoville, A. Wilkinson, J. Pluid, M. Edgehouse, and approval of protocols from USU IACUC. We thank J. Vindum and M. Koo (CAS), and J. Campbell and C. Franklin (UTA) for help with the curation of specimens. For assistance in the laboratory we thank M. Matocq, A. Runck, and M. Andrews and E. O'Leary-Jepsen (ISU MRCF), and K. Kruse and C. Osborne (UNR NGC). We appreciate statistical aid from L. Latta, K. Stewart and P. Murphy, and useful discussions with M. Matocq, E. O'Neill, and the UVa and USU Herp Groups, as well as helpful comments from J. Storz and three anonymous reviewers. For providing scientific collecting permits we thank California Department of Fish and Game. This work was supported by NSF grants to MEP, EDB Jr and EDB III (DEB-0315172, DEB-0922251), a NGS grant (7531-03) to EDB Jr and EDB III, and a USU School of Graduate Studies Dissertation Fellowship to CRF.

  • Received April 8, 2010.
  • Accepted May 14, 2010.


View Abstract