Inferring the history of divergence between species in a framework that permits the presence of gene flow has been crucial for characterizing the “gray zone” of speciation, which is the period of time where lineages have diverged but have not yet achieved strict reproductive isolation. However, estimates of both divergence times and rates of gene flow often ignore spatial information, for example when considering the location and width of hybrid zones with respect to changes in environment between lineages. To connect phylogeographic estimates of lineage structure, migration, historical demography, and timing of divergence with hybrid zone dynamics, we used population genomic data from the North American ratsnake complex (Pantherophis obsoletus). We examined the spatial context of diversification by linking migration and timing of divergence to the location and widths of hybrid zones. Artificial neural network approaches were applied to understand how landscape features and past climate have influenced population genetic structure among these lineages. We found that rates of migration between lineages were associated with the overall width of hybrid zones. Timing of divergence was not related to migration rate or hybrid zone width across species pairs but may be related to the number of alleles weakly introgressing through hybrid zones. This research underscores how incomplete reproductive isolation can be better understood by considering differential allelic introgression and the effects of historical and contemporary landscape features on the formation of lineages as well as overall genomic estimates of migration rates through time.

Supporting Information Material 1: Supplemental Methods and Results

Data S1. All Pantherophis obsoletus samples with locality information sequenced.

Data S2. R code for estimating groupings, admixture and locus clines, and RDA and ANN estimates of ecological and spatial predictors of genetic variation.

Data S3. VCF files for use with R code.

Data S4. Locality and sample data formatted for use with R.

Data S5. Unfiltered Fasta files.

Geographic Localities

R scripts

Supporting Information Material 1: Additional Methods

Table S1 – Priors and summary stats for PipeMaster; estimates of migration, timing of divergence, and historical demography

Table S2 – PHRAPL Results (open with R read.table() or Excel)

Fig. S1 – Delimiting populations using sparse nonnegative matrix factorization (SNMF) and spatial Principal Component analysis (sPCA) showing the location of four geographically distinct lineages.

Fig. S2 – DAPC and TESS estimates of population structure over K = 2 to 5.

Fig. S3 – PipeMaster model fit of observed data to simulations under IM, IMD, and IMD-LGM  models for PC1-PC10.

Fig. S4 – Expansion times estimated from PipeMaster for all adjacent taxon pairs.

Fig. S5– Results from PipeMaster after removing hybrids showing estimates and directionality of gene flow (A), divergence times between groups and taxon pairs and the Mississippi River (P.alleghaniensis/ P.spiloides vs P.obsoletus/P.bairdi) (B), and changes in population sizes over time (C).

Fig. S6 – Prediction of cline width under neutrality (black line through time) given divergence times for each species pair: A) Pantherophis alleghaniensis/P. spiloides, B) P. spiloides/P. obsoletus, and C) P. obsoletus/P. bairdi. Colored vertical lines represent times of divergence between species pairs and horizontal colored lines indicate actual cline widths.