Hybridization can have a profound negative effect on population fitness when species exhibit divergence in adaptive traits. The Streamside salamander, Ambystoma barbouri, and the Smallmouth salamander, A. texanum, are closely related species differentiated primarily by breeding habitat and reproduction-related traits, but previous work suggests patterns of hybridization and introgression between them. Here we investigate whether 1) hybridization occurs between the two species in laboratory settings, 2) their divergent reproduction-related traits are retained under a common treatment, and 3) hybrid offspring exhibit similar fitness as their non-hybrid counterparts. We bred conspecific and heterospecific pairs of A. barbouri and A. texanum and raised hybrid and non-hybrid offspring to metamorphosis. Reproduction-related traits such as oviposition location and clutch size remained well differentiated between the two species. Heterospecific pairs hybridized readily, at similar rates to conspecific pairs. The hybrid larvae generally exhibited a maternal effect and, less frequently, intermediate phenotypes with respect to measured traits. Hybrid larvae did not exhibit reduced fitness as measured by survival to metamorphosis. Our results suggest that traits differentiating A. barbouri and A. texanum are likely genetically determined despite the lack of reproductive isolation between them. This suggests that the generally parapatric distribution of the two species, often paired with abrupt transition in traits, might be driven by selection on these traits in local habitats. Further, the maternal effect and consequent lack of intermediate phenotypes observed for several traits likely reduces the negative effect of hybridization in locally adapted populations.

The southern US and northern Mexico serve as an ideal region to test alternative hypotheses regarding biotic diversification. Genomic data can now be combined with sophisticated computational models to quantify the impacts of paleoclimate change, geographic features, and habitat heterogeneity on spatial patterns of genetic diversity. In this study we combine thousands of genotyping-by-sequencing (GBS) loci with mtDNA sequences (ND1) from the Texas Horned Lizard (Phrynosoma cornutum) to quantify relative support for different catalysts of diversification. Phylogenetic and clustering analyses of the GBS data indicate support for at least three primary populations with evidence of recent admixture. The spatial distribution of populations appears concordant with habitat type, with desert populations in Arizona and New Mexico showing the largest genetic divergence. The mtDNA data also support a divergent desert population, but other relationships differ and suggest mtDNA introgression. Genotype-environmental association analyses support divergence along environmental axes. Demographic analyses support a model of allopatric divergence during the Pleistocene followed by secondary contact and gene flow. These results are consistent with inferred paleo-species distribution models. Our results also indicate that caution is warranted when fitting a multispecies coalescent model without introgression to populations that have exchanged genes throughout their diversification history. In sum, our results support allopatric divergence due to Pleistocene climate change, which was followed by secondary contact and widespread genomic introgression. Results also suggest that populations are continuing to diverge along habitat gradients. Finally, the strong evidence of admixture, gene flow, and mtDNA introgression among populations suggests that P. cornutum should be considered a single widespread species.