The repeated occurrence of similar phenotypes in independent lineages (i.e., parallel evolution) in response to similar ecological conditions can provide compelling insights into the process of adaptive evolution. An intriguing question is to what extent repeated phenotypic changes are underlain by repeated changes at the genomic level and whether patterns of genomic divergence differ with the geographic context in which populations evolve. Here, we combine genomic, morphological and ecological datasets to investigate the genomic signatures of divergence across populations of threespine stickleback (Gasterosteus aculeatus) that adapted to contrasting ecological niches (benthic or limnetic) in either sympatry or allopatry. We found that genome-wide differentiation (F_ST) was an order of magnitude higher and substantially more repeatable for sympatric benthic and limnetic specialists compared to allopatric populations with similar levels of ecological divergence. We identified genomic regions consistently differentiated between sympatric ecotypes that were also differentiated between or associated with benthic vs. limnetic niche in allopatric populations. These candidate regions were enriched on three chromosomes known to be involved in the benthic-limnetic divergence of threespine stickleback. Some candidate regions overlapped with QTL for body shape and trophic traits such as gill raker number, traits that strongly differ between benthic and limnetic ecotypes. In sum, our study shows that magnitude and repeatability of genomic signatures of ecological divergence in threespine stickleback highly depend on the geographic context. The identified candidate regions provide starting points to identify functionally important genes for the adaptation to benthic and limnetic niches.

We reanalyzed previously collected morphological, ecological and genomic data of wild-caught stickleback from British Columbia, Canada, obtained from several independent datasets. The genotype data for the genome-wide association (GWA) mapping of lake size and the F_ST analyses among small and large lakes were generated using ddRAD sequencing (Peterson, Weber, Kay, Fisher, & Hoekstra, 2012) for solitary stickleback populations from Vancouver Island. For 21 solitary populations, we obtained data on the proportion of benthic diet and gill raker numbers (a key trophic trait). Both of these datasets are from Bolnick and Ballare (2020); the genomic dataset was produced by Stuart et al. (2017). Genotype data used for F_ST analyses from three sympatric benthic-limnetic species pairs (Paxton, Priest and Little Quarry Lakes) were generated by Samuk et al. (2017) using the genotyping-by-sequencing method. Diet and gill raker number data for the sympatric species pairs from Paxton and Priest Lakes were obtained from Schluter and McPhail (1992). QTL data on different sympatric species pairs were obtained from several previous studies (Arnegard et al., 2014; Conte et al., 2015; Malek, Boughman, Dworkin, & Peichel, 2012).