Population Structure
To assess population structure in each species, we explored the genome-wide SNP data using STRUCTURE v2.3.4 (Pritchard et al.2000). We ran STRUCTURE for K values 1 to 10 with 5 replicates per K , where each replicate is a different sample of unlinked SNPs, subsampled from the same linked SNP dataset. We ran STRUCTURE for 500,000 generations and discarded the first 10% as burn-in. The data were modeled assuming admixture and correlated allele frequencies between populations (Falush et al. 2003), while all other parameters were kept as their default. For the sake of comparison to results reported by Fernandez et al. (2021), Structure Harvester (Earl & vonHoldt 2012) was then used to assess K using the Evanno method (Evanno et al. 2005). We acknowledge that many empirical studies have applied Evanno’s K and interpreted the results as a measure of model fit, but the popularity of a given methods should not be mistaken for its appropriateness. Evanno’s K lacks properties that a parameter in evolutionary genetics should ideally possess, for example it is extremely susceptible to uneven sampling (Puechmaille 2016) and reproducible inference is challenging (Novembre 2016). For these reasons, we follow Pritchard et al . (2000) in treating STRUCTURE as a tool for data exploration rather than phylogeographic inference. Our inferences are derived from the results of the model-based analyses described below.