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.