Discussion

Geographic context of divergence and parallelism

While sympatric benthic and limnetic species as well as allopatric solitary populations showed substantial variation in morphology and diet (Matthews, Marchinko, Bolnick, & Mazumder, 2010; Schluter & McPhail, 1992), the magnitude and patterns of genomic differentiation of benthic and limnetic ecotypes differed with geographic context (sympatry vs. allopatry). Sympatric benthic and limnetic species pairs show strong evidence of reproductive isolation (Rundle, Nagel, Wenrick Boughman, & Schluter, 2000). In line with these observations, we found that genomic differentiation between species was high in all three lakes (Fig. S3). Despite the reduced opportunity for gene flow afforded by allopatry, genomic differentiation between solitary populations from small and large lakes was an order of magnitude lower than that found for the sympatric species pairs. The strong divergence of benthic and limnetic pairs in sympatry suggests that these ecotypes may be further along the speciation continuum; perhaps these populations are moving towards a genome-wide phase of divergence, sometimes called ‘genome-wide congealing’ (Feder et al., 2014). In contrast, low levels of genomic divergence as observed for the solitary populations may indicate a more ‘genic phase’ of adaptation, with divergence limited to fewer loci (Feder et al., 2014). Linked selection (both genetic/genome hitchhiking and background selection) may also contribute to the observed higher degree of differentiation in the sympatric comparisons (Flaxman, Feder, & Nosil, 2013; Nosil, Funk, & Ortiz-Barrientos, 2009). These disparate patterns could be influenced by the process of reinforcement, where selection against hybrids increases reproductive isolation between lineages occurring in sympatry (e.g., Noor, 1999). Accordingly, a previous study on benthic and limnetic stickleback revealed female preference for mates of the same ecotype only in sympatric populations but not in allopatric ones (Rundle & Schluter, 1998).
Patterns of genome-wide differentiation were very consistent across sympatric species pairs (Fig. 2 and Fig. S3), and 23.3% of outlier regions were repeatedly differentiated. This is remarkable given that these species pairs have evolved independently in different lakes within the last 12,000 years (Taylor & McPhail, 1999). A study by Jones et al. surveying the same three lakes showed similar levels of repeatability, but was based on a lower number of genomic loci (Jones, Chan, et al., 2012). The repeated divergence of the same genomic regions for these sympatric species pairs suggests that a genetic constraint or bias may favor certain variants, loci, and/or regions during adaptation (reviewed in Bolnick et al., 2018). Given the young age of these species, standing genetic variation is expected to be more important for adaptation rather than de novo mutations. In stickleback, adaptation from standing genetic variation is thought to have been important for rapid adaptation to freshwater habitats (Colosimo et al., 2005; Jones, Chan, et al., 2012). The observed high degree of parallelism for the sympatric species pairs may indicate that adaptation from standing genetic variation has also been important for the repeated adaptation to benthic and limnetic niches (Jones, Chan, et al., 2012). Previous work on the sympatric benthic-limnetic pairs has also shown that divergence is heavily biased towards regions of the genome with suppressed recombination (Samuk et al., 2017); this bias may contribute to the similarity of genome-wide patterns, as only adaptive loci found within regions of low recombination may be able to diverge substantially.
The solitary populations and sympatric species pairs differ to a similar extent in trophic ecology and morphology (Fig. 1), yet the patterns of genomic parallelism are remarkably different. In contrast to the strong evidence of genome-wide parallelism found for the sympatric species pairs, there was no evidence that the same genomic regions were consistently differentiated (i.e., high FST) across solitary populations from different lakes exhibiting benthic-limnetic divergence (Fig. 2 and Fig. S2). Stickleback lake-stream pairs from Canada and Europe show similarly low levels of parallelism (Feulner et al., 2015; Rennison et al., 2019). The absence of parallelism for these solitary populations potentially means that standing genetic variation is not central to benthic-limnetic divergence in allopatry or that genetic variation may be more limited in some populations (Stuart et al., 2017), perhaps due to bottlenecks during colonization of upstream habitats. Also, it is possible that the marine ancestors were more heterogeneous than typically assumed (e.g., Stuart et al., 2017) and so different watersheds may have had genetically distinct founders. Population-specific effects of gene flow and genetic drift may also contribute to low levels of observed parallelism (Fitzpatrick, Torres-Dowdall, Reznick, Ghalambor, & Funk, 2014; Stuart et al., 2017). It would be interesting to have estimates of the demographic history of the stickleback populations investigated here, e.g., on the size of the founder populations. Genetic drift is stronger in small populations, which will reduce levels of parallelism across populations (MacPherson & Nuismer, 2017; Szendro, Franke, de Visser, & Krug, 2013). Further, if founder populations were small, there is also a high chance that the genetic diversity that selection could act upon was distinct across lakes, presumably limiting parallelism in genome-wide patterns of differentiation. Heterogeneity of the recombination landscape also doesn’t seem to play an important role when populations are diverging in allopatry (Samuk et al., 2017), which might contribute to different signatures of adaptation between the allopatric and sympatric comparisons.
Differences in the selective landscape among populations may also explain low levels of parallelism among solitary population pairs and between sympatric and allopatric comparisons. Theoretical work suggests that genomic parallelism decreases rapidly as the selection landscape becomes less parallel (Thompson, Osmond, & Schluter, 2019). In Trinidadian guppies, a famous example for parallel evolution, life-history traits and mortality were non-parallel between low and high predation, which could be attributed to stream-specific variation in disease and flooding (Fitzpatrick et al., 2014). Parallelism in morphological divergence of stream-lake pairs of stickleback decreases as environmental differences become non-parallel (Stuart et al., 2017). In marine-freshwater pairs of stickleback and in the ecotypes of the rough periwinkle, parallelism is also dependent on ecological similarity (Morales et al., 2019; Rennison et al., 2020). Even though all populations were adapting to similar benthic and limnetic niches, there may be cryptic environmental heterogeneity across the different lakes and within the discrete habitat categories that reduces overall parallelism in the selective landscape.

Identification of candidate regions for benthic-limnetic adaptation

Identification of repeatedly diverged genomic regions can increase our understanding of the predictability of evolutionary trajectories (Conte, Arnegard, Peichel, & Schluter, 2012; Stern & Orgogozo, 2008) and can help to detect important adaptive genes or traits as genetic drift is unlikely to result in patterns of repeated divergence (Schluter & Nagel, 1995). A handful of chromosomes appear to contribute disproportionately in the adaptation to benthic and limnetic niches. When looking at the genomic locations of candidate regions obtained from multiple analyses within this study, we detected evidence for significant clustering on chromosomes 4, 7 and 19. These chromosomes harbor multiple windows that were repeatedly differentiated across both sympatric benthic-limnetic species pairs and solitary populations from small and large lakes that differ substantially in their trophic ecology. Windows that were highly differentiated across sympatric benthic-limnetic species pairs and significantly associated with lake size across allopatric populations (double outliers) were also significantly enriched on chromosomes 4, 7 and 19. These chromosomes appear to be hot spots in threespine stickleback evolution (Conte et al., 2015; Jones, Chan, et al., 2012), and are enriched for multiple categories of ecologically relevant QTL (Peichel & Marques, 2017). Future comparative work examining the structure and evolutionary history of these chromosomes may provide insights into why these chromosomes are so central to the benthic-limnetic adaptation in this system.
To infer the functional relevance of genomic regions related to benthic-limnetic divergence, we incorporated published QTL data (Arnegard et al., 2014; Conte et al., 2015; Malek et al., 2012). We defined candidate regions as genomic windows that were either FST outliers across sympatric benthic-limnetic species pairs, significantly associated with lake size across allopatric populations or falling in both categories. The QTL that mapped to candidate regions were primarily located on chromosomes 4 and 7 (12 out of 32), again suggesting that these chromosomes are key to benthic-limnetic adaptation. Most QTL affected general body shape and five QTL affected the number of gill rakers (Table S2), a trait that has been shown to be important for the trophic ecology of these fish (Schluter, 1993). Benthic fish that feed on relatively large littoral invertebrates have fewer and shorter gill rakers than limnetic fish that feed on smaller pelagic zooplankton, both in sympatric and allopatric populations (Fig. 1; Bell & Foster, 1994; Schluter, 1993). Two more QTL affected other aspects of the trophic apparatus and feeding performance (Table S2). Note, however, that entire categories of potentially adaptive phenotypes have yet to be subjected to QTL analysis of benthic-limnetic divergence (e.g., immune traits, metabolic traits, the gut microbiota, etc.).
By integrating different types of data and looking for repeatability across sympatric and allopatric comparisons, we were able to identify 55 strong candidate regions underlying adaptation to the contrasting benthic and limnetic niches, several of which contain QTL for important morphological traits. These 55 candidate regions were FST outliers in benthic-limnetic species pairsand associated with lake size across solitary populations; hence, they are presumably targets of divergent selection. It is unlikely that the same regions were identified by chance as we used independent datasets, methods with different biases and had a fairly strict criterion for inclusion. This suggests that our candidate regions merit further molecular analyses to identify the genetic variants underlying adaptation to benthic and limnetic niches. Such investigation would reveal to what extent the parallelism in genomic regions observed here represents parallel changes in the same genes or perhaps even the same mutations. Parallelism is predicted to be more prevalent at higher levels of biological organization, as shown in the rough periwinkle (Ravinet et al., 2016) and Australian groundsel (Roda et al., 2013) were patterns of parallelism increased from SNPs to genomic regions to molecular pathways. Fine-mapping of genetic variants within these regions may provide good candidates for future gene editing studies in order to isolate their phenotypic and fitness effects.

Caveats

One limitation of our study is that we assessed parallelism at the level of genomic regions, which means that we cannot say whether the same loci or alleles underlie the observed patterns of differentiation. A strength of our approach is the integration of multiple datasets and types of data. However, a drawback of this approach is that we integrate sequence datasets generated using different methods (Genotyping-by-Sequencing vs. ddRAD sequencing). This means that some genomic regions were not surveyed in both datasets, so it is possible that we have underestimated the amount of parallelism between the sympatric and allopatric comparisons. Additionally, the use of arbitrary significance cutoffs means that we likely missed some regions that are evolving in parallel.

Conclusion

Our assessment of genomic parallelism in allopatry vs. sympatry revealed substantial differences in the magnitude of genomic divergence during adaptation to benthic and limnetic niches under different geographic contexts. Integration of independent datasets showed that parallelism is genome-wide among the three extant sympatric species pairs but very limited among solitary populations. By comparing candidate loci from the sympatric species pairs with those identified in solitary populations, we were able to identify candidate regions that might be essential for threespine sticklebacks’ ecological divergence along the benthic-limnetic axis.