Genomic divergence
To identify the extent of genetic divergence among sites, we used two approaches: 1) Principal Component Analysis (PCA) was used to explore genetic clustering of Mývatn stickleback, and 2) a model-based admixture analysis was used to determine population genetic structure by calculating the proportion of an individual’s genome that originates from different hypothetical ancestral gene pools (i.e., admixture coefficients). PCA was conducted in ADEGENET package , and as suggested by 100% of the initial PCs were retained when identifying the number of clusters. Admixture analyses were run using SNMF, which estimates admixture coefficients using non-negative matrix algorithms and makes no assumptions about drift or Hardy-Weinberg equilibrium (HWE). SNMF was run using the LEA package in R for number of ancestral populations (K) 1 – 10. Pairwise genome-wide Nei’s FSTs were calculated between all sites and habitats using VCFTOOLS. Contemporary effective population size (Ne ) was estimated using the LD-based method, as implemented in NeEstimator v2 . To remove physical linkage (as required in this method), we first thinned the SNPs across 10KB windows using VCFTOOLS. To provide confidence intervals of estimates of Ne , we calculated Ne across each chromosome independently.