Abstract

As currently understood, there are two species of foxsnakes (Eastern Foxsnake, Pantherophis gloydi Conant and Western Foxsnake, P. vulpinus Baird and Girard) that are separated by a large geographic disjunction that encompasses almost all of Michigan, eastern Indiana, and eastern Ohio. Phylogenetic analysis of mtDNA data of individuals from throughout the ranges of the two species inferred reciprocally monophyletic clades that revealed a new species boundary, the Mississippi River. The single key morphological character also shows a major difference at the river. Because the localities of the holotypes of P. gloydi and P. vulpinus are both within the new range of the eastern form, gloydi is recognized as a junior synonym of vulpinus and a new name, P. ramspotti, is erected for the western form. The estimates of divergence time and historical biogeography suggest that Pleistocene glaciation and the Mississippi River played a key role in speciation.

1. Introduction

As currently understood, the Western Foxsnake, Pantherophis vulpinus Baird and Girard, inhabits a considerable portion of the Central Lowlands Province and much of the land surrounding the western Great Lakes within the continental United States [1, 2] including the upper peninsula of Michigan and west through most of Wisconsin, southern Minnesota, Iowa, southeastern South Dakota, and eastern Nebraska. The southern portion of the distribution is represented by small populations found in northern Missouri, and the rest of the range spans east through northern Illinois and western Indiana (Figure 1).

In the current concept of the Eastern Foxsnake (Pantherophis gloydi Conant), the species occupies a smaller range and is disjunct from P. vulpinus. The distribution surrounds most of Lake Erie and the eastern coast of Lake Huron and includes the southern Ontario peninsula extending west to eastern Michigan and south to north-central Ohio (Figure 1). This portion of the distribution was recently examined [3] with microsatellites and it was found that habitat degradation has played a large part in the genetic structure of these populations.

The disjunct distribution of the putative taxa Pantherophis gloydi and P. vulpinus has intrigued and vexed systematists with the distribution, morphology, and ecology all having played a role in the changing taxonomy. The Western Foxsnake was described as Scotophis vulpinus by Baird and Girard [4] based on a specimen from Racine, Wisconsin (USNM 1624) and one from Grosse Ile, Michigan (USNM 1570). The Racine snake was designated as the holotype by the action of Baird and Girard [4, viii]. Conant [1] described an eastern race of foxsnake as Elaphe vulpinus gloydi (Eastern Foxsnake). He noted that in the eastern portion of its range vulpinus exhibited a “peculiar distribution and marked habitat preferences.” In the most recent systematic studies [5, 6] the two taxa have been placed in Pantherophis. A recent paper [7] places foxsnakes in a new genus, Mintonius, in an attempt to rescue Pituophis and Pantherophis from the consequences of paraphyly. However, Pyron and Burbrink [8] inferred reciprocally monophyletic Pituophis and Pantherophis, which rendered Mintonius a junior synonym of Pantherophis.

The subspecific status of Pantherophis vulpinus gloydi remained until it was elevated to species status [9]. The taxonomic revision was argued under the guise of the evolutionary species concept (Frost and Hillis [10]) in which the two populations were allopatric (based on current distributions) and that they were “in some way morphologically (and presumably genetically) distinct” [10]. It was assumed that because of the morphological difference in blotch count of the two subspecies and that their populations were disjunct that they were independently evolving entities. This recommendation was bolstered by data presented earlier by Harding [11] who concluded that P. gloydi and P. vulpinus occupied separate ecological niches and thus should be treated as separate taxa. Harding noted that P. vulpinus inhabits prairies, pastures, farmlands, and open woodlands, basically dry habitats, where it feeds on small mammals and occasionally bird eggs and nestlings. In contrast, P. gloydi inhabits shoreline marshes, vegetated dunes, beaches, and small islands of the Great Lakes, thus occupying a more aquatic lifestyle.

Both species of foxsnakes are characterized by black to chocolate blotches with a variable ground color of yellowish to light brown. The head color usually varies from a brown to reddish tint, and the belly is yellow with a black-checkered pattern. The subcloacal plate for both species is divided, and the scales are weakly keeled.

Pantherophis vulpinus is diagnosed by the number and size of its dorsal body blotches (33 to 51; mean = 40.9) and smaller sized blotches (3-4 scales in length) while specimens of P. gloydi possess a lower count (28–39; mean = 34.5) and larger blotch size (4–6 scales in length). The number of tail blotches, maximum length, and scale counts have also been noted as minor, yet distinctive characteristics distinguishing the two lineages [1].

In the original description of Pantherophis gloydi Conant [1] expressed opposition to separating the two populations based on blotch count stating, “an intergrading population is lacking, but the considerable overlapping of the two forms in almost all characteristics is too great to admit the consideration of gloydi and vulpinus as separate species.” In recent work the two taxa have received mixed taxonomic distinction. In the SSAR Scientific and Standard English Names volumes [12, 13] they are treated as separate species as they are in Harding [11] and others [14]. But in other field guides and paleoecological records they are considered a single species divided into two subspecies (e.g., [2, 1517]).

The objective of this study is to test the species hypotheses of foxsnakes within a phylogenetic framework. Specifically, (1) are there two species lineages or (2) is there a single lineage with a large geographic disjunction within the distribution? The cytochrome 𝑏 region of mitochondrial DNA was examined and used in a phylogenetic analysis to infer relationships between Pantherophis vulpinus and P. gloydi. Phylogenies were reconstructed and used to assess current taxonomic divisions. As an operational tool, we employed a phylogenetic species concept (lineages diagnosed by apomorphies).

2. Materials and Methods

2.1. Sample Collection

Samples for both Pantherophis vulpinus and E. gloydi were obtained from areas throughout both species respective ranges. The samples varied in type, consisting of blood, liver tissue, muscle tissue, shed skins, and scale clips. All muscle tissue, liver tissue, and scale clip samples were stored in 70% ethanol or kept frozen at −60°C. Blood samples were stored in a lysis buffer, while shed skins were kept in zip-lock bags until DNA was extracted. Additional samples of the Western Foxsnake were personally collected (permit information available upon request). Each captured individual was scale-clipped on one of the ventral scales, cleaned with ethanol, and released back into the environment. Institutional abbreviations are as listed at http://www.asih.org/codons.pdf.

The specimens collected from the Eastern Foxsnake distribution consisted of nine individuals from four populations within Ohio and the Ontario peninsula region. The Western Foxsnake specimens included individuals and populations from Nebraska, Minnesota, Iowa, Wisconsin, Illinois, and Missouri. Table 1 provides a list of all the samples used in this study including sample identification, locality data (see Figure 1 for map of localities), and collection information.

2.2. Molecular Methods

Total DNA was isolated from tissues (muscle, liver, blood, shed skins, and scale clips), using a standard phenol-chloroform protocol, a DNeasy kit (Qiagen), or Invitrogen’s “Easy DNA” kit. All DNA isolations were tested on a 1% agarose gel to visualize the sample’s degree of yield prior to PCR amplification.

The cytochrome 𝑏 gene was amplified using Advantage cDNA PCR kits (Clontech Laboratories) and universal primers [18]. Cytochrome 𝑏 primers: L14910: GACCTGTGATMTGAAAACCAYCGTTGT; H16064: CTTTGGTTTACAAGAACAATGCTTTA.

PCR amplification was conducted on a Perkin-Elmer thermocycler using the standard cycle conditions for the cytochrome 𝑏 primers. Negative controls were used to check for contamination of experimental samples.

PCR amplifications were visualized and tested for correct size of fragment on a 1% agarose gel. Several of the cytochrome 𝑏 PCR amplifications produced multiple sized bands ranging from ~1200 bp–300 bp. To solve this problem gel extractions were performed using a Gel Extraction kit (Qiagen) to isolate the desired PCR fragment. New cytochrome 𝑏 primers were also designed to avoid multiple bands (Qiagen). The specific primers were designed based on the GenBank sequence (Accession no. AF138758) from an individual of Pantherophis gloydi [18]. Cytochrome 𝑏 primer (designed, P. gloydi): Cytb2R: ATGCCCAACCAGCACATACTTC; Cytb2F: GTTATTGAGATTTTATTTTCGAG.

The PCR products were ligated into TopoTA vectors and cloned using the TopoTA cloning kit (Invitrogen). The products were transformed into TOP 10 competent E. coli cells and grown on LB plates at 37°C for 12+ hours. X-Gal was applied to the plates for proper colony selection. Colonies were picked from plates and grown in LB with 50 μg/ml kanamycin for more than 14 hours. Plasmid DNA was isolated from the competent E. coli cells using PerfectPrep Plasmid DNA kits (Eppendorf 5 Prime), cut with EcoRI restriction enzymes, and tested on a 1% agarose gel to verify the plasmid containing the desired insert. DNA was sequenced using the Sequitherm Excel II DNA Sequencing Kits (Epicentre Technologies) with M13-Forward and M13-Reverse primers. Advantage-2 DNA polymerase was substituted for the Sequitherm Excel kit polymerase to improve sequencing length and reading efficiency. The sequencing reactions were performed using cycle sequencing and loaded on a Licor 4000L for automated sequencing.

2.3. Alignment

Thirty-two cytochrome 𝑏 sequences (25 Pantherophis vulpinus individuals and eight P. gloydi individuals) were pre-edited on E-Seq version 2.0 to reduce automated sequencing errors. Forward and reverse sequences were then transferred to Sequencher 3.1 as Sqv. 2 (Sequencher v. 2) reports and imported as sequence files into the program. Contigs of the 33 foxsnake cytochrome 𝑏 sequences were created based on individual overlap between forward and reverse sequence fragments. The aligned data set was exported as an NEXUS file into PAUP* [19].

2.4. Outgroup

The phylogenetic placement of foxsnakes within Pantherophis may be with P. obsoletus [18, 20] although later Burbrink [21] investigated the phylogeography of the cornsnake and suggested that P. guttatus was the sister to the foxsnakes. As such, both P. obsoletus and P. guttatus were used as outgroups in the phylogenetic analyses.

2.5. Phylogenetic Analysis

The data were analyzed within parsimony (MP) and likelihood (ML) frameworks with PAUP*. The MP analyses were implemented with all sites equally weighted and gaps treated as missing data. The heuristic search was conducted using 10,000 random addition searches with global branch swapping (Tree Bisection-Reconnection). Each random addition search began with a random tree and the fit of the data to the tree was evaluated with the descriptive statistics consistency index [CI; [22]], retention index, and rescaled consistency index [RI, RCI; [23]]. The robustness of the hypothesis was estimated using the decay index [24] and a nonparametric bootstrap (1000 replicates) [25]. Character states were optimized under the conditions of delayed transformation. Under the maximum likelihood optimality criterion, the most appropriate model of evolution was selected using jModeltest v0.1.1 [26, 27] with the hierarchical likelihood ratio test (hLRT) and the Akaike Information Criterion (AIC; [28]). Under hLRT the HKY model was selected, and under the AIC the HKY+G was selected. Given the arguments made in Posada and Buckley [29] with regard to hLRT and AIC model selection, the data were analyzed under HKY+G parameters.

2.6. Timing of Divergence

We used jModeltest v0.1.1 [26, 27] with the AIC to estimate the most appropriate models of evolution and prior parameters for our analyses. The rate or substitution and date of gene divergence were estimated using BEAST v1.4.8 [30] on the BioHPC v1.4.8 parallel computing cluster at the Cornell University Computational Biology Service Unit (http://biohpc.org/default.aspx). All BEAST analyses were run for 107 Markov Chain Monte Carlo (MCMC) chains that were sampled every 1000 iterations. Log files were analyzed in Tracer v1.4.1 [31], and effective sample sizes (ESSs) were used to evaluate the estimates of posterior distributions. The first 10% of all analyses were discarded for burn-in. A summary of the output trees was generated with TreeAnnotator v1.4.8 [30], discarding the first 1000 trees for burn-in, and then analyzed with FigTree v1.2.1 [32].

To estimate the rate of substitution of cytochrome 𝑏 in foxsnakes, we aligned our sequence data (including P. guttatus) with Bogertophis rosaliae and Pseudoelaphe flavirufa sequences gleaned from GenBank (DQ902102 and DQ902109, resp.). Sequences were aligned by eye and substitution rate was calibrated using the earliest known North American Pantherophis fossil (16 mya; [6, 17]). We designated this as the split between Pantherophis and non-Pantherophis and used a LogNormal prior with a standard deviation of 1. All BEAST analyses were run using the Yule Process tree prior.

Divergence date was estimated using a strict clock model because preliminary analyses revealed clock-like evolution (ucld.stdev: 0.537; [33]). Simulations were run with the constant size prior, and two substitution rates were tested: the universal cytochrome 𝑏 rate (1% per my; [34]), and the mean substitution rate calibrated in BEAST. The analysis yielding the highest posterior probability (PP) density was then selected, and three additional runs were conducted. All iterations were then pooled using Logcombiner v1.4.8 (Drummond and Rambaut, 2007).

3. Results

3.1. Sequence Data

The cytochrome 𝑏 region of mtDNA in both species of foxsnake was composed of a 1116 bp region, but was trimmed to 1108 bp for alignment (Genbank accession numbers FJ267654–FJ267687; aligned data set available upon request). The eight base pairs were cut from the beginning of the sequence, which contained the start codon ATG and were followed by TGTTC. This region was conserved between Pantherophis vulpinus and P. gloydi specimens and did not contain any phylogenetically informative characters that would otherwise influence relationships between the two species. The cytochrome 𝑏 region of foxsnakes is similar in size (1116 bp) and structure to that of its relatives, P. obsoletus and P. guttatus [18] and other colubroid snake species [35].

The nucleotide composition of the region in foxsnakes is dominated primarily by adenine and thymine, with a low frequency of guanine. The base frequencies for thymine (T) ranged from 0.305 to 0.308 with a mean of 0.3065 while adenine (A) reached slightly higher frequencies of 0.327–0.331 with a mean of 0.329. The frequencies for cytosine (C) ranged from 0.258 to 0.261 at a mean of 0.2595 and guanine, the lowest, 0.104–0.106 (mean = 0.1048). The high frequencies of adenine and thymine are comparable to the outgroup taxa P. obsoletus (accession no. AF28361) which have base frequencies of 0.293 for As and 0.329 for Ts. The low frequency of guanine was also observed in the P. obsoletus sequence (G = 0.111).

3.2. Phylogenetic Relationships

The maximum parsimony and maximum likelihood analyses inferred the exact same clades so only the parsimony tree will be discussed further. The MP analysis yielded six most-parsimonious trees (strict consensus, Figure 2) of 222 steps, and these were found on a single island hit 10,000 times. There were 143 variable, but uninformative characters, 60 phylogenetically informative characters, and 905 constant characters. The consistency index was 0.99 when all characters were included in the analysis. When uninformative sites were excluded, the CI = 0.96, RI = 0.98, RCI = 0.97.

Two distinct clades within the foxsnakes were identified; (1) an eastern clade, comprised of most foxsnakes from east of the Mississippi River, including eastern Pantherophis vulpinus and all P. gloydi individuals (95% bootstrap proportion; decay index = 3) and (2) a western clade, comprised of most foxsnakes from west of the Mississippi River (100% bootstrap proportion; decay index = 7). Five individuals (E11a, E9b, E8a, E8b, and W5b) were identified whose haplotypes were found in opposite clades. An individual (W5b) from Muscatine County, 1A, geographically located west of the Mississippi River was grouped with foxsnakes in the Western Clade. Individuals from southern Wisconsin (E11a) and Western Wisconsin (E8a, E8b, and E9b) located east of the Mississippi River were grouped with foxsnakes in the Western Clade.

There were a total of 11 distinct haplotypes found in the cytochrome 𝑏 gene of the 33 ingroup foxsnake specimens used in the analysis. The haplotypes are shared between the following taxa and localities: (A) Buffalo County, WI (E9b), Pepin County, WI (E8a and E8b); (B) Rock County, MN (W2b), Warren County, IA (W4), Muscatine County, IA (W5a) (C) Holt County, MO (W3b), Dakota County, NB (W1); (D) Erie County, OH (E18), Ottawa County, OH (E21); (E) Dekalb County, IL (E16), Haldimond-Norfolk County, Ontario (E20b); (F) Green County, WI (E12a), Iowa County, WI (E10); (G) Putnam County, IL (E13), Essex County, Ontario (E19), Haldimond-Norfold County, Ontario (E20a, E20c, E20d, and E20e), Green County, WI (E12b), Lasalle County, IL (E14), Kane County, IL (E15), Muscatine County, IA (W5b). Unique haplotypes were found in Wabasha Co., MN (W6), Vilas Co., WI (E7), Rock Co., WI (E11a, 𝑏 ), and Dewitt Co., IL (E17) (Figure 3).

The pairwise genetic distances between the two (eastern versus western) major clades, as well as intraclade distances, were characterized based on individual comparisons (Table 2). The interclade comparisons of genetic distances ranged from 0.903 to 1.0444% and were much greater than intraclade distances: 0–0.271% (eastern clade) and 0–0.361% (western clade).

The results for sequence divergence at first, second, and third codon positions indicated an expected pattern of substitution rate. At first codon positions there were two informative, 15 variable, and 353 constant characters. Second condon positions included two informative, 12 variable, and 355 constant characters. Third positions contained 11 of the 15 informative characters with 93 variable and 265 constant. In the interclade comparison of the most diverged sequences 16 sites of the 1108 were different with a transition-transversion ratio of 3.0 (12/4). Ten of the divergent characters were in third codon postions with a ts/tv ratio of 4.0 (8/2).

3.3. Divergence Times

The best-fit model of evolution estimated by jModeltest was GTR + I. In all BEAST analyses, ESSs were high (>300), indicating increased reliability estimates. Calibration of substitution rate resulted in a mean substitution rate of 0.58% per million years (95% highest posterior density range: 0.28%–0.93%). The universal cytochrome 𝑏 rate yielded a higher PP than the calibrated rate (−1712.704 and −1730.983, resp.) and was thus selected for subsequent iterations. Pooled estimates of divergence time between Pantherophis gloydi and P. vulpinus resulted in a most recent common ancestor between 369,100 and 1,183,000 ya (mean: 742,800 ya; PP = −1712.527; ln L = −179.573).

4. Discussion

4.1. Phylogenetic Inference

The results of the phylogenetic analyses indicate the presence of two distinct clades. The Eastern Clade is diagnosed by four unambiguous synapomorphies and contains most of the individuals of Pantherophis vulpinus from east of the Mississippi River and all the individuals of P. gloydi from populations in Ontario, Ohio, and two Lake Erie islands (Kelleys Island, OH and Pelee Island, Ontario). The Western Clade is diagnosed by seven unambiguous synapomorphies (Table 3) and is composed mostly of individuals of P. vulpinus from west of the Mississippi River (Figure 2).

The two clades diagnosed in the analyses do not correspond to the currently accepted concepts of foxsnake species nor with the currently recognized distributions. The hypothesis of two widely disjunct taxa is rejected. In addition, there is no evidence that foxsnakes form a single species represented by a single lineage. In contrast, the results indicate two genetically distinct lineages (species) and reveal historical relationships previously unrecognized.

Five of the 33 individuals were found in clades that did not support an east-west separation at the Mississippi River. Individuals (E9b, E8a, E8b, and E7) of Pantherophis vulpinus were exceptions to the general pattern in that although they are distributed east of the Mississippi their haplotypes are related to members of the Western Clade. Another individual (W5b) exhibited the opposite, appearing in the Eastern Clade but was found west of the river (Figure 2). We hypothesize that these five individuals indicate that the barrier has been in the past or is currently leaky. Leaky refers to migration of either form across the barrier, but different processes may account for the distributions of these haplotypes. One is hybridization, although there have been no reports of hybrid individuals. A second process may be incomplete lineage sorting leading to the retention of an ancestral polymorphism and yet another may be trapping of individuals because of river meandering. While either of these may explain the location of the different haplotypes across the barrier, neither of them can reject the hypothesis that two species are present. The ability to interbreed is plesiomorphic and thus not surprising across species boundaries (e.g., [36]), and incomplete lineage sorting represents marks of a historical connection but not the rejection of lineages evolving on unique trajectories.

Some individuals of P. vulpinus found in separate clades, yet that were geographically close, exhibited high genetic divergence. For example, an individual from Wabasha County, Minnesota (W6-Western Clade), compared to an individual (E9a-Western Clade) from just 50 km east in Buffalo County, Wisconsin, exhibited one of the largest pairwise distances in the study (1.1%). This is in comparison to the low genetic difference (zero) exhibited between an individual from a population in Ontario, Canada (E20c) and one from Muscatine, Iowa (W5b), 1000 km away.

Harding [11] recommended that the populations on either side of the disjunction should be considered separate taxa because they are ecologically distinct. The results of the present study suggest that although the populations east of the disjunction exhibit different ecology, the ecology alone has not been sufficient to move them in genetically distinct trajectories. Although in Van Valen’s [37] original formulation of the ecological species concept he claimed that ecology keeps species apart more than reproduction, it appears that with foxsnakes that has not been the case. Or, it has not been the case long enough.

4.2. Divergence Time Estimation and Speciation

The divergence time estimates between the Eastern and Western Clades yielded a range of 369,100–1,183,000 ya (mean = 742,800 ya) which nearly encompasses the entire set of Pleistocene glaciation events in North America [38, 39]. With the increasing sophistication of analyses of molecular phylogeographic and phylogenetic data the debate on the importance of Pleistocene speciation [40] has increased, with some studies rejecting the hypothesis (e.g., [41]), and claiming speciation rates were no different during the Pleistocene relative to earlier periods in the Tertiary [42]. More studies corroborate the hypothesis that Pleistocene glaciation events played a significant role in speciation (e.g., [4348]); however Avise et al. [44] found few herpetofaunal pairs in support of the hypothesis. Pyron and Burbrink [8] examined divergence times among lampropeltine snakes and found that most divergences predated the Pleistocene except for two species pairs of Pantherophis. For the species pair of Pantherophis in this study the divergence time estimates strongly suggest Pleistocene glaciation as the driver of speciation.

Given estimated divergence times, the two clades as distinct lineages, and the patterns of genetic divergence relative to geographical space, we infer foxsnakes underwent allopatric speciation with the Pleistocene glacial lobes and the Mississippi River acting as vicariant barriers. There is evidence indicating that the Mississippi River has played an important role in the evolutionary history and speciation of several closely related species in North America (e.g., [18, 21, 4954]). In many of these studies species boundaries correspond to the lower Mississippi embayment during the Pleistocene. This is the first study on snakes providing evidence that the upper Mississippi River (but south of the glaciers) served as a similar geographic barrier in the evolution of the more northern foxsnake species. During the interglacial periods rivers underwent alluviation, carving deep trenches into river beds, and this caused the width of rivers, including the Mississippi, to increase beyond current flood plains [55], making the Mississippi a more formidable barrier to biotic expansion than it is currently. Typically, where the river appears to act as barrier, the distributions of species are disjunct. Foxsnakes and other species of Pantherophis are unusual cases in that their range is continuous across the Mississippi River. This sympatry suggests that the Mississippi River was an active boundary to gene flow and may have played a role in the past but perhaps no longer acts as an effective barrier (at least as one possible explanation for the sympatry of mtDNA haplotypes). Thus we conclude that it was the tandem of barriers, glaciers, and the Mississippi River that sundered a once continuous distribution and created conditions for allopatric speciation.

The presence of four western haplotypes east of the Mississippi River indicates dispersal across or around the river in Minnesota and northern Wisconsin. The presence of a single eastern haplotype in Iowa west of the Mississippi River indicates that dispersal also has occurred across the river in the opposite direction. Again, these data suggest that the Mississippi River may be no longer an absolute barrier to dispersal. The extent of dispersal affecting gene flow has yet to be investigated, but it is assumed that some genetic exchange is occurring (based on specimens from Buffalo Co., Wisconsin that occur in opposite clades). It is worth noting that it is well documented that mtDNA haplotypes cross what appear to be established species boundaries [36], so at minimum it appears that female genomes are distributed across the barrier.

An alternative hypothesis to leakage (dispersal) across the barrier is that the river is a functioning barrier and the distribution of western haplotypes east of the river represents an old capture event. Rivers change course over time, and the meander could have captured some individuals on the “wrong” side. Without recombination of the mtDNA, these haplotypes would persist regardless of hybridization with the eastern populations.

The DNA and the geography indicate two lineages, but what about the morphology ? Conant [1] regarded blotch count as the significant character diagnosing his two lineages. Interestingly, the blotch count difference between the two lineages was a mean of 6.5. If the geographic boundary between the lineages is moved west to the Mississippi River, the difference is still 6.0, essentially the same value considered diagnostic by Conant [1, page 7], “The chief diagnostic characteristic of the two subspecies of vulpinus is in the number of dorsal blotches, counting from the head to a point directly above the anus.” With the shift in the distribution of the species boundary west, dorsal body blotch count remains the key diagnostic morphological trait.

4.3. Taxonomy

The recognition of new species boundaries (Figure 4) for these lineages creates taxonomic problems. The type specimen (ANSP 21650) for Pantherophis gloydi is from Lucas County, Ohio. Conant [1] designated a neotype for P. vulpinus from Racine, Wisconsin (USNM 9969). Both of these localities are within the distribution of the eastern lineage, and because vulpinus is the older name, the eastern lineage takes the specific epithet P. vulpinus. Pantherophis gloydi is reduced to synonymy. This action leaves the western lineage unnamed. A third name is possibly available, Elaphis rubriceps Duméril et al. [56], from “l’Amerique du Nord”. This name has been treated as a synonym of vulpinus (e.g., [57]). Conant [1] concluded that the type was lost, and although the description of rubriceps appears to be of a specimen of P. vulpinus, he could not determine its provenance and did not use the name. In the original description the authors refer to Castelnau as the origin of the type specimen. The Alexandre Westphal-Castelnau collection was originally housed in Montpellier, France but since moved to the Muséum National d'Histoire Naturelle in Paris. However, in the 1870 catalogue of that collection there is no listing of Elaphis rubriceps. Although some of the Castelnau collection remains unregistered in the Paris museum, “it seems unlikely that the type could be in that collection and not figuring in the catalogue of 1870 from that collection. Thus the more likely answer to your question is that the type is lost, destroyed, or still in the unregistered Castelnau collection but that hypothesis seems unlikely.” (Ivan Ineich, pers. comm. in agreement with Conant [1]). Given this information, we regard Elaphis rubriceps as an unavailable name for the western form. For the western form we propose the name:

4.3.1. Pantherophis ramspotti, sp. nov.

(Western Foxsnake)

Holotype
National Museum of Natural History, USNM 578514 (original number JRP 1395), adult female gravid with seven eggs collected in Warren Co. (R24W, T76N, Sec. 30), Iowa, 21 June 2002 by Jeff Parmelee.

Paratypes
National Museum of Natural History, USNM 578515 (original number JRP 1449), Marion Co. (R21W, T75N, Sec. 25), Iowa, 5 June 05; USNM 578516 (original number JRP 1447), Warren Co. (R22W, T76N, Sec. 23), Iowa, 5 June 05; USNM 578517 (original number JRP 1448), Warren Co. (R22W, T76N, Sec. 21), Iowa, 5 June 05; USNM 578518 (original number JRP 1450), Warren Co. (R24W, T76N, Sec. 25), Iowa, 18 May 05; SLU 1036 (original number JRP 1558), Madison Co. (R27W, T77N, Sec. 36), Iowa, 15 September 02; USNM 578519 (original number JRP 1555), Warren Co. (R24W, T77N, Sec. 33), Iowa, 27 May 03.

Diagnosis
A relatively large Pantherophis, most similar and closely related to P. vulpinus. P. ramspotti differs in blotch count, distribution, and cyt- 𝑏 nucleotides. P. ramspotti has an average dorsal body blotch count anterior to the cloaca of 42.8 whereas P. vulpinus has an average of 36.8 dorsal body blotches. P. ramspotti is generally distributed west of the Mississippi River and possesses seven unique nucleotides at positions 121 (T), 191 (G), 284 (C), 679 (A), 689 (G), 707 (C), and 830 (C).

Description of the Type Specimen
Rostral wider than high, notched below. Nasal divided with nasal opening approximately centered between the two parts. Loreal wider than high in a somewhat trapezoidal shape. One preocular, twice a high as wide. Two post oculars. Temporals right: 2-3-3, left 2-3-4 with the most dorsal scale in the second row fused with the corresponding scale in the third row and in the third row a scale is divided making the third row with 4 scales but in a 1-2-1 dorsal to ventral pattern. Internasals two, one-half to one-third the size of the two prefrontals. Frontal higher than wide, broadening anteriorly and forming a pointed tip posteriorly where it contacts the parietals. Supraoculars longer than wide and at their widest less than twice as long as wide. Paired parietals smooth, less than one and a half times longer than the frontal, regular, and symmetrical posteriorly. Eight supralabials, 4+5 in contact with the orbit, and the seventh the largest. Infralabials 11 on each side, with 5 the largest on the right and 6 the largest on the left. One pair anterior chin shields in contact and larger than posterior pair which are separated by pair of smaller scales anteriorly and a row of 4 smaller scales posteriorly. Dorsal scale rows 27-25-21. First 9 rows smooth (opposite seventh ventral) grading from weakly keeled to moderately keeled beyond the cloaca. Dorsal scales with 2 apical pits at tips. Ventrals 206. Anal plate divided. Subcaudals 56 in two rows terminated by a long pointed claw-like tip. Dorsal ground color of specimen (in alcohol, six years old) light olive brown anteriorly grading darker towards cloaca and becoming lighter again towards tail tip. Color also grades to lighter from dorsal midline laterally towards ventrals with the two most ventral lateral scales distinctly lighter. The first six cm behind the head gives the appearance of narrow light stripping laterally. Scales variably stippled in black with perimeter appearing cream over anterior third of body then becoming brown. Some scales along dorsal midline appear discolored white. Dorsal blotches dark brown with perimeter often black. Dorsal blotches 44, with two irregular, one smaller than normal, and three smallish blotches separated by 1-2 scales counted as one. Tail blotches 13. Top of head hazel-olive brown with dark band across posterior half of prefrontals reaching to the dorsal anterior portion of the eyes. Another band extends from posterior of the eyes at an angle terminating at the two rear supralabials. On the right side a dark band extends from below the eye on supralabials 4+5, on the left the band extends across 2-5. Parietals and frontal with distinct dark markings. Looking from the anterior, the frontal has a broad triangle with the point anterior. The parietals have a broad “M” shape and two spots in contact with the frontal at the two posterior tips of the triangle. There is also a short dark stripe along the medial line of the parietal scale contact posteriorly. The upper and lower labials are yellow-cream with the sutures tinged with dusky gray. Chin and throat yellow-cream. Temporals same as basic head color but becoming darker posteriorly with color extending into a dorsal blotch outlined in black on each side. Venter yellow-cream with 4 rows of black blotches, with most lateral blotches extending into the first 2-3 scale rows and either alternate with dorsolateral blotches or connect with them (on the right side 31 do not and 28 connect, on the left side 2 do not and 31 connect). SVL = 74.5 cm, TL = 13.4.

Paratype Variation
Temporals range from 2-3-4 to 2-4-5. One specimen with left temporal arrangement of 2-4-1-5 (SLU 1036). Infralabials sometimes 10. Chin shields variably divided by 1-3 rows of scales. Dorsal scale rows 25-25-21 and 24-25-21 dorsal keeling variably starts 8–11 scales behind head. Ventrals range from 191 to 208 and subcaudals 60 to 65. Dorsal blotches 39-46 and tail blotches 17–19.

Range
West of the Mississippi River in southeastern Minnesota, most of Iowa, southeastern corner of South Dakota, eastern Nebraska, northwestern and northeastern corners of east-central Missouri.

Etymology
The specific epithet ramspotti is in remembrance of the late aspiring herpetologist Joseph Ramspott, who was from Nebraska, a state included in the western form’s range. He was a graduate student at Southeastern Louisiana University working in the Crother-White lab when he passed away in 2004.

Nomenclatural Comment on P. vulpinus
In the course of the present study, the issue arose as to which specimen should serve as the name bearing type of P. vulpinus. Conant [1] concluded that the Baird and Girard holotype (USNM 1624) was lost and designated a neotype (USNM 9969). When one of us (BIC) requested to examine the Baird and Girard specimens it was pointed out that there was a confusing notation in the catalogue for USNM 7269,

“Ledger has “Type B&G Cat. Serp. P. 75” in remarks but it is crossed out, and “topotype” has been added. All are original or old annotations”.


We examined USNM 7269 and compared it to Baird and Girard’s original type description and found that it did not match. Two other topotypic specimens from the original series from Racine were sent one each to the Natural History Museum, London (USNM 1625) and the Hungarian Natural History Museum, Budapest (USNM 1626), and one was possibly the missing holotype. The Hungary collection was destroyed in the revolt of 1956 but after examination the specimen in London (BMNH 1861.8.15.41) was determined not to be the missing holotype. Therefore the name bearing specimen for P. vulpinus is Conant’s neotype, and USNM 7269 is not the type but a topotype as noted in the ledger.

4.3.2. Pantherophis vulpinus Baird and Girard

(Eastern Foxsnake)

Diagnosis
As described in [1].
Range
All populations east of the Mississippi River including Wisconsin, Illinois, eastern Missouri, Indiana, Michigan, Ohio, and Ontario, Canada.

Remarks
As noted above, with the new understanding of the distribution of these two forms of foxsnakes, the holotype for the name P. vulpinus (USNM 9969) remains the same, although now the name applies to the populations of foxsnakes from east of the Mississippi River and to the populations formerly included in P. gloydi.

Historical Biogeography
Although there is some debate (as previously noted), much of the current distribution of vertebrates in North America, especially the herpetofauna, has been in some way shaped by the glaciation events during the Pleistocene [17]. The Pleistocene is characterized by several phases, indicated by glacial and interglacial phases. During the glacial phases, ice sheets extended as far as southern Illinois, Indiana, and Ohio during the Illinoian phase and through the central regions of those states in the Wisconsinan phase (Figure 5). The climate during the glacial maxima was apparently cold in the north but more moderate in the central and southern portions of the United States. The glaciations were detrimental to many vertebrate species, and one would expect ectotherms especially, but in fact while families, genera, and species of mammals went extinct, no families, no genera, and only two species of snakes are known to have become extinct during glaciation [17]. Instead snakes, including foxsnakes, appeared to make significant range adjustments.
We concur with Holman [17] that foxsnakes, possibly as part of a continuous distribution of a single species, exhibited ice age population adjustments and retreated to southern refugia during periods of glacial extent. The southern refugia hypothesis has been proposed for several amphibian and reptile taxa and involves the southward movement of northern populations during glacial advances (e.g., [17, 58, 59]). The fossil evidence of foxsnakes corroborates this hypothesis providing records that indicate the presence of P. vulpinus (sensu lato [1]) in areas throughout southeastern and south central United States during the Pleistocene [16, 17]. However, the fossil record also shows that foxsnakes were not only south, but also had northerly distributions east and west of the main ice sheets (Figure 5).
After the last ice age at the end of the Wisconsinan, approximately 20,000 ybp, foxsnakes were divided by the Mississippi River. Presumably, speciation was well underway by then because the eastern and western distributions were separated by the ice sheets and the Mississippi River. The estimated dates and the low genetic divergence (0.903–1.444%) between the clades east and west of the river support the recent separation. The postglacial distribution of foxsnakes expanded north, east, and west into the current northern distributions. Holman [60] categorized foxsnakes as secondary invaders (species able to exist in mixed coniferous-broadleaf forests) in the reoccupation of formerly glaciated regions. The pattern of reoccupation east of the Mississippi River is similar to that of Ambystoma maculatum (spotted Salamander), in that there are two Midwestern tracks, one up Indiana, Illinois and west of Michigan and the other eastward through eastern Ohio, Michigan, and into Ontario [61]. It is clear from the phylogenetic pattern and the extremely low genetic distances that these eastern tracks remained but did not diverge to east and west of Michigan until late. At minimum we can say that foxsnakes recolonized the once glaciated north through multiple fronts like other taxa (e.g., [62, 63]). However, foxsnake recolonization was different because their current distribution is only in previously glaciated regions. The current distribution overlain with the maximum extents of glaciation gives the distinct impression that foxsnakes simply filled in the new habitat made available by the retreat of the glaciers (Figure 5). We hypothesize that by the time of the recolonization two species had already evolved.
A rapid northward expansion by foxsnakes is suggested by the low intraclade genetic variability. This pattern of low genetic diversity with increasing latitude has been observed in plants [64], insects [65], salamanders [66], mammals [67], other animals [68], and more recently in snakes ([69], but see [63]) for Thamnophis sirtalis) and turtles [70]. Further fossil evidence of P. vulpinus in Saginaw Bay, Michigan approximately 4000–1000 ybp [71] is in agreement with the idea that the northward expansion must have been rapid, occurring between 10,000 ybp (end of Wisconsinan) and 1000 ypb. Currently, the geographic distribution of foxsnakes is represented by two species, Pantherophis vulpinus and P. ramspotti. The lack of genetic variability within the disjunct populations of P. vulpinus poses two hypotheses: (1) P. vulpinus occurs in the intermediate area within eastern Indiana, western Ohio, and central Michigan but has not been documented, or (2) the geographic disjunction within P. vulpinus is real, but has occurred so recently that the populations have yet to diverge (e.g., as in Nerodia floridana, [71]). Herpetologists have collected extensively in the intermediate area over the last 100 years, and no foxsnakes have been found. Therefore, it appears that the disjunction is a real and a recent event, possibly because of the expansion of agriculture and other human impacts (Holman, pers. comm.).

5. Conclusions

The molecular data employed in this study diagnose two clades and indicate that the previously accepted taxonomy is inconsistent with the geographic distribution of the clades. One clade, Pantherophis vulpinus (Eastern Foxsnake), consists of foxsnakes found east of the Mississippi River including western Indiana, Illinois, Wisconsin, portions of eastern Minnesota, the upper peninsula of Michigan, and the disjunct eastern distribution of northern Ohio, eastern Michigan (lower peninsula), and Ontario surrounding Lake Erie. The other clade, Pantherophis ramspotti (Western Foxsnake), is distributed west of the Mississippi River in Iowa, eastern Nebraska, South Dakota, Missouri, and portions of Minnesota west of the Mississippi River. We anticipate that further data will corroborate the Mississippi River as the species boundary and the specific epithets of foxsnakes will settle into stability.

Acknowledgments

Numerous people contributed to this study. F. Burbrink, F. Durbian, B. Hay, R. King, J. LeClere, E. McCumber, J. Parmelee, C. Phillips, and R. Siegel all kindly contributed specimens/tissues. I. Ineich (MNHN), L. Ford, and D. Kizirian (AMNH) are thanked for generously answering museum queries. Roy McDiarmid, Steve Gotte and James Poindexter (USNM), and Colin McCarthy (BMNH) are thanked for helping with specimen loans. J. Parmelee kindly contributed specimens for the type series of ramspotti. L. Andrews and M. Brecheen graciously helped rework certain sections and tables. F. Burbrink commented on the paper. Collection permits: Nebraska no. 226, Minnesota no. 10458, and Illinois to David W. Gardner. This work was partially funded by Southeastern Louisiana University’s faculty development grant program and by NSF DEB-9207751.