Phylogeography of the Rough Greensnake, Opheodrys aestivus (Squamata: Colubridae), Using Multilocus Sanger Sequence and Genomic ddRADseq Data

Abstract. The Rough Greensnake, Opheodrys aestivus, is a moderately sized, semi-arboreal snake broadly distributed throughout eastern North America. Although numerous taxa with similar distributions have been shown to comprise multiple species, O. aestivus has not been examined in a detailed phylogeographic context. Here, we use Sanger-sequence data (one mitochondrial, three nuclear loci) for samples from throughout the distribution of O. aestivus to elucidate phylogeographic patterns in this species. We combine this with ddRADseq data for a subset of samples to test patterns on a more genomically comprehensive scale. In both data sets, we find strong support for three deeply divergent clades within O. aestivus: peninsular Florida, central Texas, and a Main clade comprising the rest of the distribution. Estimates of divergence time suggest that the central Texas and Main clades diverged ∼1.3 million years ago (mya), while the peninsular Florida clade diverged from other lineages ∼2.9 mya, and these lineages diverged from the sister taxon, O. vernalis, ∼6.4 mya. The divergence of peninsular Florida or central Texas populations is not unique among squamates, nor is low levels of divergence from the Atlantic coast to eastern Texas, but this combination of patterns is unusual and yields important insight into the biogeography of North American biota. Further, our approach helps illustrate how dense geographic sampling with limited genomic sequencing can be used as a guide for the selection of samples to test phylogeographic patterns comprehensively. We conclude that elevating O. a. carinatus to species status may better describe the diversity of this genus.

Assessing phylogeographic structure within species is a critical aspect for making robust, fine-scale inferences regarding the evolutionary history of populations and taxa (Avise, 2000;Hickerson et al., 2010;Riddle, 2016;Rissler, 2016;Kumar and Kumar, 2018). Failure to account for such structure can bias downstream analyses and inferences, and may also impede conservation efforts (Coyne and Orr, 2004;Isaac et al., 2004;Bickford et al., 2007;Edwards et al., 2007;Bortolus, 2008;Hickerson et al., 2010;Kajtoch et al., 2016;Médail and Baumel, 2018;Jackson and Cook, 2020). Treating multiple distinct populations as a single entity can result in inflated estimates of genetic diversity and population size estimates, while oversplitting single lineages may have the opposite result. Understanding the phylogeographic patterns within taxa can also further help elucidate how species may have responded to historic climatic fluctuations (e.g., Pleistocene glacial cycles), and thus help to predict how taxa may respond to future climate change (Provan and Bennett, 2008;Row et al., 2011;D'Amen et al., 2013;DiLeo et al., 2013;Papadopoulou and Knowles, 2016;Mascarenhas et al., 2019;Nadeau and Urban, 2019;Luna-Aranguré and Vá zquez-Domínguez, 2020).
North America is an intriguing region for phylogeographic inference because many taxa are broadly distributed across the continent, and recent glaciations during the Pleistocene have likely restricted species to refugia in the southern portions of their current ranges (Swenson and Howard, 2005;Waltari et al., 2007;Provan and Bennett, 2008;Polfus et al., 2017;Lait and Hebert, 2018). Further, multiple geographic features have been shown to be important barriers to gene flow among populations, including the Mississippi and Apalachicola rivers. Indeed, several taxa previously thought to represent single, broadly distributed taxa in this region have, upon more comprehensive (often molecular) investigation, been shown instead to be multiple distinct and geographically segregated species (Burbrink et al., 2000;Burbrink, 2001;Leaché and Reeder, 2002;Lemmon et al., 2007;Crother et al., 2011;Ruane et al., 2014;Barrowclough et al., 2019;Guyer et al., 2020;Jackson and Cook, 2020).
Among these taxa widely distributed across eastern North America is the Rough Greensnake, Opheodrys aestivus, a moderately sized (maximum snout-vent length = 735 mm; Walley and Plummer, 2000) semi-arboreal species with a broad distribution throughout much of eastern North America, ranging from the Pine Barrens of southern New Jersey in the northeast, south to the Florida Keys, and west to northeastern Mexico and the Edwards Plateau of central Texas (Walley and Plummer, 2000;Powell et al., 2016). In the most comprehensive taxonomic review of the species, Grobman (1984) used scalation data from throughout the range and recognized four subspecies: O. a. carinatus from peninsular Florida, O. a. conanti from Virginia barrier islands, O. a. majalis from the western portion of the range (primarily Texas, Oklahoma, western Arkansas, Missouri, and Kansas), and O. a. aestivus from the rest of the distribution. Subsequent authors, however, have recognized only O. a. aestivus and O. a. carinatus, arguing that O. a. conanti, O. a. majalis, and O. a. aestivus represent the arbitrarily divided portions of clinal variation across the range of the species (Frost and Hillis, 1990;Frost et al., 1992;Walley and Plummer, 2000). Whether any of these subspecies represent distinct, evolutionary lineages has yet to be tested.
Here, we use multilocus Sanger-sequence data from samples throughout the distribution of O. aestivus to examine phylogeographic structure and infer the timing of divergence among lineages, while simultaneously determining whether any of the historically recognized subspecies are in fact supported as distinct lineages of Opheodrys. We supplement Sanger-sequence data with genomic data collected via double-digest restriction associated DNA sequencing (ddRADseq) for a subset of samples spanning the distribution and all major clades recovered in the Sanger data set to provide a more genomically comprehensive test of the lineages and putative species boundaries.
(ThermoFisher Scientific, Inc.), and sequenced in both directions using amplification primers by Macrogen, USA (Brooklyn, New York, USA). Complementary sequences were then assembled in Geneious ver. 6.1.8 (https://www.geneious.com; Kearse et al., 2012), and sequences were aligned using Muscle (Edgar, 2004). Alignments were visually inspected to check for errors, and protein coding genes were translated to amino acid sequences in Geneious to check for premature stop codons. We used PHASE ver. 2.1.1 (Stephens et al., 2001;Stephens and Scheet, 2005) to estimate the most probable pair of alleles for nuclear sequences with multiple heterozygous sites. For each locus, PHASE analyses were run for 1,000 iterations, sampling every iteration, following 200 iterations of burn-in. We replicated PHASE analyses for each locus three times to ensure results were consistent, using a different random starting seed for each run.
Sanger-Data Phylogeography.-We used maximum likelihood (ML) and Bayesian inference methods to infer the phylogenetic FIG. 1. Sampling localities for Opheodrys aestivus included in this study, with the distribution of the species based on IUCN Red List data (Hammerson et al., 2007)  GACYTTGTGRACTTCYACRTAATCCAT relationships among populations of O. aestivus. Analyses were conducted on (a) each phased locus individually, with the mitochondrial data partitioned by codon position, (b) the three unphased nuclear loci concatenated and partitioned by locus, and (c) a concatenated analysis of mitochondrial and unphased nuclear data, partitioned by mitochondrial codon position and nuclear locus. We estimated the best fit model of nucleotide substitution for each partition in jModeltest 2.0 (Posada, 2008) using the corrected Akaike information criterion (cAIC). We then used RAxML ver. 8.2.12 (Stamatakis, 2007) (Huelsenbeck and Ronquist, 2001;Ronquist and Huelsenbeck, 2003;Ronquist et al., 2012) using the best-fit models of nucleotide substitution for each partition. Two independent runs of 12 million iterations were conducted, sampling every 1,000 iterations, the first 2,000 samples of which were discarded as burn-in. Parameter trends and effective sample sizes (ESSs) were examined in Tracer ver. 1.6 (Rambaut and Drummond, 2007) to ensure runs had converged at stationarity, and the posterior had been sufficiently sampled. We examined population structure and identified putative lineages within the genus Opheodrys for use in multispecies coalescent analyses in Structure ver. 2.3.2 (Pritchard et al., 2000;Falush et al., 2003). Structure analyses were run assuming the possibility of admixture among populations; 20 replicates were run for each value of K, with each run consisting of a burn-in of 50,000 iterations, followed by 100,000 iterations of sampling with a thinning interval of 10. The optimal value of K was then determined by examining the Ln P(D) plot, and by calculating DK (Evanno et al., 2005) using the web interface of CLUMPAK (Kopelman et al., 2015) and Structure Harvester (Earl and VonHoldt, 2012). Clusters identified in initial analyses were then re-analyzed individually because Structure may detect only higher levels of population structure and fail to detect finer scale populations (Evanno et al., 2005).
Species Tree and Divergence Time Estimation.-We used *BEAST v. 2.6.0 to estimate the multispecies coalescent phylogeny for lineages with strong support as distinct lineages in concatenated analyses, as well as to estimate the timing of divergence among these lineages and between the two species of Opheodrys (Heled and Drummond, 2010;Bouckaert et al., 2014). Mitochondrial data were partitioned by codon position for model of nucleotide substitution and clock model; but, constrained to a single gene genealogy, nuclear loci were unlinked across substitution model, clock model, and gene genealogies. All partitions were set to an uncorrelated log-normal relaxed clock model, and to nucleotide substitution models following the best-fit model as estimated in jModeltest. Analyses were run for 250 million generations, sampling every 10,000 generations, the first 10% of which were discarded as burn-in. Convergence was determined by examining all parameter plots and ESSs in Tracer. Priors on fossil calibration points were placed on nodes in the species tree because the fossil record reflects the history of species rather than that of allelic lineages, and divergences in any gene genealogies or a concatenated phylogeny necessarily predate divergences in the species tree, thus resulting in an overestimate of divergence time (McCormack et al., 2011). Fossil calibrations were given log-normal prior distributions with mean ages and soft 95% prior distribution boundaries following Pyron and Burbrink (2009a) and based on information in Holman (2000) and Head et al. (2016). We used two fossil records to calibrate the phylogeny of Opheodrys: (1) the divergence between Pantherophis guttatus and Pituophis melanoleucus was given a mean of 16 million years (mya) with a prior credibility interval of 11.9-21.0 mya based on the oldest known fossil Pantherophis (Head et al., 2016;Tucker et al., 2014), and (2) the divergence between Cemophora coccinea and Lampropeltis getula was given a mean of 13.75 mya with a prior credibility interval of 8.4-21.3 mya based on the oldest known kingsnake, Lampropeltis similis, from the medial Barstovian age, Miocene epoch.
Demographic Analyses.-To examine the demographic history of O. aestivus, and specifically to test the importance of migration among resulting lineages, we used IMa3 (Hey andNielsen, 2004, 2007;Hey, 2010;Hey et al., 2018) to fit isolation with migration models to the Sanger sequence data. Complexity of the model and difficulties achieving suitable mixing when including all Opheodrys samples dictated that IMa3 analyses were restricted to O. aestivus, with all O. vernalis samples excluded. We tested the fit of five possible demographic scenarios that vary in their history of migration among resulting lineages (see Results below for details on the included lineages: ''Main,'' ''Texas,'' and ''Florida''): (1) no migration, (2) migration restricted to between the Main and Texas lineages (i.e., no migration between Florida and other populations, currently or historically), (3) migration only between sister lineages (i.e., between Main and Texas lineages, and between Florida and the common ancestor of Main + Texas), (4) migration between geographically adjacent populations (i.e., no migration between Texas and Florida lineages), and (5) full isolation with migration model (migration permitted among all pairs of populations).
Following a series of preliminary analyses to identify optimal upper limits on uniform prior distributions for IMa3 analyses, we set the upper limit on divergence time priors as 7.0 and 15.0 for the divergences between the Texas and Main lineages, and Texas + Main and Florida, respectively. Prior limits on population sizes were set to 5.0 for the Texas lineage, and 20.0 for all other current and ancestral populations. Prior limits on migration rates were set to 2.0 between all current population pairs, and 5.0 for migration between the Florida lineage and the common ancestral population of the Texas + Main lineages. Two replicate analyses were run, each consisting of 40 chains with geometric heating parameters 0.98 and 0.65, following the authors' recommendations for small data sets under medium heating to achieve suitable mixing. Each chain was run for 20 million generations, sampling every 100 generations, following a burn-in period of 1 million generations, yielding a posterior sample of 200,000 samples. All posterior samples were then used to calculate joint posterior densities for comparison of nested models using likelihood ratio tests (Hey and Nielsen, 2007).
Genomic ddRAD Sequencing.-To further evaluate putative species boundaries and clades recovered in the Sanger sequencing analyses in a more genomically comprehensive context, we also used a double-digest restriction associated DNA sequencing (ddRADseq) for a small subset of the samples, representing all clades identified from Sanger-sequence analyses. A subsample of 13 samples of O. aestivus, spanning the distribution of the species and all lineages recovered in Sanger-data analyses, along with a single sample the sister species, O. vernalis, were sent to Floragenex, Inc. (Portland, Oregon, USA) for ddRADseq library preparation and sequencing (Appendix), and were included as part of 95 total samples sequenced on a single lane of an Illumina HiSeq 2500 platform. Resultant sequences were demultiplexed, stacks formed, and single nucleotide polymorphisms (SNPs) called in iPyRAD (Eaton and Overcast, 2020). Reads with over five low-quality base calls (Q < 20) and loci recovered from fewer than seven samples (50%) were discarded.
Phylogenomic Analyses.-We used Plink v. 2.0 (Chang et al., 2015) to calculate the number of private SNP alleles for each lineage recovered from Sanger sequence data, as well as to calculate pairwise Fst (Hudson et al., 1992) among these lineages. To visualize population genetic structure based on the SNP data, we used the adegenet v. 2.1.3 (Jombart, 2008;Jombart and Ahmed, 2011) package in Program R to conduct principal components analysis (PCA).
Both maximum likelihood and Bayesian inference were used to estimate the phylogeny of these samples from the ddRADseq dataset. RAxML was used to estimate the ML phylogeny using the full sequences, with 50 search replicates to find the ML phylogeny using the GTR+G model for each partition. Other parameters were left as default values, and 1,000 nonparametric bootstrap replicates were conducted to assess branch support. Bayesian analyses were conducted in Mr.Bayes, and were restricted to include only SNPs because of computational restrictions. Two independent runs of 4 chains were run for 12 million iterations each, sampling every 1,000 generations, the first 2,000 samples of which were discarded as burn-in. Stationarity and convergence were assessed via examining parameter trends and ESSs in Tracer.
We used SNAPP (SNP and AFLP Package for Phylogenetic analysis; Bryant et al., 2012) to estimate the species tree for Opheodrys using the genomic data, assuming the three clades found within O. aestivus (see results) as distinct lineages. Analyses were run for 6 million generations, sampling every 1,000 generations, the first 10% of which were discarded as burn-in. Parameter trends and ESSs were examined in Tracer to ensure chains had reached stationarity and sufficiently sampled the posterior.
Concatenated phylogenetic analyses ( Fig. 2A)  Multispecies coalescent analyses in *BEAST yielded similar results to the concatenated phylogenetic analyses, with strong support (PP = 1.00) for the monophyly of the genus Opheodrys, and for the monophyly of O. aestivus s. l. (BPP = 0.93). However, relationships among the three major clades within O. aestivus s. l., were unresolved (Fig. 3). Credibility intervals in divergence time estimates were relatively large (Fig. 3), likely as a result of the limited number of fossil calibrations and the lack of an informative fossil calibration within Opheodrys. However, we recovered an estimated divergence between O. aestivus s. l. and FIG. 3. Calibrated species tree for Opheodrys aestivus from Bayesian multispecies coalescent analyses of Sanger data in *BEAST. Fossil calibration points are indicated with asterisks. Values above branches indicate mean estimated node age in millions of years before present (mya), bars on nodes correspond to 95% credibility intervals. Values below branches correspond to BPP branch support. Colors of symbols at branch tips correspond to major clades in Figures 1 and 2. O. vernalis of 6.43 (3.36 -9.83) mya, reflecting a late Miocene divergence. Within O. aestivus s. l., we recovered an estimated divergence of 1.34 (0.32-2.82) mya between the Main and Texas clades, reflecting a divergence in the early Pleistocene or near the Pliocene-Pleistocene border, and of 2.94 (0.61-6.09) mya between these clades and the Florida clade (Fig. 3), reflecting an earlier divergence in the Pliocene or early Pleistocene.
Coalescent analysis of isolation with migration models in IMa3 found the best fit demographic model to include migration between geographically adjacent taxa (i.e., between Texas and Main clades, between Main and Florida clades, and between Florida and Texas+Main clades) or no migration among lineages (Table 2). In the adjacent migration model that includes migration between the Florida and Main clades and in the full model including all migration parameters, the migration rate between the Florida clade and other clades is much lower than that between the Texas and Main clades (Table 2).
We recovered 20,691 loci (9,970-19,382/sample) via ddRADseq that passed filters, including 81,258 SNPs. Of these, 33,983 SNPs were recovered from all lineages identified from the Sanger sequence data, including 4,159 SNPs with private alleles restricted to the Florida clade, 10,583 SNPs with private alleles restricted to the Main clade, 2,437 SNPs with private alleles restricted to the Texas clade, and 9,268 SNPs with private alleles restricted to O. vernalis (Table 3). The PCA of the ddRADseq data supported substantial differentiation among the Sanger sequence-based clades. The first principal component (28.9% of variation) primarily showed differentiation between O. aestivus and O. vernalis, while the second and third principal components (12.4% and 10.5% of variation, respectively) revealed strong differentiation of the Florida and Texas clades from each other and from the O. aestivus Main clade and O. vernalis (Fig. 4).
The data set included 18,582 SNPs that were independent (i.e., from different loci) and biallelic, and thus suitable for, and used in, SNAPP species-tree analyses. Raw ddRADseq reads and alignments of filtered loci are available via Dryad (Repository number: 10.5061/dryad.brv15dv8w).
Concatenated phylogenomic analyses of ddRADseq data strongly corroborate results of Sanger data-based analyses (Fig. 5). The monophyly of O. aestivus was strongly supported (MLBS = 100, BPP = 1.00), as was that of the Main clade (MLBS = 100, BPP = 1.00), and for the sister relationship between the Main and Texas clades (MLBS = 100, BPP = 1.00). Species-tree analyses in SNAPP corroborated other analyses with respect to high support for a sister relationship between the Main and Texas clades (BPP = 1.0; Fig. 5). However, as in the species-tree analyses based on Sanger data, the relationships between the Main + Texas clades, the Florida clade, and O. vernalis were not well-supported (BPP = 0.62).

DISCUSSION
The divergence between the peninsular Florida populations and other populations is not surprising, given the morphological differences between these clades and the long history of documented phylogeographic breaks between peninsular Florida and the rest of continental North America in a wide variety of taxa. Similar breaks were documented in some of the earliest phylogeographic studies (Avise et al., 1987;Avise, 2000), and have been reported in a wide variety of terrestrial taxa (Avise et al., 1987;Swenson and Howard, 2005;Seal et al., 2015;Fetter and Weakley, 2019), including amphibians (Austin et al., 2002;Means et al., 2017;Barrow et al., 2018), turtles (Walker et al., 1998a,b), squamates Fontanella et al., 2008;Manthey et al., 2016;McKelvy and Burbrink, 2017), birds (Barrowclough et al., 2018(Barrowclough et al., , 2019, and mammals (Avise et al., 1987;Ellsworth et al., 1994;Jackson and Cook, 2020 Burbrink et al., 2008), or potentially much younger divergences (albeit without published date estimates, e.g., in turtles [Walker et al., 1998b] and birds [Barrowclough et al., 2019]). Combined, these shared geographic patterns of divergence make it clear that the history of peninsular Florida, either as a single event or as a result of multiple waves of divergence, is an important factor that has played a major role in the evolutionary and biogeographic history of much of North American biota. The lack of structure within the Main O. aestivus clade, however, is less expected, because we found little to no divergence among populations from the Texas Gulf Coast, east to the Florida panhandle, and north to southern New Jersey (Figs. 1, 2). Most squamates show substantial population structure in this region, at least across major barriers or ecological transitions such as the Mississippi River and/or as habitat changes from eastern deciduous forest to grasslands (e.g., Coluber constrictor ; Lampropeltis getula complex [Pyron and Burbrink, 2009b;Krysko et al., 2017]; Lampropeltis triangulum complex [Ruane et al., 2014]; Pantherophis guttatus complex [Burbrink, 2002;Myers et al., 2020]; Pantherophis obsoletus complex [Burbrink et al., 2000[Burbrink et al., , 2021Burbrink, 2001]; Sceloporus undulatus complex [Leaché and Reeder, 2002]). Many squamates also show more fine-scale phylogeographic structure within this region (e.g., Diadophis punctatus ; Scincella lateralis [Jackson and Austin, 2012]). However, our observed pattern of phylogeographic structure is not unique to O. aestivus. The Scarletsnake, Cemophora coccinea, similarly shows limited population structure across the region from western Louisiana to southern New Jersey (Weinell and Austin, 2017). Thus, our data suggest that while many species may have persisted in, and expanded from, multiple Pleistocene refugia, others, including C. coccinea and O. aestivus, may have recolonized northern regions primarily from single refugia in the postglacial period. Additional research should focus on the ecological similarities between those taxa that have apparently expanded from a single refugium following Pleistocene glacial retreat versus those that apparently expanded from multiple refugia, or test whether these differences are merely a result of stochasticity.
We found evidence for an additional clade, perhaps even a distinct species, that appears restricted to the Edwards Plateau region of Texas and may have become isolated during Pleistocene glacial cycles, based on our estimated divergence time of 1.99 (0.30-2.62) mya. The Edwards Plateau is a unique ecosystem that occurs at the confluence of the Gulf Coast, Great Plains, and Chihuahuan Desert biomes, and represents the northeastern-most portion of the distributions of several primarily Mexican species (e.g., Thamnophis cyrtopsis, Gerrhonotus infernalis; Webb, 1980;Powell et al., 2016). Numerous other taxa with broad distributions across eastern North America similarly have clades restricted to this region, including Pantherophis bairdi in the P. obsoletus complex (Burbrink et al., 2000;Burbrink, 2001), Agkistrodon laticinctus in the A. contortrix complex (Burbrink and Guiher, 2015), and some clades of Scincella lateralis (Jackson and Austin, 2010), among others (Baird et al., 2006;Rodriguez et al., 2012;Thomson et al., 2018). Although estimates of the divergence times of these other lineages restricted to the Edwards Plateau are limited, these data support the importance of this region as a possible Pleistocene refugium for numerous disparate taxa, and its importance in driving the generation of diversity in taxa currently broadly distributed across eastern North America.
Taxonomic Treatment.-Our analyses recover a deep and strongly supported divergence between the populations of O. aestivus s. l. from peninsular Florida and all other O. aestivus populations, with our estimates suggesting this lineage diverged from the rest of O. aestivus s. l. nearly 3 mya, during the late Pliocene (Figs. 2, 3). Further, demographic analyses of isolation with migration models find that migration between this peninsular Florida population and other populations of O. aestivus is either excluded from the model or very low (Table  2). This clade is represented by a limited sample size (Florida clade n = 3), but these samples span the northern and southern extremes of the Florida peninsula (Fig. 1), and the northern-most sample in this clade is <100 km from samples in the Main clade (mean cytb uncorrected p-distance = 0.099). This is further supported via previously published morphological studies with geographic sampling that appears to correspond with distinctness of this lineage. Grobman (1984) recognized the subspecies O. a. aestivus, with a type locality of 1 mile west of Parksville, McCormick Co., South Carolina (well within the distribution of the Main clade we recovered), as including the majority of the eastern North American populations of O. aestivus, but found that peninsular Florida populations had more strongly keeled scales, and that the keels extended laterally to the third dorsal scale row, although it is worth noting that squamate scales are known to vary within species across their geographic ranges (e.g., Natrix tessellata; Mebert, 2011). Based on his results, Grobman (1984) described these peninsular Florida populations as the subspecies O. a. carinatus, the type locality of which is Archibold Biological Station, Highlands Co., Florida. Further, Plummer (1987) found the peninsular Florida populations to be larger in mean and maximum size relative to other populations of O. aestivus. Combining evidence from our results and these previously published studies, we suggest that the peninsular Florida populations be elevated to species status, encompassing the former subspecies O. a. carinatus and to be recognized as Opheodrys carinatus Grobman 1984 stat. nov. We also find the central Texas populations from the Edwards Plateau as distinct and diverging from the Main clade of O. aestivus near the Pliocene-Pleistocene boundary (Figs. 2, 3). It is possible this is the same entity as the previously described subspecies O. a. majalis, in part, as recognized by Grobman (1984). However, the type locality of majalis is in Indianola, Texas, on the Texas Gulf Coast, a different ecoregion to the Edwards Plateau, to which our data suggest this clade may be restricted, and our other Gulf Coast sampling from further east are within the Main O. aestivus clade (Figs. 1, 2). Further, previous morphological studies have found far less distinctive differences between these populations and the rest of O. aestivus, relying largely on broadly overlapping differences in the number of ventral and subcaudal scales (Grobman, 1984). As such, we refrain from elevating the central Texas clade to any specific status pending denser molecular sampling that can better refine the distribution of this clade, determine its status as a distinct species, and elucidate whether this clade perhaps represents majalis as distinct or a wholly undescribed lineage. At this time, we recommend that no subspecies be recognized for O. aestivus because we found no evidence that they (O. a. majalis, O. a. conanti) represent any clearly recognizable evolutionary entities, and suggest that O. vernalis, O. aestivus, and O. carinatus best characterize the diversity found within Opheodrys based on currently available information.
Conclusions.-Our data show that, despite having previously been considered a single broadly distributed species, O. aestivus is likely a complex of at least two species: a peninsular Florida species; a broadly distributed eastern North American species; and, potentially, a third species in central Texas. Our results display an unusual combination of patterns found in other taxa. The divergence of peninsular Florida has long been a documented pattern, and several taxa or lineages are restricted to central Texas, yet few other squamates show such low levels of divergence from east Texas to the Florida panhandle and north to New Jersey.
Our approach also illustrates how dense geographic sampling can be used as a guide for selecting samples for more comprehensive genomic sequencing as a cost-effective means of combining intense geographic and genomic sampling. Finally, we recover identical phylogenies and similar support values and relative branch lengths based on Sanger and ddRADseq data, suggesting that while there are certainly situations where genomic data is necessary, particularly in taxa with short internodes or introgression, a small number of loci may, in some cases, be used to resolve phylogeographic or phylogenetic questions as effectively as genomic sampling.
APPENDIX.-Collection localities and GenBank accession numbers for samples included in this study. Samples in bold indicate inclusion in both Sanger and ddRADseq datasets. Str. Label indicates code used to label structure plots, with samples organized from northeast to southwest within each clade. Collection codes abbreviated as follows: CAS: California Academy of Science, FTB: Frank T. Burbrink