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.

Methods
Electronic 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.

针对分布范围广泛的物种开展的诸多系统地理学（phylogeographic）研究均发现，其存在与栖息地变化及基因流（gene flow）物理屏障相关的遗传-地理结构（genetic-geographic structure）。此类研究可能得出种群结构（population structure）缺失、支系结构（lineage structure）暗示独特物种已被发现的结论，或是提出在进行物种界定（species delimitation）前需开展更多研究。研究者通常不会权衡错误界定物种或未完成物种界定的相对风险，且往往不会借助更多数据对这些问题展开更深入的重新分析。借助基因组（genomic）数据与更完善的建模能力（modeling capabilities），我们可通过以下方式更清晰地完成物种界定：解析与生物地理学（biogeography）、支系间迁移相关的物种形成（speciation）成因，明确杂交带（hybrid zone）相对于亲本支系（parental lineage）生态的分布位置，以及类群（taxon）间基因的差异性渐渗（introgression）。本研究利用基因组规模数据，对3种新北区乳蛇（Lampropeltis elapsoides、L. triangulum及L. gentilis）的起源展开分析，以更深入理解此前基于较小规模遗传数据集界定的上述类群的分化历程。那些拒绝“纯距离隔离（isolation by distance）”假说、转而支持“环境驱动的基因流降低”的分析方法明确表明，这3个支系仍应被视为独立物种。研究结果凸显了姐妹类群（sister-taxa）间显著的环境变化——这类变化源于大平原（Great Plains, GP）至新北区东部（Eastern Nearctic, ENA）森林区域的栖息地转变。该区域以爬行动物和两栖动物类群更替著称，但鲜有系统地理学研究探讨该区域的环境-遗传结构。研究表明，在大平原/新北区东部交界地带相遇的L. triangulum与L. gentilis可能形成于中更新世，且在绝大多数时间里维持着部分生殖隔离，每代仅发生不足1次的迁移交换，并形成了存在位点差异性渐渗的杂交带。此外我们还发现，当L. triangulum与L. gentilis分别与更为古老的L. elapsoides接触时，均发生了有限的基因流。本研究得出结论：该属的系统发育网状演化（phylogenetic reticulation）现象，以及诸多其他类群的类似情况，在演化历程中普遍存在。此外，将生物学物种概念（biological species concept）应用于全基因组分析时，极易对“物种如何形成并在基因流存在的情况下维持独特性”这一复杂过程产生严重误解。

方法
电子补充图表
Fig. S1 六种历史种群动力学模型示意图：IM——恒定种群大小的迁移隔离模型；IMD——每个物种各发生1次种群动态变化的迁移隔离模型；IMBott——每个物种各发生2次种群动态变化、形成种群瓶颈（demographic bottleneck）历史的迁移隔离模型；IM-sc——恒定种群大小的（末次冰盛期（last glacial maximum）后）二次接触（secondary contact）迁移隔离模型；IMD-sc——每个物种各发生1次种群动态变化的（末次冰盛期后）二次接触迁移隔离模型；IMBott-sc——每个物种各发生2次种群动态变化、形成种群瓶颈历史的（末次冰盛期后）二次接触迁移隔离模型。
Fig. S2 利用TESS3r对K=4时的种群聚类估算结果，显示的类群分别为L. annulata（紫色）、L. gentilis（绿色）、L. triangulum（红色）及L. elapsoides（蓝色）。
Fig. S3 基于不同GBS序列组装过滤策略下的最低贝叶斯信息准则（BIC）值估算种群结构的DAPC结果，各策略对应的个体平均缺失率、个体数及位点数分别为17%~48%、129~159、137~3391。
Fig. S4 基于不同GBS序列组装过滤策略下的最低贝叶斯信息准则（BIC）值估算种群结构的TESS3r结果，各策略对应的个体平均缺失率、个体数及位点数分别为17%~48%、129~159、137~3391。
Fig. S5 基于祖先系数的类群归属预测结果：在路易斯安那州密西西比河南端的3个类群中，相较于TESS3r估算的均匀分布系数（黄色），L. gentilis（蓝色）、L. triangulum（绿色）及L. elapsoides（红色）的归属支持度显著更高。
Fig. S6 主成分空间下观测数据与不同历史种群动力学模型的拟合度。
Fig. S7 针对符合中性渐变群宽度阈值的位点：L. triangulum与L. elapsoides间的40个位点、L. triangulum与L. gentilis间的9个位点，结果展示（A）两组类群的DAPC空间预测结果，（B）两组类群的判别函数分布，（C）TESS3r估算的祖先系数分布。
Fig. S8 两个物种类对（L. gentilis × L. triangulum与L. triangulum × L. elapsoides）的每位点渐变群宽度与Fst值估算结果。

电子补充数据
数据D1 本研究使用的样本及采样位点列表。
数据D2 基于ipyrad为Lampropeltis triangulum、L. gentilis及L. elapsoides生成的FASTA格式文件。
数据D3 基于ipyrad为Lampropeltis triangulum、L. gentilis及L. elapsoides生成的VCF格式文件。
数据D4 用于估算支系结构、基于混合比例计算渐变群宽度、估算位点点渐变群宽度，以及冗余分析（RDA）与人工神经网络（ANN）以探究空间与生态因素对Lampropeltis triangulum、L. gentilis及L. elapsoides遗传结构影响的R代码。
数据D5 利用PipeMaster为Lampropeltis triangulum、L. gentilis及L. elapsoides模拟历史种群动力学参数时使用的先验分布。