Introduction
Understanding the role of biogeography in areas that generate abundant diversity is of great interest to researchers of evolution and ecology. Finding the origins of diversity are also integral to conservation practicesthat prioritize areas in need of protection, as well as cataloging the Earth’s biodiversity (Dirzo & Raven, 2003). Biogeographically varied regions with high endemism can illuminate diversification processes by examining the ways in which these environments generate their greater rates of biodiversity. This biodiversity is well-documented in areas such as the Mississippi River Embayment (Bradley, 2010; Burbrink 2008; McKelvy 2017) and the Dry Diagonal in Brazil (Fonseca 2018; Oliviera 2017; Oliviera 2015), which have transitional ecotones with wide-ranging generalist taxa. However, understanding processes at play for generating community-level diversity can be inhibited by unrecognized cryptic species (Bickford et al, 2006), thereby making molecular species delimitation important for assessing underlying diversity of heterogeneous biogeographical regions.
The Florida peninsula is a hotspot of biodiversity within North America. It ios the only southeastern state to have endemic species in each of its counties (Estill and Cruzan, 2001).This is likely due its subtropical climates, peninsular form combined with distance from other land masses in the south, and complex geological history, including isolation during the Pliocene and Pleistocene from historically higher sea levels (Reese, 2013; Lodge, 2005; James, 1961). Southern Florida’s vascular plants are just over 60% tropical, with an additional 9% endemic to Florida, while most of the freshwater and terrestrial vertebrates that live in SouthernFlorida mostly temperate North America in origin (Lodge, 2005). This creates an interesting mix of temperate fauna and tropical and subtropical flora species within Florida’s ecosystem. 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 northern boundary of the Florida peninsula (Webb, 1990). These conditions may have generated species with unique morphological characteristics found only in ecoregions specific to Florida, such as those found in the Florida Sand Pine Scrub (Ricketts & Imhoff, 2003). It is estimated that 40 – 60% of the Pine Scrub species are endemic (Webb, 1990). The number of endemic plant species is estimated to be as high as 385, or 11% of total fauna, with at least 22 listed as endangered or threatened by the State of Florida (Myers and Ewer, 1990). Endemic vertebrates found in the Florida Sand Pine Scrub are the Florida Scrub Lizard (Sceloporous woodi) Florida Sand Skink (Plestiodon reynoldsi), the blue-tailed mole skink (Plestiodonegregius lividus) and the Florida Scrub Jay (Aphelocoma coerulescens; Myers and Ewer, 1990). Additionally, there are several taxa with independent phylogeographic lineages found only in peninsular Florida, such as the White-Tailed Deer (Odocoileus virginianus) (Ellsworth et al., 1994), the Oldfield Mouse (Peromyscus polionotus) (Avise et al., 1984; Degner 2007; Van Zant and Wooten 2007), and several turtle species (Walker and Avise, 1998). Studies of snakes in Florida have shown similar populations divergences. A multi-locus study by Burbrink & Guiher (2015) showed significant ecological niche differences between the populations of the Florida Cottonmouth (Agkistrodon conanti) and the Northern Cottonmouth (Agkistrodon piscivorus), with a divergence dated to 2.5-1.9 million years ago (MYA) ago. Krysko et al (2016) found support for divergence in the Eastern Indigo Snake (Drymarchon couperi) in Florida (but see Folt 2018). Additionally, population divergence in Florida is also supported in the yellow-bellied kingsnake (Lampropeltis calligaster) (McKelvy 2017).
The North American racer (Coluber constrictor) 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 mtDNA identified six lineages across North America. The most deeply divergent lineage was found in peninsular Florida (Burbrink et al, 2008), which was dated to be in the late Miocene/Early Pliocene, about 6 MYA. This study relied on a single mitochondrial DNA (mtDNA) locus to infer phylogeographic history of this taxon. It may be that these studies only reflect the history of the locus and not that of the species as a whole, as gene trees reflect local optimums and have the same or older divergences than species trees (Degnan and Rosenberg, 2009). Also, gene tree histories within the same individual may differ from each other due to incomplete lineage sorting, horizontal gene transfer, gene duplication and branch length heterogeneity (Carstens and Knowles, 2007; Edwards 2009). Myers et al (2017) investigated the genus Coluber to resolve the taxonomic nomenclature between Coluber and the sister taxon Masticophis. Multi-locus coalescent analyses showed deep divergence in the species tree between peninsular Florida and eastern North America populations of Coluber.
The processes responsible for driving divergence in C. constrictor is unclear but two candidate hypotheses are likely possible: 1) elevated seas levels may have isolated taxa, and 2) ecological speciation from continental US, driven by Florida’s unique ecosystems may have reduced gene flow (Rundle & Nosil, 2005). Knowing the timing of divergence, assessing gene flow, and discovering dynamic contact zones may help understand possible causes of divergence into unique niches. Here, we examine lineage divergence of C. constrictor between peninsular Florida and the eastern continental North America populations using coalescent-based methods to ascertain the phylogeographic history. Coalescent-based approaches provide a testable method to understand 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, 2014). 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 (Dupius, 2012) which may obscure population delineations. Parameters such as effective population size and divergence times can also be simultaneously estimated (Yang & Rannala, 2010).
First, given potential problems with single locus studies, we determine if the mtDNA gene tree is in accordance with multiple nuclear DNA (nDNA) gene trees. Next, we look at evidence of gene flow between populations in peninsular Florida and the lineage in eastern continental United States. And finally, we determine if the timing of divergence in these lineages could be due to vicariant events and/or differences in ecological niches. The results show significant genetic divergence between these two populations, and add to the growing body of data that show Florida to be a biodiversity hotspot for generating species diversity.
Materials and Methods
Taxon Sampling and Sequence Data
We obtained 113 tissue samples of C. constrictor from museum collections and field sampling, geographically ranging from south Florida to Virginia, and west to southern Louisiana (FIG. XX, Appendix 1). DNA was extracted using Qiagen DNeasy kits from shed skin, liver, or muscle tissue. One mtDNA locus (CYTB) and three protein-coding nuclear DNA loci (nDNA; CMOS, R35, NT3), 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, H14910 and THRSN2 primers (Burbrink et al., 2008) were used, and H19410 for cycle sequence reactions, in one direction. Primers for CMOS were S77 and S78 (Lawson, 2005), and for R35 were GeneF and GeneR (Vidal & Hedges, 2005), and for NT3 were F3 and R4 (Noonan & Chippindale, 2006). 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. DnaSP v5 (Librado & Rozas, 2009) was used to calculate polymorphic sites (S), haplotype number (h), haplotype diversity (Hd)diversity (π), Tajima’s D and it’s P value for each locus.
Population Assignment and Speciation
To assess the number of populations in C. constrictor in the Southeastern United States, we used STRUCTURE 2.3.2 (Pritchard et al, 2000a), which assigns individuals to inferred populations using a Bayesian clustering algorithm, assuming Hardy-Weinberg and linkage equilibrium. We explored a range of number of populations (K) from 1 to 4 and running 5 replicate analysis for each. Each independent run implemented 100,000 generations following a burn in of 10,000 generations, assuming a non–linkage model and uncorrelated allele frequencies. We chose the best value of K based on the rate of changes in the log-probability of data between successive K values, ∆K (Evanno et al, 2005), using STRUCTURE HARVESTER v0.6.9 (Earl & vonHoldt, 2012).
Bayesian Phylogenetics and Phylogeography (BPPv3; Yang and Rannala, 2010;Yang 2015) was used to determine if the populations recognized by STRUCTURE represent distinct species. This method uses multiple independent loci in a coalescent framework with the species phylogeny, assessing lineage sorting due to ancestral polymorphism. Fine-tuning parameters were set to the recommended values of swapping rates ranging between 0.30 and 0.70. All loci were included in the analysis, and three combinations of priors for ancestral populations size (θ) and root age (τ0) were assessed 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.
Gene Tree Analysis and Species Tree Estimation
Phylogeographic structure for each gene was estimated using maximum likelihood in RaxmlGUI 1.1 (Silvestro and Michalak, 2011). All trees were rooted with the coachwhip (Masticophis flagellum), the sistertaxon to the species C. constrictor (Myers et al, 2017). Nucleotide substitution was modeled with GTRGAMMA as assessed by jModelTest (Darriba et al, 2012), and 1000 nonparametric bootstrap replicates were run to assess nodal support. Clades with values > 75% were considered well supported (Felsenstein, 2004).
We ran *BEAST using the mtDNA gene CYTB. This analysis requires a priori assignment of individual alleles to a species before estimating the relationship, and we therefore made assignments based on STRUCTURE results. All *BEAST analyses were run for 100,000,000 generations and sampled every 10,000 generations. We assessed convergence of the MCMC runs and ESS values (>200) using TRACER. Analyses with nDNA did not converge despite repeated attempts, therefore mtDNA alone was used. The first 10% of sampled genealogies were discarded as burn-in.
Incomplete Lineage Sorting
We calculated the genealogical sorting index (gsi) for each gene (Cummings et al, 2008) which gives the degree of ancestry of the given populations based on a gene tree. The result of this index is based on a scale of 0 – 1, with 0 indicating complete non-exclusivity, and 1 indicating monophyly of the gene. The maximum likelihood gene trees from Raxml for the two populations 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.
Historical Demography
We built Bayesian skyline plots (Drummond et al, 2005) for each lineage in BEAST 1.8.0. We used CYTB only for the same reasons as previously stated; lack of convergence with any analysis using nDNA. We used an mtDNA substitution rate of 9.8 x 10^-9 substitutions/site/year to calibrate the molecular clock. We ran 1.0 x 10-8 generations sampling every 5,000 generations. Parameter convergence, stationarity and ESS values (>200) were visually assessed using TRACER. (FIG XX)
Estimation of Divergence Dates and Migration
To estimate divergence time, ancestral population sizes, migration, and gene flow we used the program IMa2 (Hey and Nielson, 2007; Hey 2010) on the mtDNA locus, and on 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). We 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, with a generation time of two years (Ernst and Ernst, 2003), and a substitution model of HKY (Hasegawa et al, 1985) as the closest model fit this program offers, 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 = 40, m1, m2 = 1, and τ = 6. 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 were tested against the five parameters of the full model, such that Ne1 ≠ Ne2 ≠ NeA, and m1 ≠ m2, and m1 and m2 were evaluated where effective migration = 0. The 16 nested models tested that some or all of these parameters were equivalent to the full model using the 2LLR test from a chi-square distribution.
Our previous analyses showed C. constrictor is composed of two major clades with low levels of migration between them. Our demographic analyses showed a population expansion for both lineages. We used an approximate Bayesian computation (ABC) approach (Beaumont 2010) to test the modeling of likely scenarios.
The seven scenarios are as follows. Panmixia, as a control scenario for comparison (P), Isolation without migration (Is), Isolation without migration and demographic change in the continental population (IsDcont); the demographic change here and in the following scenarios being initial contraction after a population split, then expansion; Isolation without migration and demographic change in both populations (IsD) Isolation with migration (Im), Isolation with migration and demographic change in the continental population (ImDcont), Isolation with migration and demographic change in both populations (ImD).
We performed 1 000 000 simulations for each model using MSABC (Pavlidis et al 2010). We used an R script to sample parameters from prior distributions and call MSABC. Parameter prior distributions were based on results of *BEAST, BSP and IMa2 analyses. We implemented priors in demographic quantities and transformed the ms-scaled parameters using equations from the ms manual (Hudson 2002). We called MSABC one time for each of the four genes, using the same number of samples and length of loci. For the mitochondrial gene, we sued one-fourth of the sampled population size. Five summary statistics for each gene were used: nucleotide diversity and Tajima’s D for each populations, and the FST between populations. We transformed the observed sequence data into ms-like files using the fas2ms perl script (from MSABC package) and calculated the same summary statistics for each locus. To estimate posterior probabilities and model support (‘postpr’ function), we used R package ‘abc’ (Csilléry et al 2012). We set the tolerance to 0.0001 and implemented both the nonlinear neural network regression and rejection regression methods. We summarized simulated model fit to the observed data using a principal component analysis (PCA) calculated in R.
Ecological Niche Modeling
To construct Ecological Niche Models (ENMs) we used 19 bioclim variables describing variation in precipitation and temperature from the Worldclim data (Hijmans et al., 2005) taken at 30-s resolution. Using the locality points from our sampled individuals, we ran 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 considered identical if the observed ENMs are not any more different than randomly drawn pairs of the analyzed samples.
We 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, using the I statistic (Warren et al, 2010).
Results
DNA polymorphism
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). All 113 samples were amplified and sequenced for CMOS, 108 for CYTB, 109 for NT3, and 111 for R35, for a final data matrix with 97% completion. DnaSP analysis showed lower allelic diversity in the Florida population for three out of four of the genes, and greater diversity in the eastern continental population except for R35 (Table 1). All genetic sequences have been deposited in GenBank (Accession numbers: YYY-ZZZ).
Incomplete Lineage Sorting and Population Structure
Gene trees were inferred with maximum likelihood in RaxML, with blue branches indicating eastern continental individuals, and red indicating Florida individuals, as assessed by Structurama (FIGS XX).
The gsi analyses shows that all loci were significant with respect to sorting for the Florida lineage, where only CYTB and R35 are significant for the eastern continental lineage (Table XX). In Florida, CYTB approaches monophyly (0.93), followed by CMOS (0.74) and R35 (0.42), with NT3 being the least monophyletic (0.25).
STRUCTURE detected two population clusters (k=2), with a break at northern peninsular Florida. (FIGS X.X)
Speciation and Divergence
Every run of BPP delimited the two populations, Florida and the eastern continental, into separate species with a probability of 1 under all three models of divergence.
Species tree estimation and migration
IMa2 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. Sequences were trimmed to delete the non-informative base pairs at the ends, and some sequences were discarded as not informative enough once trimmed for this purpose, with a total of 1871 base pairs being used, and a final base pair length of 531 for CMOS, 403 for CYTB, 466 for NT3, and 471 for R35. This gave a final data matrix that was 84% of the original sequences, with the lowest number of individuals used for any gene being 90 for CYTB. The mean times across all runs of the most recent common ancestor were as follows: CMOS – 8.6 MYA, CYTB – 7.6 MYA, NT3 – 3.2 MYA, and R35 – 6.0 MYA. The estimate of divergence is consistent across runs (1.4 Ma, 95% CI = 0.83-2.6 Ma) dating it to the Lower Pleistocene.
The IM model with the lowest AIC was Model 2, which showed different ancestral and extant population sizes with equal migration between the two taxa (m1 = 0.23, 95% CI = 0.06-0.48; m2 = 0.23, 95% CI = 0.10-0.41). The model with the second lowest AIC was the Full Model, with unequal population sizes and unequal migration rates. The ΔAIC was 1.99. The predicted effective population size of the ancestral population (NeA; 228,000; 95% CI = 140,000-340,000) was significantly smaller than the Florida population (NeF; 490,000; 95% CI = 349,000-673,000) with the eastern continental lineage being the largest (NeC; 728,000; 95% CI = 560,000-934,000). It also showed similar migration rates between the extant populations with 1.68 gene copies per generation going into Florida from the continental population, and 2.56 in the reverse direction.
*BEAST analyses for CYTB showed the gene tree divergence to date to 5.17 MYA with a 95% HPD interval of 3.8-6.7 MYA, and the species tree divergence dating to 137,000 years ago with a 95% HPD interval of 0.065-2.14 MYA.
Approximate Bayesian computation analysis showed the highest support for the scenario of Isolation with Migration with demographic change in both populations. Summary statistics of our observed data were within the bounds of the summary statistics for demographics changes simulations in the PCA predictive plots, confirming good model fitting. (Table XX)
ABC: gene copies per generation going into Florida from the continental population is 2.13, while the reverse direction shows 1.34 (Table XX).
Both the Florida and continental populations show consistent large ancestral population sizes with a strong contraction and population expansion in the late Pleistocene (See Figures 5a and 5b).
Ecological Niche Modeling
The results from I (0.40) and Schoener’s D (0.18) statistics were both similarly significant, so we used the statistic I values for calculations to avoid redundancy and for brevity, as it tends to have the greater variation of the two (Glor & Warren, 2010). The Identity Test showed that the ENMs are not identical (t = 322.7, df = 99, P value < 2.2e-16, 95 percent confidence interval: 0.95, 0.96), calculated by using I. The range break blob test supports these results (t = 15.6, df = 99, P value < 2.2e-16, 95 percent confidence interval: 0.56, 0.61), using the same observed I statistic. This suggests the observed climates have greater difference between Florida and the eastern continental United States than would be expected by randomized geographic boundaries.
Discussion
The biogeography of Florida has repeatedly been shown to be a species generator since its emergence from the Eocene, with theoretical pressures ranging from glacial refugia to edaphic and climatic differentiation from the rest of the continent. This study also shows the signal of that influence on C. constrictor. There are two distinct, divergent populations; the first is in peninsular Florida, and the second is the adjacent population in the remaining range of the continent. The contact zone of these two populations occurred at the junction of the Florida peninsula with the continent (Fig. XXX). The divergence time of these populations is confidently within the Pleistocene, although IMa2 and *BEAST have differing results as to when during the Pleistocene this occurs.
We found significant differences in the ecological niches of the eastern continental and Florida populations. This is unsurprising as there are several climatic differences between the temperate climate of the continent and the sub-tropical climate of Florida. The peninsula has high humidity from being surrounded by sea and ocean on three sides, and the dry/wet seasonality of Florida is markedly distinct compared to the temperate seasons of the continent (Myers and Ewel, 1990). Florida rarely experiences frost; Tallahassee in the northernmost part of the state had the most days of freezing between 1930 and 1979 – 31 days, whereas Homestead in the southern portion of the state only had 2 (Myers and Ewel, 1990).
A confounding effect regarding the long time intervals between the species and gene tree divergences seen in both the IMa2 and BEAST analyses could be the impact of the large ancestral effective population size as addressed by Edwards and Beerli (2000), where Τ – τ = t; T is time in generations since the gene divergence occurred, τ is the population divergence, and t is the difference with an expected value of 2Ne generations of the ancestral population. Substituting our BEAST analysis results, this gives us (5,170,000/2) – (137,000/2) = 2,516,500. For comparison, expected t from ABC analysis is 2(947,438) = 1,894,876. While not completely congruent, this is well within the 95% CI: 336,434-3,404,476. This highlights the impact a large ancestral population size can have on the timing of divergence, completely relevant in considering the wide continental range and generalist nature of C. constrictor.
It has been suggested that Florida endemics are only slightly morphologically different from continental species due to recent divergences (James, 1961). Even so, there is a significant difference in the morphology of C.constrictor. 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 interspecies resource partitioning with C. flagellum, as C. constrictor was found to be smaller when C. Flagellum was present (Steen et al, 2013). There could be a confounding factor such as the effect of Florida’s ecosystem, as this effect was not found in other areas of the country where the two snakes are sympatric.
Historically, mtDNA has been the primary locus for phylogeographic studies and often still provides the basic structure for subsequent studies (Burbrink et al, 2008). Few researchers have followed up these primary studies within a multilocus coalescent framework to retest previous outcomes (Myers et al, 2013). We found here that the population structure of the mtDNA analysis by Burbrink et al, 2008 holds up well under multilocus analyses although for snakes, though this is not always the case (Ruane et al, 2014; Myers et al, 2013). The timing of divergence of the gene trees are also in accordance with the previous study, although it the species tree divergence is much younger than the gene trees. Our analyses reveal that gene flow seems to be a contributing factor in the evolutionary history of the C. constrictor, testing of which has not been done before.
It is likely that diversification in racers in the southeastern United States is a multi-part process: Initially populations were likely isolated by vicariant events associated with glacial minima and maxima, then differing ecological niches occupied by these populations after receding sea levels may have reinforced the divergence event as the independent populations once again came into contact (Nosil et al, 2008). Several studies have shown the influence of a “peninsular effect” on Florida populations, in which the peninsular is a site of refuge during glacial maxima, after which it becomes a source of population expansions (Estill and Cruzan, 2001; Fetter and Weakley, 2019). Additionally, the changing ecological niches of Florida may act as a population filter, which would require adaptive characteristics of an expanding population (Fetter and Weakley, 2019).
Additional studies on a genomic scale should be done as it could be crucial in elucidating the true evolutionary history of this species in Florida and across its range of continental North America. This would add to the growing body of knowledge of how biota respond to biogeographical influences over time. Many areas of conservation research focus on rare and endangered species, but generalists and wide-ranging taxa are important in setting evolutionary baseline responses to the impact of biogeography on species history.