We collected multiple specimens of garter snakes of the genus Thamnophis from the highlands of Mexico between 2015 and 2022. We photographed all live snakes, including dorsal, lateral, and ventral profiles, and euthanized them with pentobarbital. We took tissue samples from muscle or liver upon death and preserved them in 96% ethanol. We fixed specimens in 10% formalin and transferred them to 70% ethanol for permanent storage.

All collected materials were deposited at the Instituto de Investigaciones sobre los Recursos Naturales (INIRENA) of the Universidad Michoacana de San Nicolás de Hidalgo (UMSNH) in Morelia, Mexico; the Museo de Zoología, Facultad de Ciencias (MZFC) of the Universidad Nacional Autónoma de México (UNAM) in Mexico City; and the Facultad de Estudios Superiores, Zaragoza (MZFZ) of the Universidad Nacional Autónoma de México (UNAM), also in Mexico City. Museum acronyms throughout the text follow Sabaj (2020). Specimen numbers for all material examined are provided in Appendices 1, 2. We were not able to measure type specimens of previously described taxa, so we used measurements of the type specimens provided in original descriptions and other published literature (Cope 1861, 1866; Jan 1863; Boulenger 1893; Smith 1942; Walker 1955; Rossman 1969; Rossman and Burbrink 2005). Measurements for T. errans were taken from Webb (1976), whereas measurements from Thamnophis sumichrasti Cope, 1866, and Thamnophis mendax Walker, 1955, were taken from their original descriptions.

Distribution maps were generated based on the GBIF database (www.gbif.org), which includes both museum records and distribution records from the citizen scientist platform iNaturalist (inaturalist.org). iNaturalist records were curated by us to assure that no misidentifications were present before generating the maps and are current up through January 2024. Observations that could not be positively identified were removed.

The mountains of central Jalisco have numerous other species of Thamnophis that occur in sympatry or near sympatry with the species described herein. These include Thamnophis copei (Dugés, 1879), Thamnophis cyrtopsis (Kennicott, 1860), Thamnophis eques (Reuss, 1834), Thamnophis melanogaster (Peters, 1864), and Thamnophis pulchrilatus (Cope, 1885). While these species are not closely related to the species described herein, we include them in our comparatives within the species description to aid in the identification of individuals in the field.

Abbreviations used in the text and tables are as follows: snout–vent length (SVL) , tail length (TL) , total length (TotL) , head length (HL) , head width (HW) , eye diameter (ED) , rostral height (RH) , rostral width (RW) , internasal length (INL) , internasal width (INW) , prefrontal length (PFL) , prefrontal width (PFW) , frontal length (FL) , maximum anterior frontal width (MAFW) , maximum posterior frontal width (MPFW) , parietal width (PW) , parietal length (PL) , loreal length (LL) , loreal height (LH) , mental length (ML) , mental width (MW) , anterior chin shield length (ACSL) , posterior chin shield length (PCSL) , internasal suture length (INK) , prefrontal suture length (PFK) , internasal rostral contact (INR) , rostral nasal contact (NR) , distance from frontal to snout (DFS) , muzzle length (MZL) , number of labials in contact with anterior chin shields (LCAC) , ventral scales (VS) (Dowling (1951) method and explained in Rossman et al. (1996) , subcaudal scales (SC) , dorsal scale rows at one head length behind parietal (DSRA) , dorsal scale rows at midbody (DSRM) , dorsal scale rows at one head length prior to anal scale (DSRT) , supralabials (SL) , infralabials (IL) , preoculars (PRO) , postoculars (POO) , maxilliary teeth (MT).

Scale measurements were taken in the following manner: HL, distance from the tip of the snout to the posterior border of the parietal scales; HW, distance taken at the posterior edge of the jaw; RH, the distance of the rostral scale from the median point of the mouth to the vertex formed by the internasal suture; RW, distance of the rostral scale measured between each suture formed by the 1^st supralabial and prenasal; INL, distance of the internasal scale from its anterior border with the prenasal and rostral to its posterior border with the prefrontal; INW, distance from the prenasal to the medial suture between each internasal; PFL, distance from the postnasal, nasal, and internasal borders back to the border between the frontal and supraocular; PFW, distance from the postnasal to the median suture between the prefrontals; FL maximum length of the frontal shield; distance from the posterior part of the prefrontals to the medial union between the parietals; MAFW, the maximum width of the anterior portion of the frontal scale measured between each vertex formed by the prefrontal and supraocular scales; MPFW, the maximum width of the posterior portion of the frontal scale measured between each vertex formed by the parietal and supraocular scales, PW, distance from the union of the postocular and anterior temporal scales to the median suture between the parietals; PL, distance from the border between the supraocular and frontal to the posteriormost point of each parietal scale; LL, maximum length from the upper anterior border of the nasal to the lower posterior border with the preocular and supralabial; LH, distance from the supralabial to the union with the prefrontal and preocular; ML, taken from the medial point of the mouth to the posterior end of the mental scale, where the first pair of infralabials meet each other; MW, measured along the border of the mouth from one infralabial border to the other; ACSL, distance from the suture between the first and second infralabial posteriorly to the median suture between posterior chinshields; PCSL, maximum length from the union with the anterior chinshield and infralabial to the posterior border of the posterior chinshield; SC, counted on the left and right sides, with the first subcaudal scale interpreted as the first scale posterior to the cloaca that was not counted by the anal scale; Dorsal scales were counted at one head length behind the posterior edge of the parietals, at midbody, and at one head length before the anterior border of the anal scale. Muzzle length was defined as the combined length of the internasal suture (INK) and pre-frontal suture (PFK); muzzle shape was calculated by dividing INR by NR.

All laboratory procedures were carried out at the UNAM FES-Zaragoza in Mexico City. We used a standard ammonium acetate protocol (Fetzner 1999) to extract genomic DNA from liver or muscle tissue. We then sequenced two mitochondrial loci: Cytochrome b (Cytb) and NADH dehydrogenase subunit 4 (ND4). For Cytb, we used Gludg-L (Palumbi 1996) as the forward primer and ATRCB3 (Harvey et al. 2000) as the reverse primer. For ND4, we employed ND4 and ND4_Leu (Arevalo et al. 1994) as the forward and reverse primers, respectively. Both loci were amplified using a standard polymerase chain reaction (PCR) protocol: an initial denaturation at 95 °C for 3 minutes and 30 seconds, followed by 35 cycles of denaturation at 95 °C for 30 seconds, annealing at 53 °C for 30 seconds, extension at 72 °C for 1 minute, and a final extension at 72 °C for 15 minutes, with a terminal hold at 10 °C. The PCR products were purified using a polyethylene glycol method (Lis 1980) and sequenced by Macrogen Korea (Standard-Seq of Macrogen Inc.).

Raw chromatograms were trimmed and edited using Geneious v. 2023.1 (Biomatters Ltd., Auckland, NZ). To infer phylogenetic relationships from the new samples, we included additional sequences of the genus Thamnophis as well as two outgroups obtained from GenBank. All new sequences were deposited in GenBank (Appendix 2).

Each gene was aligned separately in MAFFT version 7 (Katoh et al. 2017) using the Q-INS-I option. The alignments were then concatenated using FASconCAT v.1.04 (Kück and Longo 2014). The final alignment comprised 1828 base pairs (1104 bp for Cytb and 724 bp for ND4), including sequences from 81 representatives of the genus Thamnophis and individuals each of “Adelophis” foxi Rossman & Blaney, 1968 and “Adelophis” copei (Dugès, 1879), as well as one individual of Nerodia erythrogaster (Forster, 1771) and Tropidoclonion lineatum (Hallowell, 1856) as outgroups. We calculated pairwise genetic distances in the mitochondrial, Cytb, and ND4 genes using MEGA X software (Kumar et al. 2018).

Maximum likelihood (ML) analysis of the concatenated dataset was performed using IQ-TREE (Nguyen et al. 2015) using the IQ-TREE web server (Trifinopoulos et al. 2016). We employed an auto-substitution model and conducted 1,000 bootstrap replicates for support assessment. Separate ML analyses were also performed for each mitochondrial gene. The resulting topologies are presented in Suppl. materials 1, 2.

For Bayesian phylogenetic inference (BI), we initially used PartitionFinder v1.1.1 (Lanfear et al. 2012) to determine the most suitable model of partitions and nucleotide evolution for each locus using the Bayesian information criterion (BIC). The identified optimal partitions were: GTR + I + gamma for the first codon position of both Cytb & ND4, HKY + I + gamma for the second codon position of Cytb, HKY + gamma for the second codon position of ND4, and GTR + gamma for the third codon position of both Cytb & ND4. Our dataset was organized accordingly by locus and codon position. We then conducted BI using Mr. Bayes v3.2.2 (Ronquist et al. 2012) on the CIPRES science gateway server (Miller et al. 2011). This analysis involved four runs, each with 10 million generations and a sampling interval of 1,000 generations, incorporating three heated chains and one cold chain. Convergence was assessed using Tracer v1.6 (Rambaut et al. 2015), focusing on likelihood and parameter estimate overlaps, effective sample sizes, and the potential scale reduction factor (PSRF). Convergence was achieved within 200,000 generations, allowing us to discard the initial 25% of each run as burn-in. The results of these runs were combined using TreeAnnotator 2.7.4. (Bouckaert et al. 2019) to create a concatenated tree and visualized with FigTree v1.4.2 (Rambaut 2014).

Considering that the topologies from both ML and BI analyses were almost identical, we have included only the maximum likelihood phylogeny in this paper, with the Bayesian Inference phylogeny provided as Suppl. material 3.

In our ML phylogeny (Fig. 1), we find concordance with prior evolutionary hypotheses, including those derived from mitochondrial DNA (mtDNA) studies (de Queiroz et al. 2002), those combining mtDNA and nuclear DNA in their analyses (McVay et al. 2015, Deepak et al. 2022), as well as analyses by Hallas et al. (2022), which were based on ddRADseq data. The robustness of our phylogenetic nodes is largely high, with most nodes garnering greater than 95% bootstrap support. However, several internal nodes exhibited lower support values, so we collapsed nodes below the 50% bootstrap support threshold.

Figure 1. 

Concatenated maximum likelihood inference of the phylogenetic relationships of Thamnophis and closely related Natricines based on the mitochondrial genes Cytb and ND4.

Our analysis corroborates the delineation of three primary clades within Thamnophis, aligning with previous studies. The first clade, “Ribbon Snakes,” is composed of Thamnophis sirtalis (Linnaeus, 1758), Thamnophis proximus (Say, 1823), and Thamnophis saurita (Linnaeus, 1766), collectively forming a lineage sister to the remaining members of Thamnophis. The other sampled Thamnophis segregate into two major clades. The first of these, referenced as the “Widespread Clade” by de Queiroz et al. (2002), predominantly comprises species from the USA and northern Mexico. Our data support the close phylogenetic relationship between Thamnophis fulvus (Bocourt, 1893) and Thamnophis chrysocephalus (Cope, 1885), albeit with notable genetic divergence among T. chrysocephalus populations from Guerrero compared to those from Oaxaca and Veracruz. Further, the clade encompassing T. fulvus and T. chrysocephalus formed a sister lineage to a cluster of species endemic to the USA and northern Mexico, with Thamnophis cyrtopsis forming a sister relationship to a clade comprising species from the USA and Baja California, inclusive of “T. aff. pulchrilatus” and an individual matching the original description of Thamnophis vicinus Smith, 1942 (currently a junior synonym of T. cyrtopsis), which is notable as this population has never been studied phylogenetically.

In the other major clade, referenced as the “Mexican Clade,” our results recover an early split between T. nigronuchalis Thompson, 1957 + T. rufipunctatus (Cope, 1875), sister to all other lineages. This split is followed by the divergence of T. copei and T. melanogaster. Remarkably, this study represents the first inclusion of T. copei in a molecular phylogeny, revealing an unexpected non-sister relationship to T. foxi, which was its sole congener in the now invalid genus Adelophis (see below). The new species reported herein formed a polytomy with T. mendax, T. sumichrasti, and T. scalaris. That group was recovered as the sister clade to a group that includes Thamnophis exsul Rossman, 1969; T. errans + T. scaliger; and a group comprising T. foxi alongside T. bogerti Rossman & Burbrink, 2005; T. conanti Rossman & Burbrink, 2005; and T. lineri Rossman & Burbrink, 2005. However, it is noteworthy that the support for this latter grouping is very low, and T. bogerti, T. conanti, and T. lineri do not form monophyletic groups.

Finally, our BI analysis exhibited a high degree of similarity to those relationships obtained through maximum likelihood, with an almost identical topology. The key distinction lies in the variation of support values for certain groups. In the BI analysis, T. ahumadai was recovered as a sister to T. mendax + T. sumichrasti + T. scalaris, but with low support (posterior probability = 0.61). The BI tree, elucidating these differences, is presented in Suppl. material 3.

Our molecular results support the hypothesis that populations formerly assigned to Thamnophis scalaris from Jalisco, Mexico, belong to an undescribed species. We analyzed molecular samples from two isolated highland populations (Sierra Cacoma, Sierra de Tapalpa; Appendix 2) of “T. scalaris” from Jalisco. Both are each other’s closest relatives, with 100% bootstrap support. These populations fall within the “Mexican Clade” of our phylogenetic tree (Fig. 1). Together, these two populations are sister to a clade comprising T. scalaris, T. sumichrasti, and T. mendax, albeit with low bootstrap support (67%). The two Jalisco populations have similar genetic distances to their closest relatives. These genetic distances are 0.04–0.05 (ND4) and 0.03 (Cytb) to Thamnophis bogerti (as understood herein, see below), which appears to be their closest relative according to genetic distances. The two Jalisco populations also have genetic distances of 0.04–0.06 (ND4) and 0.03–0.04 (Cytb) from the superficially similar T. scalaris; 0.05–0.07 (ND4) and 0.04–0.05 (Cytb) from geographically proximate populations of T. errans; and 0.04–0.05 (ND4) and 0.04 (Cytb) from T. exsul. In comparison, the T. scalaris populations analyzed herein (Morelos, Estado de México, Querétaro, Puebla, and Veracruz) have intraspecific genetic distances of 0.00–0.02 (ND4) and 0.00–0.01 (Cytb). Within the T. scalaris populations analyzed, it should be noted that one specimen (AEVB 0104) from La Joya, Acajete, Veracruz had a ND4 genetic distance of 0.02–0.04 from other non-Veracuz specimens of T. scalaris. Similarly, another specimen (UOGV 3932) from nearby Nogales, Veracruz, had a higher Cytb genetic distance (0.02) from all other specimens of T. scalaris. These genetic distances are still low when compared to the interspecific genetic distances that we recovered between other species, but they are nonetheless interesting and may show evidence of hybridization with nearby populations of T. bogerti (see below).

Based on eleven genetic samples of Thamnophis scalaris collected across its range and an additional five samples of Thamnophis scaliger from several localities in central Mexico, we found that both T. scalaris and T. scaliger are monophyletic and do not represent sister species (Fig. 1). De Queiroz et al. (2002) had confusing results for the relationship between T. scalaris and T. scaliger based on the placement of their “T. scaliger 2” within T. scalaris. Their “T. scaliger 2” was based on LSUMZ 42638 (Accession Number: AF420189), which is supposedly a T. scaliger from the vicinity of Villa Victoria, Estado de México. We compared this sample to five T. scaliger samples and nine T. scalaris samples, all identified in the field by us. This sample was consistently recovered in a clade containing morphologically verified T. scalaris in all of our analyses. While we were not able to examine the specimen, we hypothesize that this genetic sample was misidentified as a T. scaliger by de Quieroz et al. (2002) and consequently re-identify it as T. scalaris.

The genus Adelophis Dugès in Cope, 1879, was originally described as a close relative of the North American genus Tropidoclonion Cope, 1860. The first authors to suggest a close relationship between Adelophis and Thamnophis were Dunn (1931) and Rossman and Blaney (1968), who suggested that Adelophis was derived from Thamnophis. Later, de Queiroz et al. (2002) obtained molecular sequences from a specimen of Adelophis foxi (LSUMZ 40846) and found that it grouped within Thamnophis, with its closest relatives being T. melanogaster and T. validus. De Queiroz et al. (2002) suggested that their sample of Adelophis may have been mixed up with another sample; however, they assessed that the genetic distances of the analyzed sequences were too far removed from any of the other analyzed species to represent any of them. These authors had samples from all valid Mexican species of Thamnophis at that time except for two (Thamnophis postremus Smith, 1942, and Thamnophis rossmani Conant, 2000) and cautioned against making any taxonomical changes pending further sampling of “A. foxi” and the similar “A. copei.” McVay et al. (2015) obtained additional sequences from the same specimen of “A. foxi” and confirmed the results obtained by de Queiroz et al. (2002), thus dispelling the possibility of PCR contamination. Nonetheless, these authors also refrained from making any taxonomical changes because the sequenced specimen was the same individual used earlier by de Queiroz et al. (2002) and was thus subject to the same caveats as the prior study. Later, Hallas et al. (2022) found Adelophis to be placed within Thamnophis across three molecular datasets (mtDNA, nDNA, and ddRADseq) and formally sank the genus Adelophis into Thamnophis. They cautioned that “some might refrain from formal changes to A. copei until that taxon can be suitably evaluated in a phylogenetic analysis.” However, Hallas et al. (2022) considered the putative sister relationship between “A. copei” and “A. foxi” (Rossman and Blaney 1968; Rossman and Wallach 1987) to be sufficient evidence refuting the argument that “A. copei” could be nested outside of Thamnophis. In this study, we sequenced samples from a specimen of Thamnophis copei for the first time and confirmed that T. copei falls within the “southern clade” of Thamnophis like T. foxi (de Queiroz et al. 2002; McVay et al. 2015; Hallas et al. 2022). This confirms that the placement of Thamnophis copei is correct and that Adelophis should be subsumed into Thamnophis. Surprisingly, however, our results suggest that the morphologically similar T. copei and T. foxi are not sisters to one another, though both are consistently nested within the “Mexican clade” of Thamnophis. More extensive sequence data is needed to determine the exact phylogenetic relationship between these two species and closely related species. It is noteworthy that T. foxi has not been collected since the 1970s, despite several recent collection attempts.

Thamnophis godmani was described from Omiltemi and “Amula” in Guerrero (Günther 1886) and was later considered to be distributed in the Sierra Madre del Sur and Sierra Madre Oriental from central Guerrero and central Veracruz (respectively) to the Isthmus of Tehuantepec (Rossman et al. 1996). Thamnophis godmani has historically been recognized as a close relative of T. scalaris, and Smith (1942) even considered T. godmani to be a subspecies of the former taxon. The specific identity of T. godmani was reaffirmed when more specimens and morphological data became available (Rossman in Varkey 1979). Later, Rossman and Burbrink (2005) conducted a multivariate study of morphological characters in various populations of T. godmani and described three species: T. bogerti, T. conanti, and T. lineri.

Rossman and Burbrink (2005) provided subtle morphological and mensural characters to justify the description of T. bogerti, T. conanti, and T. lineri as species distinct from T. godmani. While the morphological data provided does suggest that T. godmani may be specifically distinct from the other three, the morphological differences are trivial once data from T. bogerti, T. conanti, and T. lineri are compared to one another. When compared to the intraspecies variation of other widespread highland Thamnophis such as T. scaliger, T. scalaris, T. errans, T. chrysocephalus, and T. sumichrasti, the subtle geographic variation found across T. bogerti and its relatives over a relatively broad range suggests that these taxa likely represent a single evolutionary species unit. Morphological and mensural differences for the three species were presented in the description only in text; however, here we compare them in Table 2.

We evaluated the status of these species (T. bogerti, T. conanti, and T. lineri) in our phylogeny, and our results indicate that T. bogerti is paraphyletic with respect to T. conanti and T. lineri (Fig. 1, Suppl. materials 1–3). Additionally, we found very low levels of genetic divergence between these species. Genetic distances in the mitochondrial gene ND4 were as low as 0.002 between T. bogerti and T. conanti, whereas both of these species show a genetic distance of 0.04–0.05 with their closest relative, T. scalaris. Genetic distances in the Cytb between T. bogerti, T. conanti, and T. lineri ranged from 0.00–0.01, whereas all three of these “species” had distances between 0.03–0.04 from their closest relative, T. ahumadai. (Suppl. material 4). These results support the tree topologies of McVay et al. (2015) and Hallas et al. (2022). Unfortunately, we did not have any genetic material of T. godmani to compare with, and recent analyses that did include T. godmani (Hallas et al. 2022) are based on a sample of T. bogerti from Oaxaca (J. Campbell, pers. comm). Thamnophis godmani is not known to occur in Oaxaca, and all samples of “T. godmani” found on Genbank are based on specimens that were collected before T. bogerti was described (Rossman and Burbrink 2005). Consequently, our phylogenetic results indicate that T. bogerti, T. conanti, and T. lineri belong to the same evolutionary unit. As these three species were described in the same paper, we invoke Article 24.2 of the International Code of Zoological Nomenclature (ICZN 1999) and suggest that T. lineri and T. conanti are junior synonyms of T. bogerti, the latter of which was alphabetically and sequentially described first in Rossman and Burbrink (2005). The specific relationship between T. godmani and T. bogerti remains to be tested; however, the isolated nature of the Guerrero populations of this species may indeed prove to be of a specific nature.

Rossman and Burbrink (2005), in their review of the populations formerly assigned to T. godmani, stated that “T. godmani occurs in at least four discrete geographic areas that are effectively separated at the present time by habitat disjunctions unsuitable for these residents of montane pine-oak forests (1768–3048 m).” However, these sky islands inhabited by the snakes formerly assigned to T. godmani have not been proven to act as biogeographic barriers for other reptiles (Bryson et al. 2011; Palacios-Aguilar et al. 2021). Furthermore, T. bogerti is frequently collected as low as 1100–1200 m in humid environs in Oaxaca, and oak woodland in the Sierra Juárez and Sierra Zongolica ranges down to at least 1500 m. Thus, these perceived “isolated mountain ranges” are not barrier-isolated for an elevation-adapted snake like T. bogerti.

An unidentified population of Thamnophis occurs in the vicinity of Pinal de Amoles in the Sierra Gorda region of Querétaro (Rossman and Lara-Gongora 1997). These specimens (UMMZ 105415–416) are cataloged as T. scalaris but originate from well outside of the known range of this species. Rossman and Lara-Gongora (1997) discussed these specimens in detail, comparing them to T. scalaris, T. scaliger, and T. exsul, and concluded that “until fresh material becomes available, it would seem prudent to defer judgment of the identity of these Querétaro specimens.” We collected two specimens (INIRENA 2937–38, Fig. 8) of this population in 2020 and included them in our analyses.

Figure 8. 

Photographs of Thamnophis scalaris from Querétaro. Female (INIRENA 2937) (A, B) and male (INIRENA 2938) (C, D) both from 1 km N of Pinal de Amoles, Municipio de Pinal de Amoles, Querétaro.

Morphologically, the two specimens seem to match T. scalaris relatively well. INIRENA 2937 is an adult female and has 17-17-17 dorsal scale rows, whereas INIRENA 2938 is a subadult male and has 19-17-17 dorsal scales. Supralabials are seven on both specimens; infralabials are 9/10 on INIRENA 2937 and 10/10 on INIRENA 2938. Ventrals 141 on both specimens, subcaudals 52 (INIRENA 2937) and 57 (INIRENA 2938), respectively. The subcaudal count is lower than the ranges given for eastern T. scalaris by Rossman and Lara-Gongora (1997), which were 53–69 (females) and 64–85 (males). The tail lengths are short for T. scalaris, with TL/TotL 19.8% on INIRENA 2937 and 21.4% on INIRENA 2938. Rossman and Lara (1997) gave TL/TotL ranges of 20.2–25.0% (females) and 25.9–33.2% (males). Rossman and Lara (1997) noted that the hemipenes on UMMZ 104415 were nine subcaudals long. The hemipenes on INIRENA 2938 are six subcaudals long, although they may not have been everted to their fullest extent. While the tail lengths and subcaudal counts are very low for eastern T. scalaris and the hemipenes are short, in all other aspects, these specimens seem to match T. scalaris reasonably well.

Our phylogenetic analysis suggests that these specimens are nested with T. scalaris. Genetic distances for ND4 between this population and nearby Veracruz T. scalaris range from 0.00–0.01. In comparison, the genetic distances of ND4 between the Queretaro specimens and T. scaliger range from 0.02–0.06, between T. exsul 0.04 and T. errans 0.05. The genetic distance of Cytb in these specimens to nearby T. scalaris in Puebla and northern Veracruz ranges from 0.00–0.01. Comparatively, the genetic distances of Cytb in the Queretaro specimen to specimens of T. scaliger are 0.04–0.07, and for T. exsul, they are 0.04. Based on these results, we formally assign the Pinal de Amoles (Sierra Gorda) population to T. scalaris and confirm the presence of T. scalaris in Querétaro.