Genetic Analysis and Captive Breeding Program Design for the Eastern Massasauga Sistrurus catenatus catenatus

Conservation genetic analyses have clear management implications in the design and implementation of captive breeding and reintroduction programs (Kozfkay et al. 2008; Sharma et al. 2009; Tzika et al. 2009; Witzenberger and Hochkirch 2011). The eastern massasauga Sistrurus catenatus catenatus (Figure 1), a subspecies of the massasauga rattlesnake, ranges across the North American Midwest and is threatened by habitat destruction, fragmentation, and persecution (Szymanski 1998; Rouse and Willson 2002). Except in Michigan, where it is a species of special concern, the eastern massasauga is listed as threatened or endangered within each state or province in which it occurs (Illinois, Indiana, Iowa, Michigan, Minnesota, New York, Ohio, Pennsylvania, Wisconsin, Ontario; Levell 1998; Rouse and Willson 2002) and is a candidate for federal listing under the U.S. Endangered Species Act (ESA 1973, as amended).

In 2006, the Association of Zoos and Aquariums (AZA) established a Species Survival Plan (SSP®) aimed at increasing the size and preserving genetic diversity of the captive population as a possible source in future reintroduction and augmentation (Earnhardt et al. 2011). Key elements of the SSP include development of a studbook detailing the pedigree of the managed population and development of a science-based breeding and transfer plan. Achieving these objectives has been hampered by the unknown origin of some SSP animals and incomplete information on the population genetic structure of natural populations.

Previous genetic analyses of the eastern massasauga using microsatellite DNA loci demonstrated that population subdivision occurs over short distances (<25 km; Gibbs et al. 1997; Andre 2003; Kropiewnicki 2008; Chiucchi and Gibbs 2010). In addition, analysis of nuclear and mitochondrial DNA sequences provided evidence that the eastern massasauga warrants elevation to full species status as distinct from the western massasauga S. c. tergeminus and desert massasauga S. c. edwardsii (Kubatko and Gibbs 2010; Kubatko et al. 2011). Incidentally, because species are ranked higher than subspecies and populations, such a taxonomic change would elevate the eastern massasauga priority rank under the ESA (ESA 1973). Genetic analyses have also provided clarification of the taxonomic identity of extant populations in northwestern Missouri as western massasaugas (Gibbs et al. 2010; Gerard et al. 2011; Kubatko et al. 2011). Previously, these snakes had been variably classified as S. c. catenatus, S. c. tergeminus, or hybrids between the two (Evans and Gloyd 1948; Szymanski 1998). However, a range-wide analysis of the eastern massasauga, including representatives of most extant populations and coupled with parallel analysis of captive snakes, has not been conducted. We utilized mitochondrial DNA sequence data to identify geographic structure among wild populations, clarify geographic boundaries among taxonomic units, and characterize the genetic make-up of captive individuals held at AZA institutions as part of the captive breeding program. We synthesize our results with those of other studies and discuss their implications for captive breeding of the eastern massasauga.

Tissue and DNA samples from living animals across the range of the eastern massasauga were provided by colleagues and SSP participating institutions (Figure 2; Table S1, Supplemental Material; Ray 2009). We analyzed 186 snakes: 109 wild snakes from 34 U.S. and Canadian counties and districts, all 52 members of the SSP breeding program (23 of known and 29 of unknown origin), 18 other captives of unknown origin that have since died or are otherwise excluded from breeding, and 7 outgroup representatives. Outgroup samples include western and desert massasaugas (n  =  7) and a single sample from Pottawattamie County, Iowa, a disjunct population whose taxonomic status was unclear at the start of this study (Table S1, Supplemental Material). Genomic DNA was extracted using DNeasy Blood and Tissue Kits and Puregene Core Kits (Qiagen, Inc., Valencia, CA) following manufacturer protocols and used as template for polymerase chain reaction (PRC) amplification of the mitochondrial gene NADH dehydrogenase subunit II (ND2). Reactions were performed using Failsafe PCR systems premix H (Epicentre Technologies Corp., Madison, WI) with primers CE2330 (5′-CTA ATA AAG CTT TCG GGC CCA TAC-3′) and CE2331 (5′-TTC TAC TTA AGG CTT TGA AGG C-3′) from Janzen et al. (2002).

We began amplification reactions with a 2-min denaturing step at 94°C, followed by 40 cycles of amplification using a 30-s 94°C denaturation step, a 30-s 50°C annealing step, and a 90-s 72°C elongation step. After final elongation for 7 min at 72°C, we purified amplicons using ExoSAP-IT (USB Corp., Cleveland, OH) and sequenced them using a commercial service (Macrogen, Inc., Seoul, South Korea). Sequencing reactions used amplification primer CE2330 and an internal primer, L5238s (5′-ACT TGA CAG AAA ATT GCC CCC-3′) modified from the L5238 primer (de Queiroz et al. 2002).

We inspected, trimmed, and aligned sequences using Geneious v5 (Drummond et al. 2009). We used the computer program TCS 1.21 (Clement et al. 2000) to create an unrooted statistical parsimony network connected at a minimum 95% significance level. We generated gene trees using maximum likelihood and maximum parsimony. The program jModelTest (Posada 2008) identified the GTR+G evolutionary model as most appropriate for the data set. We conducted maximum likelihood analyses in MEGA4 (Tamura et al. 2007) with 1,000 bootstrap replicates. We conducted maximum parsimony analyses in PAUP* 4.0b10 (Swofford 2002) with 1,000 bootstrap replicates and 100 random addition sequence replicates per bootstrap replicate. A 3–base-pair (bp) deletion observed in some haplotypes (see Results) was treated as both a 1-bp and a 3-bp substitution in maximum likelihood analyses, and as a 1-bp and a 3-bp gap using fifth-state coding in maximum parsimony analyses (Giribet and Wheeler 1999; Ogden and Rosenberg 2007; Simmons et al. 2007). We calculated nucleotide and haplotype diversity, tests for selective neutrality (Tajima's D, Fu and Li's F, Fu and Li's D), and tests for population expansion (Fu's F, Ramos-Onsins and Rozas' R²) in DnaSP 5 (Librado and Rozas 2009). We calculated net sequence divergence and transition/transversion ratios using MEGA5 (Tamura et al. 2007). We quantified the degree of exclusive ancestry of geographic ND2 subunits using the genealogical sorting index (gsi; Cummings et al. 2008) with 10,000 permutations. This index quantifies the historical relationship among groups without the requirement of reciprocal monophyly; here, this is an advantage given the recency of the hypothesized range expansion of the eastern massasauga into formerly glaciated areas. We further explored phylogeographic structure by comparing N_ST, an ordered measure of genetic differentiation that incorporates distances among haplotypes, and G_ST, an unordered measure of genetic differentiation, using Permut (Pons and Petit 1996) with 1,000 permutations. When haplotypes within a geographic ND2 subunit are more similar in sequence than expected given their overall frequency, N_ST exceeds G_ST, thus providing a signal of geographic subdivision (Pons and Petit 1996).

We analyzed 922 bp of ND2 from 186 massasaugas and identified 23 unique haplotypes exhibiting 133 synonymously variable bases (Table S1, Supplemental Material, GenBank Accession No. GQ359794–GQ359809, GQ359813–GQ359814, JQ609659–JQ609663). Five haplotypes contained the same 3-bp deletion, located within the reading frame. Whereas the other haplotypes contained a tandem repeat (CCTCCT), these five haplotypes contained only a single copy (CCT), suggestive of a deletion due to slipped-strand mis-pairing (Graur and Li 2000). The program MEGA5 estimated the transition/transversion bias to be R  =  4.43. Unrooted statistical parsimony network, maximum likelihood, and maximum parsimony analyses all identified the eastern massasauga haplotypes as divergent from western and desert massasauga haplotypes (Figure 3). Maximum parsimony analyses yielded the same most parsimonious trees regardless of whether the deletion was treated as a 1-bp or a 3-bp gap (tree lengths  =  144 and 146 steps, respectively). The maximum parsimony majority rule consensus tree was topologically congruent to the trees obtained from maximum likelihood analyses, regardless of whether the deletion was treated as a 1-bp or 3-bp substitution (Figure 4). Net sequence divergence equaled 14% between eastern massasaugas and western+desert massasaugas. Western massasaugas and desert massasaugas differed by just 1%. Haplotypes of a wild-caught snake from Pottawattomie County, Iowa, and four captive snakes lacking locality data (two law-enforcement confiscations and two from a private collection) were similar to western and desert massasauga haplotypes. As a result, the four captive snakes have been excluded from the eastern massasauga SSP.

We identified 18 haplotypes among the 179 eastern massasaugas included in our analysis. An unrooted statistical parsimony network revealed three ND2 subunits (Figure 3), corresponding to separate geographic regions (Figure 2): (1) a western ND2 subunit in Iowa, Wisconsin, and Illinois (n  =  42 animals of known origin, number of haplotypes H  =  5, haplotype diversity h  =  0.656, nucleotide diversity π  =  0.0010); (2) a central ND2 subunit in Indiana, southern and central Michigan, Ohio, and far southwestern Ontario (n  =  56 and 34 animals of known and unknown origin, respectively; H  =  8, h  =  0.357, π  =  0.0005); and (3) an eastern ND2 subunit in New York, Pennsylvania, northern Michigan, and remaining portions of Ontario (n  =  34 and 13 animals of known and unknown origin, respectively; H  =  5, h  =  0.412, π  =  0.0005). Net sequence divergence was low, equaling 0.5% between western and central ND2 subunits, 0.8% between western and eastern ND2 subunits, and 0.4% between central and eastern ND2 subunits. Western and central ND2 subunits were separated from each other by a minimum of three mutational steps and from the eastern ND2 subunit by the 3-bp deletion. Maximum likelihood and maximum parsimony gene trees provided moderate support for the monophyly of the central+eastern ND2 subunits (bootstrap value  =  75–77, depending on how the deletion was scored) and variable support for the monophyly of the eastern ND2 subunit (bootstrap value  =  61–96, depending on how the deletion was scored; Figure 4).

Divergence between eastern and western+desert massasaugas and division of the eastern massasauga into western, central, and eastern ND2 subunits was supported by gsi values that were all equal to 1.0 and significantly different from 0.0 (Figure 4). Thus, despite the presence of polytomys and sometimes only moderate bootstrap support using maximum likelihood and maximum parsimony, evidence for exclusive ancestry of eastern massasaugas from western+desert massasaugas, of the central and eastern ND2 subunit from the western ND2 subunit, and of the eastern ND2 subunit from the central ND2 subunit, is strong (Figure 4). Similarly, within the eastern massasauga, ordered genetic differentiation among ND2 subunits was significantly greater than unordered genetic differentiation (N_ST  =  0.82, G_ST  =  0.56, P < 0.001), indicating that individuals within ND2 subunits share haplotypes that are more similar in sequence than would be expected in the absence of geographic structure.

Demographic analysis indicates population stability in the western ND2 subunit (Fu's F  =  −0.63, P  =  0.185; Ramos-Onsins and Rozas' R²  =  0.12, P  =  0.512) but population expansion in the central (F  =  −6.14, P  =  0.002; R²  =  0.03, P  =  0.002) and possibly eastern ND2 subunits (F  =  −2.44, P  =  0.059; R²  =  0.06, P  =  0.056). Tests for selection (Tajima's D, Fu and Li's F, and Fu and Li's D) were all nonsignificant.

Of 52 snakes included in the eastern massasauga SSP breeding plan as of January 2012, 29 were of unknown origin. Our analysis allowed us to assign these animals unambiguously to one of the three mtDNA ND2 subunits identified above. These 52 snakes included 12, 26, and 14 western, central, and eastern ND2 subunit representatives, respectively. Three of five western ND2 subunit haplotypes, three of eight central ND2 subunit haplotypes, and two of five eastern ND2 subunit haplotypes were present among them.

Our results support the conclusions of Kubatko et al. (2011) regarding the taxonomic distinctness of the eastern massasauga. Although we analyzed a single mitochondrial gene (vs. 18 nuclear genes and 1 mitochondrial gene in Kubatko et al. 2011), our sampling of the eastern massasauga clade was more extensive, including 132 snakes from 34 U.S. and Canadian counties and districts (vs. 9 eastern massasaugas in Kubatko et al. 2011). Furthermore, our inclusion of samples from Bremer County in northeastern Iowa (n  =  5; belonging to the eastern massasauga clade) and Pottawattamie County in southwestern Iowa (n  =  1; belonging to the western+desert massasauga clade) further clarifies the geographic distribution of these taxa. Placement of the Pottawattamie County sample in the western+desert massasauga clade is supported by independent analysis of two other mitochondrial genes (H. L. Gibbs, The Ohio State University, personal communication). As with populations in western Missouri (Gibbs et al. 2010; Gerard et al. 2011; Kubatko et al. 2011; Wooten and Gibbs 2011), our results indicate that Pottawattamie County, Iowa, should be excluded from eastern massasauga management plans.

Our study demonstrates the existence of three geographic ND2 subunits within the eastern massasauga that were not evident in previous analyses that were largely based on sequence data from nuclear genes (Kubatko et al. 2011). Kubatko et al. (2011) included one mitochondrial marker in their analyses (296 bp of ATP 6 and 8); however, low levels of variation rendered it uninformative with respect to genetic structure within eastern massasaugas (figure S2.a in Kubatko et al. 2011). Larger sample size and more complete geographic coverage enhanced our ability to detect geographic structure. Our sampling encompassed 34 unique locations (compared with 9 locations in Kubatko et al. 2011), and included coverage of northern Michigan (upper and lower peninsula), Indiana, Iowa, and the north side of Lake Erie, which provided data critical in resolving geographic ND2 subunits. In addition, much of the current distribution of the eastern massasauga was ice-covered during the Wisconsinan glacial maximum approximately 18,000 years before present, resulting in a relatively short time-frame over which geographic ND2 subunits could arise. The four-fold lower effective population size of mitochondrial DNA makes this marker more likely than nuclear sequences to reveal geographic structure (Zink and Barrowclough 2008; Barrowclough and Zink 2009).

We find the deletion event that distinguishes our eastern ND2 subunit to be especially noteworthy. This rare genomic change provides phylogenetic information (e.g., Janečka et al. 2007), which, in the present case, is congruent with geography, results from a simple directional mechanism (slipped-strand mis-pairing resulting in a deletion from the ancestral sequence of two tandem repeats), and is unlikely to reflect homoplasy. We observed no other intraspecific deletions among approximately 200 ND2 sequences from four other snake species, nor have intraspecific ND2 deletions been reported in other snakes (Janzen et al. 2002; Placyk et al. 2007), fish (Ptacek and Breden 1998), or birds (Birks and Edwards 2002; Drovetski et al. 2004; Johnson et al. 2005; Moyle and Marks 2006; Pasquet et al. 2006; Rheindt et al. 2008; Jønsson et al. 2010). Based on glacial history and inferred directionality of the deletion event, we interpret the geographic pattern seen in our ND2 sequence data as resulting from range expansion from unglaciated (Iowa, Wisconsin, southern Illinois) into formerly glaciated regions. Demographic analysis suggests that post-glacial range expansion is associated with demographic stability in the western (unglaciated) ND2 subunit and population expansion in the central and possibly eastern ND2 subunits. A lack of samples from northwestern Indiana, where the eastern massasauga is now rare or absent, makes it difficult to determine where the boundary between the western and central ND2 subunits lies. However, boundaries between the central and eastern ND2 subunits fall between collection sites separated by about 32 km in Crawford and Kalkaska counties in northern Michigan, and about 64 km in Ashtabula County, Ohio, and Butler and Venango counties, Pennsylvania.

Genetic analyses of the eastern massasauga from a variety of studies provide different degrees of geographic resolution consistent with the marker type used. Nuclear DNA sequence data from Kubatko et al. (2011) suggest a general lack of geographic structure and only low levels of variation within the eastern massasauga. Mitochondrial sequence data from this study revealed an intermediate level of geographic structure, with three weakly differentiated ND2 subunits. Finally, microsatellite DNA data (Gibbs et al. 1997; Andre 2003; Kropiewnicki 2008; Chiucchi and Gibbs 2010) reveal strong differentiation even at a fine geographic scale, suggesting that local populations are relatively isolated and that gene flow is generally rare.

Although variation in geographic resolution among marker types is expected, these contrasting patterns create a challenge to individuals and agencies tasked with appropriately managing at-risk species. For example, management of captive breeding programs requires balancing demographic and genetic objectives. Captive populations founded by a few individuals and maintained at low numbers rapidly lose genetic variation through stochastic processes and frequently show evidence of inbreeding depression (Allendorf and Luikart 2007). For the eastern massasauga SSP, which includes just 52 animals and is dispersed across 22 zoos, rapid population increase might best be achieved by treating captives as a single unit to immediately maximize the number of breeding pairs; this strategy is supported by the apparent lack of geographic structure based on nuclear DNA sequences. In contrast, maintaining patterns of mitochondrial genetic variation present in wild populations (e.g., for possible reintroduction or augmentation) would require managing the three ND2 subunits separately and may necessitate delayed breeding of some captive individuals. Furthermore, because the number of snakes per ND2 subunit is small, acquisition of additional captives may be desirable. Finally, maintaining patterns of microsatellite DNA variation found in wild populations, along with potential local adaptation, would require managing animals from each source population separately. However, such a strategy would require acquisition of significant numbers of additional snakes (e.g., recent microsatellite DNA analyses identify 29 genetically discrete populations; Kropiewnicki 2008, Chiucchi and Gibbs 2010), possibly putting source populations at greater risk. Furthermore, limited space and resources at participating zoological institutions make it difficult to maintain enough captive individuals to support such an approach.

The benefits of dividing the captive population into genetic subunits must be weighed against practical considerations for numbers of individuals maintained and the need for rapid captive population growth to minimize loss of genetic variation and adaptation to captivity (Robert 2009; Williams and Hoffman 2009). For the eastern massasauga, evidence for long-term population isolation and local adaptation is equivocal. Chiucchi and Gibbs (2010) report low rates of gene flow on both contemporary and historical time scales, making local adaptation more likely. In contrast, Michigan populations (which have lower numbers of unique alleles and fewer unique alleles at high frequency than do populations elsewhere) apparently became isolated more recently, making local adaptation less likely (Kropiewnicki 2008). These practical and genetic considerations led the AZA to adopt a strategy where the three ND2 subunits we identified are being maintained separately during initial breeding efforts (Earnhardt et al. 2011). Importantly, the eastern massasauga SSP is updated annually, allowing changes in breeding design to be incorporated as additional information becomes available. Investigations of functional genes, such as those contributing to immune function (Hughes 1991; Madsen and Ujvari 2006), could be especially useful in identifying the geographic scale of local adaptation and natural subunits for captive breeding. Such genes may also prove useful as monitoring tools should augmentation and reintroduction become part of future eastern massasauga management (e.g., Swanson et al. 2006; Smith and Hughes 2008). Regardless of the breeding design adopted, increasing the number of individuals and diversity of haplotypes among captive individuals through cooperation with organizations outside of the AZA SSP is important. For instance, breeding loans from non-AZA facilities or donations of confiscated animals from law enforcement agencies and other sources could be assigned to ND2 subunits using the tools developed here, thus improving long-term prospects for genetic health and, ultimately, the persistence of the species.

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any supplemental material. Queries should be directed to the corresponding author for the article.

Table S1. Geographic locality in North America, haplotype, and NADH dehydrogenase II (ND2) subunit membership of wild and captive massasaugas Sistrurus catenatus. Haplotype numbers correspond to symbols in Figures 2–4. ND2 subunits of the eastern massasauga S. c. catenatus are identified as W (western ND2 subunit), C (central ND2 subunit), and E (eastern ND2 subunit). Subunit membership of two SSP Breeding Plan animals (22892, 107716) was assigned based on county of origin; haplotypes will be determined after tissue samples become available. Snakes included in the SSP Breeding Plan are indicated with a +.

Found at DOI: http://dx.doi.org/10.3996/032012-JFWM-026.S1(353 KB DOC).

We thank M. Andre, C. Berg, H. L. Gibbs, B. Jellen, R. S. King, B. Kingsbury, R. Knutson, S. Lougheed, J. Marshall, C. Phillips, M. Redmer, C. Smith, T. VanDeWall, D. Wynn, Audubon Zoo, Buffalo Zoo, Columbus Zoo, Detroit Zoo, Fort Wayne Zoo, Houston Zoo, Indianapolis Zoo, John Ball Zoo, Lincoln Park Zoo, Milwaukee County Zoo, Notebaert Museum, Oklahoma City Zoo, Potter Park Zoo, Racine Zoo, Seneca Park Zoo, Toledo Zoo, Toronto Zoo, and Wildlife Discovery Center for samples; Illinois Department of Natural Resources, Lincoln Park Zoo, and Northern Illinois University for funding; and B. Ball for preparing figures. We gratefully acknowledge input provided by J. Earnhardt, H. L. Gibbs, F. Janzen, L. Kubatko, H. Reinert, M. Sovic, two anonymous reviewers, and the Subject Editor to improve the manuscript.

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.