Demographic history modeling
We used the Generalized Phylogenetic Coalescent Sampler (G-PhoCS) version 1.2.3 (Gronau et al. 2011) to estimate divergence times, effective population sizes and gene flow. Because of the computationally intensive nature of this analysis, we restricted our data set to a maximum of 10 individuals per subspecies, selecting one sample per location. We ran G-PhoCS using the standard MCMC settings described in Gronau et al. (2011) and default parameters with 75,000 burn-in generations and 750,000 additional sampling generations. We assessed adequate mixing and convergence using Tracer. Divergence times and effective population sizes were converted from mutation scale to generations (T) and individuals (Ne) respectively, by assuming an average mutation rate of 10-9 mutations per bp per generation (Kumar and Subramanian 2002). The model implemented in G-PhoCS is conditioned upon a given phylogenetic topology. Thus, we ran G-PhoCS using the topology of the nuclear genomic tree (which is concordant with the Bayesian and MP mitochondrial trees, see the Results section below).