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.