Mechanisms of speciation in Coluber constrictor
Abstract
Areas with high species endemism and unique habitats, such as Florida, can provide insight into mechanisms of speciation. Related, delimiting species is integral to understanding evolutionary processes of speciation and establishing the taxonomic status of metapopulations. Genetic analysis of Florida’s taxa has revealed several cryptic species. Identifying cryptic species has often used a single locus. Relying on a single locus for the evolutionary history of a species has been shown to give misleading results due to incomplete lineage sorting and introgression. These single locus studies require revisiting within multi locus coalescent frameworks. Here, we follow up previous research on a single locus in Coluber constrictor with multi locus coalescent analyses and ecological niche modeling. We find support for the proposal from previous studies that the population of this snake in peninsular Florida should definitively be raised to species status.
Introduction
Understanding diversity and areas where prolific speciation occurs is integral to conservation practices and prioritizing areas in need of protection, as well as cataloging the Earth’s biodiversity (Dirzo & Raven, 2003). However, awareness of such diversity can be thwarted by the existence of cryptic species (Bickford et al, 2006). Therefore, molecular species delimitation is important in uncovering underlying diversity.
Florida presents a unique opportunity for the study of the process of speciation owing to its subtropical latitudes, peninsular form, and complex geological history (Lodge, 2005), which could all be major drivers of speciation. It is considered a hotspot of endemism within the US and importantly, with no direct connection to the tropics (Reese, 2013). Temperatures in south Florida rarely fall below freezing; this is important in terms of possible lack of defined breeding seasons, as well contributing to ecological systems that continuous growth creates. Other geographical areas that have similar climate, such as parts of the Caribbean, are blocked from Florida by large bodies of water to the east, south, and west. The northern part of the state becomes a barrier to species intolerant of temperate climates, though this mostly holds true for plants, as most of the vertebrates that live there have dispersed into Florida from temperate North America (Lodge, 2005), creating a mix of temperate and subtropical species in the ecosystems. Although current sea levels act as a barrier on three sides of Florida, the geological history of Florida shows periodically higher sea levels during interglacial epochs, transforming the peninsula into a chain of islands (Lodge, 2005) and isolating insular populations from the rest of the continent via the Suwannee Channel at roughly the Florida/Georgia border (Webb, 1990). These conditions may have generated species with unique morphological characteristics found only in ecoregions specific to South Georgia and Florida, such as the Southeastern Conifer Forests and The Florida Sand Pine Scrub (Ricketts & Imhoff, 2003). Found specifically in Florida are the Florida Sand Skink (Plestiodon reynoldsi), the Florida Scrub Jay (Aphelocoma coerulescens), and the Scrub Palmetto (Sabal etonia). Additionally, there are several species with unique phylogeographic lineages found only in the Florida peninsula, such as the White-Tailed Deer (Odocoileus virginianus) (Ellworth et al., 1994), the Oldfield Mouse (Peromyscus polionotus) (Avise et al., 1984), and several turtle species (Walker and Avise, 1998), yet it is unclear if these represent unique taxa given that unique morphological differences are lacking between the Florida lineages and elsewhere.
One of the most deeply divergent lineages of any species found in the Southeastern United States is within the North American racer (Coluber constrictor), a snake representing a unique opportunity for understanding mechanisms of speciation in Florida (Burbrink et al. 2008). Additionally, it is one of the few vertebrates found across the entire continent of North America, ranging longitudinally from Canada to Guatemala (Wilson, 1978). A previous phylogeographic study using only mtDNA identified six lineages across North America, including a deeply divergent lineage found in peninsular Florida (Burbrink et al., 2008). The date of origin of this lineage from the other mtDNA clades was estimated to be in the late Miocene/Early Pliocene, about 6 MYA. This study, and several like it, including the aforementioned research on O. virginianus and P. polionotus, relied on a single mitochondrial DNA (mtDNA) locus to infer phylogeographic history of these species. Unfortunately, it may be that these studies only reflect the history of the locus and not that of the species as a whole (Carstens and Knowles, 2007). However, Myers et al (in review) investigated the genus Coluber with regards to resolving taxonomic nomenclature between Coluber and Masticophis, its sister taxon, using multi-locus coalescent analyses, and found deep divergence in the species tree between peninsular Florida and eastern North America populations. Significant ecological niche differences have also been shown in another snake in Florida, in particular the populations of the Florida Cottonmouth (Agkistrodon conanti) versus the Northern Cottonmouth (Agkistrodon piscivorus; Burbrink & Guiher, 2015).
While it seems evident that there is significant divergence in the Florida from the continental populations, what remains unclear is what drove the divergence. One possibility is that elevated seas levels isolated taxa and subsequent ecological speciation occurred. Another possibility is that the unique ecosystems and subsequent ecological niches created within Florida isolated species and reduced gene flow.
Delimiting species found in closely related taxa has been challenging to evolutionary biologists, (Carstens et al, 2013). In the last few decades, the increasing ease in obtaining genetic data has greatly expanded the ability to study speciation and identify cryptic taxa. Coalescent-based approaches provide a testable method to not only delimit taxa, but understanding what is driving population differentiation (Yang & Rannala, 2010). These methods provide a more robust framework for identifying unique taxa and have been used in many recent studies for delimitation and the inference of species trees (e.g., Myers et al 2013; Ruane et al 2013). When sampling multiple loci, coalescent-based methods can account for incongruences between gene and species histories, such as incomplete lineage sorting and gene tree heterogeneity (Dupuis 2012). In addition, parameters such as effective population size, and divergence times can also be simultaneously estimated (Yang & Rannala, 2010). Best practices suggest that coalescent-based methods be used when there is the ability to sample many individuals across the range of the species being studied, and when the loci are selection-free (Rosenberg and Nordberg, 2002).
As a new area of research, species delimitation via genetic analysis has had its difficulties in development during the maturation process of the method given vulnerabilities in techniques that only become apparent through trial and error. Being aware of these as they arise and using a combination of programs to offset the potential issues can give robust results and avoid false positives (Carstens et al., 2013). This is why a coalescent-based analysis like BPP with a population assignment program like Structurama. Ultimately, a multi-analysis approach to species delimitation is best; integrative taxonomy should be used instead of reliance on any one particular method (Padial et al, 2010).
Here, I examine lineage divergence of C. constrictor between peninsular Florida and the rest of continental North America using coalescent-based methods to address four questions. First, given potential problems with single locus studies, we determine if mtDNA provided a correct estimation of phylogeographic history. Second, is there support for a distinct species of racer in peninsular Florida? Third, is there evidence of gene flow between populations in peninsular Florida and the lineage on the continental US? Fourth, can any divergence be attributed to differences in ecological niches? We find our results show significant genetic divergence between these two populations, and add to the growing body of data that show Florida to be a distinct environment from the continental United States, as well as a biodiversity hotspot for generating species diversity.
Materials and Methods
Taxon Sampling and Sequence Data
I obtained 112 tissue samples of C. constrictor ranging from south Florida to Virginia, and west to southern Louisiana (Fig. 1, Appendix 1). DNA was extracted using Qiagen DNeasy kits from shed skin, liver, or muscle tissue. One mtDNA locus (cytb) and three nuclear DNA (nDNA; cmos, r35, nt3) loci, all protein coding genes, were amplified with GoTaq Green Master Mix (Promega©) according to manufacturer protocol. PCR products were cleaned using ExoSap-IT (USB©). The sequencing reaction consisted of 3ul DTCS (Beckman Goulter), 2 ul primer, 2 ul template and 3 ul deionized water. For PCR amplification of cytb we used the H14910 and THRSN2 primers (Burbrink et al., 2008), and H19410 for cycle sequence reactions, in one direction. Primers for cmos were S77 and S78 found in Lawson et al. (2005), for r35 were GeneF and GeneR (Vidal & Hedges, 2005), and for nt3 were F3 and R4 (Chippindale & Noonan, 2006). No gaps were present in any of the genes among the samples, and sequences were aligned in Geneious 6.1. using Geneious Read Mapper (Kearse et al., 2012). For the three nDNA, phasing of heterozygous genotypes was conducted in Phase v2.1.1 (Stephens and Donnelly, 2003), with parameters set to run for 100 iterations with a thinning interval of one, and a burn-in of 100. Phased genotypes were used for all coalescent analyses. We used DnaSP v5 (Librado & Rozas, 2009) to calculate polymorphic sites, haplotypes, and nucleotide diversity for each locus.
Gene-Tree Analysis
Phylogeographic structure for each gene was estimated using maximum likelihood in RaxmlGUI 1.1 (Silvestro and Michalak, 2011). All trees were rooted with Coluber flagellum, the sister taxon to the species C. constrictor (Myers et al, in press). Nucleotide substitution was modeled with GTRGAMMA, and 1000 nonparametric bootstrap replicates were run to assess nodal support. Clades with values > 75% were considered well supported (Felsenstein, 2004).
Speciation and Divergence
In the course of this study, I used the general lineage species concept as an evolutionary framework to delineate taxa, where species are defined as separately evolving metapopulations lineages (De Quieroz, 2007). To assess the number of populations of C. constrictor in the Southeastern US, we used Structurama 2.0 (Huelsenbeck and Andolfatto, 2007), which placed individuals into probable populations using a Dirichlet process prior for K populations with a gamma shape set at 0.1 and a scale parameter at 10.0 (Pella & Masuda, 2006). I ran Structurama four times using MCMC for 5 x 105 generations, saving every 100th sample, with a burn in of 10 x 105. The number of K populations was chosen based on the one with the highest posterior probability value.
To determine if the populations recognized by Structurama represent distinct species, Phylogenetics and Phylogeography was used (BPPv2.0; Yang and Rannala, 2010). This method employed the species phylogeny and assesses lineage sorting due to ancestral polymorphism. A guide tree is required to assesses probabilities of nodes and delimit species in BPP, however since only the existence of only a single node was examined, the guide tree is elementary. Fine-tuning parameters were set to the recommended values of swapping rates ranging between 0.30 and 0.70. All four loci and three combinations of priors for ancestral populations size (θ) and root age (τ0) to investigate possible timing of divergence and population size changes; these were set to θ~G (2, 2000) and τ0~G(2, 2000) for small ancestral population and shallow divergence, and for large ancestral population and shallow divergence, θ~G (1, 10) and τ0~G(2, 2000), and finally, for large ancestral populations and deep divergence, θ~G (1, 10) and τ0~G (1, 10). We ran BPP with these parameters three times each for 1 x 106 generations, with a burn-in of 1 x 104, and a sampling frequency of once every five generations.
Incomplete Lineage Sorting
Initial analysis showed the mitochondrial gene trees and the nuclear gene trees were not generally reciprocally monophyletic with respect to the population sorting assessed by Structurama. I therefore assessed the degree of discordance, which could be due to factors such as incomplete lineage sorting or introgression by calculating the genealogical sorting index (gsi) (Cummings et al 2008). The results of this index is based on a scale of 0 – 1, with 0 indicating complete non-exclusivity, and 1 indicating monophyly of the gene. Maximum Likelihood gene trees were used for the analyses and the index was estimated on the gsi web server (http://www.genealogicalsorting.org). P values were estimated based on 10 000 permutations.
Processes of Diversification:Estimating Divergence Dates and Migration
To estimate divergence time, ancestral population sizes, and migration we used the program IMa (Hey and Nielson, 2007) on the mtDNA locus and the fully phased data set of the nDNA sequences. The posterior probability estimates for six factors of interest for studies of population divergence were calculated: Ne1 (effective population size of population one-peninsular Florida), Ne2 (effective population size of population two-the continental populations adjacent to peninsular Florida), Nea (ancestral effective population size), m1 (migration of Floridian population into continental populations), m2 (migration of continental populations into Florida), and τ (time of divergence). I used a mtDNA (cytB) substitution rate of 1.086 x 10-5 substitutions/locus/year (95% HPD: 7.8 x 10-6 – 1.4 x 10-5) (Burbrink et al., 2011) to adjust units to years and effective population size and a generation time of two years (Fitch, 1963) with a substitution model of HKY (Hasegawa et al., 1985) applied to all genes. Inheritance scalars were set at 0.25 for cytb and at 1.0 for the three nDNA loci. To estimate priors for the final analyses, the data set was run in M mode with different seeds and geometric heating of 30 chains. The analysis was run three times on the final data set at 1 x 106 trees with the following priors: q1 = 5, m1, m2 = 2, and τ = 10. This output was used to run log-likelihood ratio tests (2LLR) to examine the significance of a series on nested models in the “L” mode of IMa. All nested models are tested against the five parameters of the full model, such that Ne1 ¹ Ne2 ¹ NeA, and m1 ¹ m2. Therefore, the 16 nested models test if some or all of these parameters are equivalent to the full model using the 2LLR test from a chi-square distribution. Also, m1 and m2 are evaluated where effective migration = 0.
I also used MIGRATE-N, which implements the island model (Beerli 2006), to test for gene flow between populations. Four models of migration were tested. The first assumes the 2 populations have equal migration rates in both directions; the second assumes there is only gene flow from the Florida population to the continental lineage; the third assumes gene flow from continental lineage to Florida; and the fourth assumes the two populations are only a single breeding population. Model fit was assessed by using Bayes Factor (Beerli & Palczewski, 2010). Additionally, I compare the results of gene flow from these two programs to the theory that the limit for gene flow in light of this species assumption is about 1 individual per 10 generations, or 1 gene copy per generation (Zhang et al, 2011).
Ecological Niche Modeling
I tested if these two lineages occupy distinctly different niches, given that significant differences in ecological niches is further evidence of speciation. To construct Ecological Niche Models (ENMs) I used 19 bioclim variables describing variation in precipitation and temperature from the Worldclim data (Hijmans et al., 2005) taken at 30-s resolution. I used ENMTools 1.3 (Warren et al., 2010) to test the hypothesis that the predicted niches of the two populations were more different than expected by chance (Warren et al., 2008). ENMs were reconstructed for both populations using MAXENT v3.3.3 with default parameter settings (Phillips, et al, 2006). Model performance was assessed by evaluating the receiver operating characteristic curve (ROC) and the area under the ROC curve (AUC) statistic. The niches of the two lineages are identical if the observed ENMs are no more different than randomly drawn pairs of samples. I also tested the hypothesis that the population split was due to an environmental barrier by using the range-break blob test. Using the same data, this test randomly draws a line through the geographic area of both populations, calculates the ENMs on either side of the simulated line, then estimates the I statistic and Schoener’s D(Warren et al., 2010) for niche divergence. If the observed statistics are not within 95% of the simulated statistics, then the geographic areas being tested can be considered statistically significant (Glor & Warren, 2010). Testing both of these hypotheses was done by comparing the observed niche overlap to a null distribution of 100 replicates simulated in ENMTools. Statistical significance was determined by using one sample t-test in R, and using the I statistic (Warren et al., 2010).
Results
Sequence Data
A total of 1886 base pairs were sequenced for analysis. No gaps were detected for cmos (531 bp), cytb (415 bp), nt3 (467 bp), or r35 (473 bp). All112 samples were amplified and sequenced for cmos, 108 for cytb, 109 for nt3, and 111 for r35. No individual was missing more than one gene sequence in any analysis. DnaSP analysis showed lower allelic diversity in the Florida population for three out of four of the genes, and greater diversity in the continental population with except for R35 (Table 1). All genetic sequences have been deposited in GenBank (Appendix 3).
Speciation and Divergence
All four runs of Structurama divided the individuals into the two groups of Florida and the continental populations (Fig 1) with an 86-90% probability. The remaining percentages identified a third population with 10-13% probability. There was no support for one or four populations.
Every run of BPP delimited the two populations, Florida and the continental, into separate species with a probability of 1 under all three models of divergence. Burnin occurred by 10,000 generations, and a sampling frequency of one out of every five generations used in the posterior probability distribution.
Incomplete Lineage Sorting
Gene trees inferred with maximum likelihood in RaxML (Figures 2-5).
The gsi analyses shows that all loci were significant with respect to sorting for Florida lineage, where only cytb and r35 are significant for the continental lineage (Table 2). In Florida, Cytb approaches monophyly (0.9263), followed by cmos (0.7364) and r35 (0.4246), with nt3 being the least monophyletic (0.2472).
Migration
IMa estimated migration, divergence times, and ancestral population size using estimates from 99,000 samples after burn-in and attained ESS values for all parameters > 136 (most were above 32,000), with burn in occurring at 100,000 generations. The mean time across all runs of the most recent common ancestor were as follows: cmos – 7.907 Ma, cytb – 7.656 Ma, nt3 – 4.181 Ma, and r35 – 5.846 Ma. The estimate of divergence is consistent across runs (1.3 Ma, 95% CI = 0.96-1.8 Ma) dating it to the Lower Pleistocene. The model with the lowest AIC showed different ancestral and extant population sizes with some migration between the two taxa (m1 = 0.25, 95% CI = 0.06-0.51; m2 = 0.26, 95% CI = 0.11-0.46). The DAIC was 1.94. The predicted effective population size of the ancestral population (NeA; 8,297,500; 95% CI = 5,154,300-12,291,000) was significantly smaller than the Florida population (NeF; 15,931,500; 95% CI = 11,384,800-21,655,600) with the continental lineage being the largest (NeC; 23,483,300; 95% CI = 18,040,100-30,217,800). It also showed equal migration between the extant populations (P < 0.005; 2LLR = -0.0006-0.0125; Appendix 4) with 1.84 gene copies per generation going into Florida from the continental population, and 2.87 in the reverse direction. MIGRATE-N 3.3.2 (Beerli, 2006) showed support for the full model, in which 1.3 gene copies per generation went into Florida from the continent, and 1.7 in the reverse direction, mirroring the results from IMa.
Ecological Niche Modeling
The results from I and Schoener’s D statistics both significant and thus I report only the statistic I values to avoid redundancy, and because it tends to have the greater variation of the two (Glor & Warren, 2010). The Identity Test showed that the ENMs are not identical, with significant P value (t = 322.7416, df = 99, P value < 2.2e-16, 95 percent confidence interval: 0.9509631, 0.9577796). This was calculated using the observed summary statistic for I (.40). The range break blob test supports these results (t = 15.5827, df = 99, P value < 2.2e-16, 95 percent confidence interval: 0.5634646, 0.6111684), using the same observed I statistic. This suggests the observed climates have greater divergence between Florida and the continental U.S. than would be expected by randomized geographic boundaries.
Discussion
Historically, mtDNA has been the primary locus for all phylogeogeographic studies and often still provides the basic structure for subsequent studies. Few have followed up these primary studies within a multiloci coalescent framework to retest previous outcomes (Myers et al 2013). I found here that the mtDNA analysis by Burbrink et al (2008) holds up well under multi-locus analyses, although for snakes, this is not always the case (Ruane et al, 2014; Myers et al, 2013).
The multilocus coalescent analyses here indicates that the Florida clade of C. constrictor is distinct from the continental lineages in North America, with the population range existing approximately at the boundary of Georgia and Florida. During the Pleistocene (17,000 to 2.8 MYA), Florida experienced fluctuating sea levels through 4 glaciation cycles (Lodge, 2005). Prior phylogeographic studies of peninsular Florida have shown divergence between populations driven by vicariant events related to glacial cycling (Soltis et al., 2006; Hewitt, 2000). Divergence time estimates for racers in this region are consistent with a Pleistocene diversification hypothesis (95% CI: 0.96-1.8 MYA).
Three out of four of the genes sequenced show lower allelic diversity in the Florida population. A bottleneck event could explain this finding, though more loci and additional methods would be necessary to confirm this. Discordance between the gene trees and species trees are most likely due to incomplete lineage sorting (Degnan and Rosenberg, 2009), though our results indicate some low levels of migration, which could also account for some discordance. This may be able to be resolved with more individuals, longer DNA sequences, or more genes. The difference in migration results between Migrate and IMa may be explained by the underlying assumptions in each program;. Migrate assumes allelic coalescence between populations is due to migration and population size whereas IMa assumes this is because of migration and the history of the populations themselves, which could account for the discrepancy (O’Meara et al., in press). The second best fit model in IMa considers that the migrations rates between populations are distinct (Appendix 4), which is similar to the Migrate results. In both cases it is clear that diversification has not occurred without geneflow, thus supporting a parapatric model of speciatipon.
The ecological speciation model states three causes of divergent selection: differences in environment (habitat structure, climate, resources, predation and competition), sexual selection, and intraspecies interaction between populations (Rundle and Nosil, 2005). There is evidence that at least the first is acting on the peninsular Florida population of C. constrictor. I found there to be significant differences in the ecological niches of the continental and Florida populations. This is very likely due to the differences between the temperate climate of the continent and the sub-tropical climate of Florida. Another factor to consider is differing morphology. Body size in Peninsular Florida for C. constrictor was found to be significantly smaller on average than the rest of Eastern North America (Auffenberg, 1955; Steen et al, 2013). This may be due to resource partitioning with M. flagellum (Steen et al, 2013) in the differing environment of Florida from the rest of the continent, or is merely corollary, considering this effect was not found in other areas of the country where the two snakes are sympatric. Auffenberg (1955) supported various subspecies of C. constrictor based on various morphological characteristics, such as coloration, number of scales, and hemipenal characteristics. Since C. constrictor is a generalist, only avoiding deserts and high elevations and latitudes, these differences are more likely due to response in variations within the environment in which these populations are located. Auffenberg (1955) shows that hemipenal basal spines are longer in C. c. priapus (hence the name), but that this taxon was found from Florida, north to North Carolina, and west to Louisiana and C. c. paludicolus, which is found in two disparate populations in the Everglades and Cape Canaveral, Florida. Other than body size, these morphological differences are at odds with the coalescent and spatial studies presented here.
Speciation in racers in the SE US is likely a multi-part process: First, populations were likely isolated by a vicariant event associated with rising sea levels followed by reinforcement due to ecological factors divergence to complete speciation (Nosil et al, 2008). I point out that it is possible that these two distinct racers are in the processes of despeciation, where resumed gene flow breaks down genetic and ecological differences between sister taxa. Research into despeciation is uncommon, but bears further study in racers using next generation sequencing data and better modeling approaches.
Taxonomic conclusions
Using coalescent methods for species delimitation, I support two cryptic lineages of racer, one in peninsular Florida and one in the continental United States. These populations of C. constrictor are likely each unique species supported by coalescent delimitation techniques, differences in ecological niche, low gene flow between delineated populations, and body size. I present the nomenclature for the new Florida lineage and the continental taxon, a species synonymy, and description of the range of the taxon.
Coluber Priapus (Dunn & Wood, 1939)
Common Name: Florida Racer
Holotype: ANSP 16111 (priapus)
Type Locality: West Palm Beach, FL.
Etymology: Named for the larger hemipenal basal spines compared to the rest of C. constrictor populations.
Synonymy: This species includes previously recognized subspecies of Coluber constrictor priapus (Dunn & Wood, 1939), Coluber constrictor helvigularis (Auffenberg, 1955), and Coluber constrictor paludicola (Auffenberg & Babitt, 1955).
Distribution. Coluber priapus is comprised of all populations found in peninsular Florida, the northern boundary of this species found in the panhandle to the west including the Georgia border to the east, and south into the upper and lower Florida Keys.
Coluber Constrictor (Linnaeus, 1758)
Common Name: North American Racer
Holotype: NHR Lin-456 (formerly MDG), a 1380 mm female (P. Kalm, Sept. 1748-Feb. 1751, via Mus. Drottn.).
Type Locality: “Canada,” (possibly in error). Probably vicinity of Philadelphia, USA fide Dunn & Wood (1939:1).
Etymology: A misnomer; C. constrictor does not constrict its prey before consuming.
Synonymy: This species includes subspecies Coluber constrictor anthicus (Cope, 1862), Coluber constrictor etheridgei (Wilson, 1970), Coluber constrictor flaviventris (Say, 1823), Coluber constrictor foxii (Baird & Girard, 1853), Coluber constrictor latrunculus (Wilson, 1970), Coluber constrictor mormon (Baird & Girard, 1852), and Coluber constrictor oaxaca (Jan, 1863). Burbrink et al 2008 found that these populations genetically sorted into groups via mtDNA analysis that geographically correlated with the following names according to Joe Collins; They also need nDNA analysis to determine if they should be raised to species status; Coluber constrictor (Eastern lineage), Coluber flaviventris (Central lineage), Coluber mormon (Western lineage), and Coluber latrunculus (Gulf Coast lineage).
Distribution: Coluber constrictor is comprised of all populations found across continental North America, south into Mexico and Guatamala, and north into Canada, excluding Florida and southern Georgia. Not found in deserts or or higher mountain elevations.
Acknowledgements
I would like to thank X. Chen, E. Myers, and S. Ruane for help with analyses and edits to this manuscript, A. McKelvy for help collecting specimens, Kenneth Krysko and the Florida Museum of Natural History for their tissue donations. I would also like to thank my thesis committee members Frank Burbrink, Eugenia Naro-Maciel, and Richard Veit. This research was supported, in part, by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grants CNS-0855217 and CNS-0958379.
References
Avise, J.C., Shapira, J.F., Daniel, S.W., Aquadro, C.F., Lansman, R.A., 1983. Mitochondrial DNA differentiation during the speciation process in Peromyscus. Mol. Bio. Evol. 1, 38-56.
Auffenberg, W., 1955. A reconsideration of the racer, Coluber constrictor, in eastern United States. Tulane Studies in Zoology. 2, 87-158.
Beerli, P., 2006. Comparison of Bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics. 22, 341-345.
Beerli, P., Palczewski, M., 2010. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics. 185, 313-326.
Bickford, D., Lohman, D.J., Sodhi, N.S., NG, P.K.L., Meier, R., Winker, K., Ingram, K.K., Das, I., 2006. Cryptic species as a window on diversity and conservation. Trends Ecol. Evol. 22, 148-155.
Burbrink, F.T., Fontanella, F., Pyron, R.A., Guiher, T.J., Jimenez, C., 2008. Phylogeography across a continent: the evolutionary and demographic history of the North American racer (Serpentes: Colubridae: Coluber constrictor). Mol. Phylogenet. Evol. 47, 274-288.
Burbrink, F.T., Yao, H., Ingrasci, M., Bryson, R.W., Guiher, T.J., 2011. Speciation at the Mogollon Rim in the Arizona Mountain Kingsnake (Lampropeltis pyromelana). Mol. Phylogenet. Evol. 60, 445-454.sdx
Burbrink, F.T., Guiher, T.J., 2014. Considering gene flow when using coalescent methods to delimit lineages of North American pitvipers of the genus Agkistrodon. Zoo. Jour. of the Lin. Soc.
Carstens, B.C., Knowles, L.L., 2007. Estimating species phylogeny from gene-tree probabilities despite incomplete lineage sorting; an example from Melanophlus grasshoppers. Syst. Bio. 56, 400-411.
Carstens, B.C., Pelletier, T.A., Reid, N.M., Satler, J.D., 2013. How to fail at species delimitation. Mol. Ecol. 22, 4369-4383.
Cummings, M.P., Neel, M.C., Shaw, K.L., 2008. A Genealogical Approach to Quantifying Lineage Divergence. Evol. 62, 2411-2422.
Degnan, J.H., Rosenberg, N.A., 2009. Gene tree discordance, phylogenetic inference and the multispecies coalescent. Trends in Ecol. and Evol. 24, 332-340.
Dirzo, R., Raven, P.H., 2003. Global state of biodiversity and loss. Annu. Rev. Environ. Resour. 23, 137-167.
Dupius, J.R., Roe, A. D., Sperling, F.A., 2012. Multi-locus species delimitation in closely related animals and fungi: one marker is not enough. Mol. Ecol. 21, 4422-4436.
Edwards, S.V., 2008. Is a new and general theory of molecular systematic emerging? Evol. 63-1, 1-19.
Ellsworth, D.L., Honeycutt, R.L., Silvy, N.J., Bickham, J.W., Klimstra, W.D., 1994. Historical biogeography and contemporary patterns of mitochondrial DNA variation in White-Tailed Deer from the southeastern United States. Evol. 48, 122-136.
Felsenstein, J., 2004. Inferring Phylogenies. Sinauer Associates, Sunderland, Mass.
Fitch, H.S.F., 1963. Natural History of the Racer Coluber constrictor. Univ. of Kansas, Lawrence. V 15.8, 351-468.
Fujita, M.K., Leaché, A.D., Burbrink, F.T., McGuire, J.A., Moritz, C., 2012. Coalescent-based species delimitation in an integrative taxonomy. Trends Ecol. Evo. 27, 480-488.
Glor, R.E., Warren, D., 2010. Testing ecological explanations for biogeographic boundaries. Evol 65, 673-683.
Hagesawa, M., Kishino, H., Yano, T.A., 1985. Dating of the human ape splitting by a molecular clock of mitochondrial-DNA. J. Mol. Evol. 22, 160-174.
Hewitt, G., 2000. The genetic legacy of the quaternary ice ages. Nature. 405, 907-913.
Hey, J., Nielsen, R., 2007. Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc. Natl. Acad. Sci. USA 104, 2785-2790.
Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A. 2005. Very high resolution interpolated climate surfaces for global land areas. J. Climatol. 25, 1965-1978.
Huelsenbeck, J.P., Andolfatto, P., 2007. Inference of population structure under a Dirichlet process model. Genetics 175, 1787-1802.
Kearse, M., Moir R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., Buxton, S., Cooper, A., Markowitz, S., Duran, C., Thierer, T., Ashton, B., Meintjes, P.L., and Drummond, A.J. (2012). Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics (2012) 28 (12): 1647-1649
Lawson, R., Slowinski, J.B., Crother, B.I., Burbrink, F.T., 2005. Phylogeny of the Colubroidea (Serpentes): New evidence from mitochondrial and nuclear genes. Mol. Phylogen. Evol. 37, 581-601.
Leaché, A.D., Fujita, M.K., 2010. Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus). Proc. R. Soc. B. 277, 3071-3077.
Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 25, 1451-1452.
Lodge, T.E., 2005. The Everglades Handbook: Understanding the Ecosystem. CRC Press.
Myers, E. A., Rodriguez-Robles, J. A., Denardo, D. F., Staub, R. E., Stropoli, A., Ruane, S., Burbrink, F. T., 2013. Multilocus phylogeographic assessment of the California Mountain Kingsnake (Lampropeltis zonata) suggests alternative patterns of diversification for the California Floristic Province. Mol. Ecol. 22, 5418-5429.
Myers, E.A., Burgoon, J., Ray, J.M., Martinez-Gomez, J.E., Matias-Ferrer, N., Mulcahy, D.G., Burbrink, F.T., 2015. Using Bayesian coalescent model testing to resolve taxonomic incongruence at the genus level in the New World snakes Coluber and Masticophis. In review.
Noonan, B.P., Chippindale, P.T., 2006. Dispersal and vicariance: the complex evolutionary history of boid snakes. Mol. Phylogen. Evol. 40, 347-358.
Nosil, P., Harmon, L.J., Seehausen, O., 2008. Ecological explanations for (incomplete) speciation). Trends in Eco. and Evo. 24, 145-156.
O’Meara, B., Jackson, N., Morales, A., Carstens, B., 2015. Phylogeographic inference using approximate likelihoods. Evol. In press.
Padial, J.M., Miralles, A., De la Riva, I., & Vences, M. (2010). Review: The future of integrative taxonomy. Front Zool , 7 : 1-14.
Pella, J., Masuda, M., 2006. The Gibbs and splitmerge sampler for population mixture analysis from genetic data with incomplete baselines. Can. J. of Fisheries and Aquatic Sciences. 63, 576-596.
Phillips, S.J., Anderson, R.P., Schapire, R.E., 2006. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231-259.
Pyron, R. A., Burbrink, F. T., Collins, G. R., de Oca, A. N. M., Vitt, L. J., Kuczynski, C. A., Wiens, J. J., 2010. The phylogeny of advances snakes (Colubroidea), with discovery of a new subfamily and comparison of support methods for likelihood trees. Mol. Phylogen. Evol. 58, 329-342.
de Quieroz, K. 2007. Species concepts and species delimitation. Syst. Bio. 56, 879-886.
Reese, J.S., Noss, R.F., Oetting, J., Hoctor, T., Volk, M., 2013. A vulnerability assessment of 300 species in Florida: Threats from sea level rise, land use, and climate change. PloS One 8, e80658.
Ricketts, T., Imhoff, M., 2003. Biodiversity, urban areas, and agriculture: Locating priority ecoregions for conservation. Conserv. Eco. 8: 1.
Rosenburg, N.A., Nordberg, M., 2002. Genealogical trees, coalescent theory and the analysis of genetic polymorphisms. Nature Rev. Genet. 3, 380-390.
Ruane, S., Bryson, R. W., Pyron, R. A., Burbrink, F. T., 2014. Coalescent Species Delimitation in Milksnakes (genus Lampropeltis) and impacts on phylogenetic comparative analyses. Sys. Bio. 63, 231-250.
Rundle, H.D., Nosil, P., 2005. Ecological speciation. Ecol. Letters. 8, 336-352.
Silvestro, D., Michalak, I., 2011. RaxMLGUI: A Graphical Front-End for RaxML. Org Divers Evol 12, 335-337.
Soltis, D.E., Morris, A.B., McLahchlan, J.S., Manos, P.S., Soltis, P.S., 2006. Comparative phylogeography of unglaciated eastern North America. Molec. Eco. 15, 4261-4293.
Steen, D. A., McClure, C. J. W., Smith, L. L., Halstead, B. J., Dodd Jr., C. K., Sutton, W. B., Lee, J. R., Baxley, D. L., Humpries, W. J., Guyer, C., 2013. The effect of coachwhip presence on body size of North American racers suggests competition between these sympatric snakes. J. of Zoo. 289, 86-93.
Stephens, M., Donnelly, P., 2003. A comparison of Bayesian methods for haplotype reconstruction from population genotype data. Am. J. Hum. Genet. 73, 1162-1169.
Vidal, N., Hedges, S.B., 2005. The phylogeny of squamates reptiles (lizards, snakes, and amphisbaenians) inferred from nine nuclear protein-coding genes. C.R. Bio. 328, 1000-1008.
Walker, D., Avise, J.C., 1998. Principles of Phylogeography as illustrated by freshwater and terrestrial turtles in the southeastern United States. Annu. Rev. Ecol. Syst. 29, 23-58.
Wallach, V., Williams, K. L., & Boundy, J. (2014). Snakes of the World: A catalogue of living and extinct species. CRC Press.
Warren, D.L., Glor, R.E., Turelli, M., 2008. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evol 62, 2868-2883.
Warren, D.L., Glor, R.E., Turelli, M. 2010. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography. 33, 607-611.
Webb, S.D., 1990. Historical Biogeography. In: Myers, R.I., Ewel, J.J., (Eds.), Ecosystems of Florida. University Central FL Press, Orlando, pp. 70-102.
Wilson, L.D., 1978. Coluber Constrictor. Catalogue of American Amphibians and Reptiles. No 218.1-4.
Yang, Z.H., Rannala, B., 2010. Bayesian species delimitation using multilocus sequence data. Proc. Natl. Acad. Sci. USA 107, 9264-9269.
Figure 1: Localities of sampled C. constrictor individuals, and population assignments
Figure 2: Cytb gene tree
Figure 3: Cmos gene tree
Figure 4: Nt3 gene tree
Figure 5: R35 gene tree
| Table 1. DnaSP analysis results for four genes of C. constrictor | |||||||
| Locus | Lineage | Length (bp) | Sample size | Polymorphic sites | Haplotypes/Alleles | Nucleotide Diversity | Tajima’s D |
| cmos | Florida | 531 | 74 | 13 | 11 | 0.00141 | -2.02507 P < 0.05 |
| Continental | 531 | 150 | 21 | 22 | 0.00297 | -1.61082 0.10 > P > 0.05 | |
| cytb | Florida | 415 | 34 | 40 | 8 | 0.01342 | -1.95934 P < 0.05 |
| Continental | 415 | 74 | 51 | 25 | 0.02540 | -1.01643 P > 0.10 | |
| nt3 | Florida | 467 | 74 | 26 | 33 | 0.00670 | -1.28555 p > 0.10 |
| Continental | 467 | 144 | 47 | 83 | 0.00866 | -1.74758 0.10 > P > 0.05 | |
| r35 | Florida | 473 | 74 | 26 | 38 | 0.01102 | -0.89931 p > 0.10 |
| Continental | 473 | 148 | 18 | 16 | 0.00249 | -1.99857 P < 0.05 |
| Table 2. Gsi Analysis results for C. constrictor | ||||||||||
| Lineage | Cytb + P value | Cmos + P value | Nt3 + P value | R35 + P value | Grouped + P value | |||||
| Florida | 0.9263 | 1e-04 | 0.7364 | 2e-04 | 0.2472 | 4e-04 | 0.4246 | 1e-04 | 0.5836 | 1e-04 |
| Eastern | 0.8558 | 1e-04 | 0.1407 | 0.1282 | 0.1086 | 0.2802 | 0.4329 | 1e-04 | 0.3845 | 1e-04 |
| Table 3. MIGRATE-N analysis results for C. constrictor | |||
| Model | Bezier lmL | Model Probability | Model Choice |
| 1 – migration both ways | -7332.25 | 1.44E-53 | 2 |
| 2 – migration F to NE | -7444.90 | 1.72E-102 | 3 |
| 3 – migration NE to F | -7210.58 | 1 | 1 |
| 4 – single population | -7411.59 | 5.04E-88 | 4 |