DOI: 10.5061/dryad.g79cnp5qm
title: Ecological divergence and the history of gene flow in the Nearctic milksnakes (Lampropeltis triangulum complex)

description: Many phylogeographic studies on species with large ranges have found
 genetic-geographic structure associated with changes in habitat and
 physical barriers preventing or reducing gene flow. These interactions
 with the geographic space, contemporary and historical climate, and
 biogeographic barriers may have complex effects on population genetic
 structure and speciation. While allopatric speciation at biogeographic
 barriers is considered the primary mechanism for generating species, more
 recently it has been shown that parapatric modes of divergence may be
 equally or even more common. With genomic data and better modeling
 capabilities we can more clearly define causes of speciation in relation
 to biogeography and migration between lineages, the location of hybrid
 zones with respect 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. Results from artificial neural networks
 show that a mix of a strong biogeographic barrier, environmental changes,
 and physical space has driven geographic structure in these taxa. These
 results underscore conspicuous environmental changes that occur as sister
 taxa diverged near the Great Plains into the forested regions of the
 Eastern Nearctic. This area has been recognized as a region for turnover
 for many vertebrate species but as we show here, the contemporary boundary
 does not isolate the sister species, L. triangulum and L. gentilis. These
 two species likely formed in the mid-Pleistocene and have remained
 partially reproductively isolated over much of this time, exchanging fewer
 than one migrant/generation with hybrid zones showing 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. Given the strong agreement between nuclear
 and mtDNA genomes, along with estimates of fundamental niche, we suggest
 that all three lineages should continue to be recognized as unique
 species. Furthermore, this work emphasizes the importance of considering
 parapatric modes of divergence and differential allelic introgression over
 a complex landscape when considering mechanisms of speciation.  
description: <b><i>Electronic
 Supplementary Text</i></b> Text S1 –
 Supplementary information providing additional detail for computational
 methods used. <b><i>Electronic
 Supplementary Figures</i></b> Fig. S1 –
 Location of samples for grouping into numbered (and colored) populations
 for addressing longitudinal allele surfing. Fig. S2 – A)
 results from DAPC showing Bayesian information criteria (BIC) against
 number of clusters, B) discriminant analysis of principal components
 (DAPC) scatter plots, and B) DAPC membership probabilities for each
 individual sample by species. Fig. S3 –
 Correlation of admixture coefficients between TESS3r and SNMF for
 <i>Lampropeltis gentilis, L. triangulum</i>, and <i>L.
 elapsoides</i>. Fig. S4 – Pie
 charts over a map of North America (partial) showing admixture
 coefficients generated from SNMF for the milksnakes examined
 here. Fig. S5 – 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 19-48%,
 129–159,137–3391, respectively. Fig S6 – 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 19-48%,
 129–159,137–3391, respectively. Fig. S7 –
 Predictions of group membership given ancestral coefficients, where
 <i>Lampropeltis gentilis</i> (green), <i>L. triangulum
 </i>(blue), <i>L. annulata</i> (purple) and <i>L.
 elapsoides</i> (red) are supported relative to those coefficients
 estimated from TESS3r equally supporting <i>L. gentilis</i>
 <i>L. triangulum, </i>or <i>L. elapsoides</i>
 (yellow) at southern end of the Mississippi River in Louisiana.
 Fig. S8 –
 The four groups of milksnakes predicted by DAPC, red =
 <i>Lampropeltis triangulum</i>, blue = <i>L.
 elapsoides</i>, green = <i>L. gentilis</i>, purple =
 <i>L. annulata, </i>and red dots with yellow halos showing
 individuals predicted as <i>L. triangulum</i> but with
 <i>L. gentilis</i> mtDNA. Fig. S9 – Map
 across all milksnakes demes showing genetic diversity estimated using
 EEMS. Fig. S10 – Fit
 of observed data to different historical demographic models over principal
 component space. Fig. S11 – Using
 40 loci between <i>Lampropeltis triangulum </i>and <i>L.
 elapsoides</i> and nine loci for <i>L. triangulum
 </i>and<i> L. gentilis </i>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. S12 –
 Estimation of cline width and
 <i>F</i><sub>ST</sub> for each locus for both
 species-pair comparisons: <i>Lampropeltis gentilis</i> x
 <i>L. triangulum </i>and <i>L. triangulum </i>x
 <i>L. elapsoides</i>. Fig. S13 –
 Location of mtDNA clades and GBS groups relative to the 50% admixture
 clines (yellow  =  cline between <i>Lampropeltis gentilis</i>
 and <i>L. triangulum</i>; pink = cline  between <i>L.
 triangulum</i> and <i>L. elapsoides</i>) estimated using
 GBS data. <b><i>Electronic
 Supplementary Data </i></b> Data. D1 – List
 of samples and localities used for this research. Data. D2 – Fasta
 files generated for <i>Lampropeltis triangulum, L. gentilis, L.
 annulata,</i> and <i>L. elapsoides</i> from
 ipyrad. Data. D3 – VCF
 file generated for <i>Lampropeltis triangulum, L.
 gentilis</i>, <i>L. annulata,</i> and <i>L.
 elapsoides</i> from ipyrad filtered. 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 <i>Lampropeltis triangulum,
 L. gentilis, L. annulata,</i> and <i>L.
 elapsoides.</i> Data. D5 –
 Priors used for simulation of historical demographic parameters in
 PipeMaster for <i>Lampropeltis triangulum, L. gentilis</i> and
 <i>L. elapsoides.</i> Data. D6 –
  Milksnake genetic assignment file results from DAPC and TESS3r.
 Data. D7 –
 Raw VertNet (http://vertnet.org) downloads, R code, and cleaned and
 partitioned localities for <i>L. triangulum, L. gentilis</i>,
 and <i>L. elapsoides </i>used in ecological niche modeling.
 Note: R scripts to run this are a separate file and not zipped with these
 datasets.