With limited sampling, geographic variation within a single species can be difficult to distinguish from interspecific variation, confounding our ability to draw accurate species boundaries. We argue that thorough sampling and analysis of contact zones between putative taxa can determine if assortative mating or selection against hybrids exists (supporting the presence of two distinct species), or alternatively if mating is random among genotypes and admixture among adjacent populations is gradual and continuous (supporting geographic variation within a single species). Here, we test two alternative hypotheses for two pairs of named taxa at contact zones within the American milksnake (Lampropeltis triangulum) complex. A prior morphological analysis found areas of gradual intergradation among named taxa, and concluded that the taxa represented geographical races of a single polytypic species. In contrast, a subsequent analysis of gene sequence data, but with limited sampling near the contact zones, hypothesized distinct boundaries between species at the contact zones. At the contact zone between proposed species L. triangulum and L. gentilis, we examined a ~700 km-wide transect across the states of Kansas and Missouri, with thorough sampling and reduced-representation genomic-level sequencing, to test the two opposing taxonomic hypotheses. Our transect analyses included examinations of population structure, fixed differences, cline-fitting, and an admixture index analysis. These analyses all supported a gradual and continuous geographic cline across a broad intergrade zone between two geographic forms of L. triangulum, thus providing strong support for a single species in this region (and no support for the recognition of L. gentilis as a distinct species). At a second contact zone between proposed species L. triangulum and L. elapsoides (but variously treated as species or subspecies by different researchers) in Kentucky and Tennessee, we re-evaluated morphological data. In this case, the contact zone analysis indicated sympatry and reproductive isolation of the two taxa, and thus strongly supported L. triangulum and L. elapsoides as distinct species. We conclude that detailed studies of contact zones, based on either genetic or morphological data, are essential for distinguishing intraspecific from interspecific variation in the case of widely and continuously distributed taxa.

(1) Sequencing method: ddRAD sequencing

(2) Bioinformatics pipeline: iPyrad

(3) Analyses:

• Phylogenetic tree (RAxML)
• Population genetics (PCA, sNMF, conStruct)
• Admixture index
• HZAR
• Fixed difference analysis

(4) Data visualization (custom scripts)

Please see README.txt for full details.