Population structure

Beagle files with the nuclear genome positions of single nucleotide polymorphisms (SNPs) were produced by ANGSD (Korneliussen et al., 2014) using the following options: -GL 1 -doGlf 2 -minMaf 0.05 -C 50 -baq 2 -remove_bads 1 -uniqueOnly 1 -SNP_pval 1e-6. Linkage disequilibrium (LD) was estimated as r2 values for all SNP pairs minimum 500 kbp apart with ngsLD v1.1.1 (Fox et al., 2019) and a LD decay curve was plotted using the script provided by (Fox et al., 2019) using 0.05% of all estimated r2 values. This indicated that a distance threshold of 2,000 bp was adequate for linkage pruning. A total of 560,735 SNPs were obtained using all samples (n=108, Table S1), and 319,236 SNPs were found to be in putatively neutral regions of the genome (see below). Admixture proportions were estimated by running NGSadmix version 32 (Skotte eKt al., 2013) for K equal 2 and 3 with 300 seed values, ensuring convergence (convergence was not reached for K = 4 and above). A principal component analysis (PCA) using the same SNP set was obtained with PCAngsd version 0.1 (Meisner & Albrechtsen, 2018).
The mitochondrial genome (mtDNA) for each individual was obtained as a consensus sequence of the reads mapped to the mitochondrial DNA sequence included in the reference genome. The individual sequences were generated in ANGSD (Korneliussen et al., 2014) using the option -doFasta 2 for sites with a minimum sequencing depth of 10X (-setMinDepth 10). An haplotype network was designed using mitochondrial SNPs with minor allele frequency >30% (total of 26 SNPs) in POPART (Leigh & Bryant, 2015) with the Median Joining Network algorithm (Bandelt at al., 1999).