Ecological divergence and the history of gene flow in the Nearctic milksnakes (Lampropeltis triangulum complex)

Many phylogeographic studies on species with large ranges have found genetic-geographic structure associated with changes in habitat and physical barriers to gene flow. These studies may conclude absence of population structure, lineage structure that indicates unique species have been discovered, or suggest more research is needed prior to delimitation. Comparative risks of delimiting species incorrectly or failing to delimit species are usually not weighed and a more detailed return to these problems with more data often does not occur. With genomic data and better modeling capabilities we can more clearly delimit species by understanding causes of speciation with respect to biogeography and migration between lineages, the location of hybrid zones in relationship to the ecology of parental lineages, and differential introgression of genes between taxa. Here we examine the origins of three Nearctic milksnakes (Lampropeltis elapsoides, L. triangulum, and L. gentilis) using genomic-scale data to better understand the diversification of these taxa previously delimited based on a smaller genetic dataset. Methods that reject pure isolation by distance in favor of environmentally driven reduction in gene flow clearly indicate that all three lineages should continue to be recognized as unique species. These results underscore conspicuous environmental changes that occur between the sister-taxa forming due to changes in habitat from the Great Plains (GP) to the forested regions of the Eastern Nearctic (ENA). This area has been recognized for turnover of reptile and amphibian species but with few phylogeographic studies examining environmental-genetic structure in this region. We show that the two species meeting at the GP/ENA, L. triangulum and L. gentilis, likely formed in the mid-Pleistocene and have maintained partial reproductive isolation over much of this time, exchanging fewer than one migrant/generation, and formed a hybrid zone with differential introgression of loci. We also show that when L. triangulum and L. gentilis are each in contact with the much older L. elapsoides, some limited gene flow has occurred. We conclude that phylogenetic reticulation in this genus, and likely for many other taxa, is common across throughout time. Furthermore, the application of the biological species concept to the whole genome will give rise to grave misunderstandings of the complexities of how species form and remain unique even in the face of gene flow.Please visit this link to view the most updated version of this data submission: https://doi.org/10.5061/dryad.g79cnp5qmElectronic Supplementary Figures Fig. S1. Illustration of the six historical demographic models: IM - Isolation with migration model with constant population size; IMD - isolation with migration model with one demographic change in each species; IMBott - isolation with migration model with two demographic changes in each species generating a demographic bottleneck history; IM-sc - Isolation with a secondary contact migration (following the last glacial maximum) model with constant population size; IMD-sc - isolation with secondary contact (following the last glacial maximum) migration model with one demographic change in each species; IMBott-sc - isolation with a secondary contact migration (following the last glacial maximum) model with two demographic changes in each species generating a demographic bottleneck history. Fig. S2. Population cluster estimation using TESS3r for K = 4, showing groupings corresponding to L. annulata (purple), L. gentilis (green), L. triangulum (red) and L. elapsoides (blue). Fig. S3. DAPC estimates of population structure considering lowest BIC over different GBS sequence assembly filtering strategies with average missing data per individual, number of individuals, and number of loci ranging from 17-48%, 129-159,137-3391, respectively. Fig S4. TESS3r estimates of population structure considering lowest BIC over different GBS sequence assembly filtering strategies with average missing data per individual, number of individuals, and number of loci ranging from 17-48%, 129-159,137-3391, respectively. Fig. S5. Predictions of group membership given ancestral coefficients, where L. gentilis (blue), L. triangulum (green) and L. elapsoides (red) are overwhelmingly supported relative to evenly distributed coefficients (yellow) estimated from TESS3r among all three taxa at southern end of the Mississippi River in Louisiana. Fig. S6. Fit of observed data to different historical demographic models over principal component space. Fig. S7. Using 40 loci between Lampropeltis triangulum and L. elapsoides and nine loci for Lampropeltis triangulum and L. gentilis falling below neutral cline widths showing (A) DAPC spatial predictions for two groups, (B) discriminant function distributions of those groups, and (C) distributions of ancestral coefficients from TESS3r. Fig. S8 – Estimation of cline width and Fst for each locus for both species-pair comparisons: L. gentilis x L. triangulum and L.triangulum x L. elapsoides. Electronic Supplementary Data Data. D1. List of samples and localities used for this research. Data. D2. Fasta files generated for Lampropeltis triangulum, L. gentilis and L. elapsoides from ipyrad. Data. D3. VCF file generated for Lampropeltis triangulum, L. gentilis and L. elapsoides from ipyrad. Data. D4. R code for estimating lineage structure, cline widths from admixture proportions, cline widths for loci, and redundancy analyses (RDA) and artificial neural networks (ANN) to examine the effect of space and ecology on genetic structure of Lampropeltis triangulum, L. gentilis and L. elapsoides. Data. D5. Priors used for simulation of historical demographic parameters in PipeMaster for Lampropeltis triangulum, L. gentilis and L. elapsoides