Results
Phenotypic and ecological variation along the
benthic-limnetic
axis
Across the 21 solitary stickleback populations, trophic ecology
(proportion of benthic food items; Fig. 1A) and morphology (gill raker
number; Fig. 1B) varied with lake size (log surface area; Table S1).
Stickleback from smaller lakes preferentially consumed benthic prey
(linear regression, R2 = 0.6117, P <
0.001) and had more gill rakers (R2 = 0.3039, P= 0.006). Sympatric benthic-limnetic species pairs from Paxton and
Priest Lake fell towards the ends of these distributions. Benthic
species mainly consumed benthic prey and had fewer gill rakers whereas
limnetic species consumed limnetic prey and had more gill rakers (Fig.
1). Including benthic-limnetic species pairs into linear regression
models substantially reduced the proportion of variation explained by
the models for diet (R2 = 0.3038, P = 0.003)
and gill raker number (R2 = 0.1415, P = 0.036).
This can be explained by the fact that the residuals of the sympatric
species pairs were consistently among the highest (Fig. S1), and
specialization appears independent of lake size.
Genomic architecture of benthic-limnetic differentiation in
sympatry and
allopatry
Mean genome-wide FST between allopatric solitary
populations from small and large lakes were much lower (0.026 - 0.059;
Fig. S2) than between benthic and limnetic species occurring in sympatry
(0.199 - 0.226; Fig. S3). These differences remained when
FST outliers were removed (solitary populations: 0.015 -
0.042; sympatric species pairs: 0.166 - 0.193). Genome-wide patterns of
differentiation were largely shared across sympatric benthic-limnetic
species pairs but not among solitary populations or between solitary
populations and benthic-limnetic species pairs; correlation coefficients
for FST values across the genome were low for the
small-large lake solitary population pairs (mean ρ = 0.023, ‘solitary’
in Fig. 2) and much higher for the benthic-limnetic species pairs (mean
ρ = 0.524, ‘benthic-limnetic’ in Fig. 2). When comparing genome-wide
patterns in FST between the small-large lake solitary
population pairs and the benthic-limnetic species pairs, correlation
coefficients were very low (mean ρ = –0.002, ‘solitary vs.
benthic-limnetic’ in Fig. 2).
Across benthic-limnetic species pairs, we identified 287
FST outlier windows based on mean normalized
FST values. These windows were unequally distributed,
with some chromosomes showing significant enrichment (P< 0.05, chromosomes 1, 4, 7 and 20) and others completely
lacking outliers (chromosomes 3, 6, 10, 14 and 15; Fig. 3A). Chromosome
7, which constitutes around 7% of the genome, had the highest
proportion of outliers (14%; Fig. S4A). For the sympatric
benthic-limnetic comparison, a total of 67 windows were repeatedly
differentiated (outlier in at least two of the three population pairs),
representing 23.3% of FST outlier windows. These
windows were spread across more than half of the chromosomes (12 out of
21; Fig. 3A and Fig. S3). In the allopatric solitary comparison for
small and large lakes, 520 windows were significantly repeatedly
differentiated (P < 0.05 after correction for multiple
testing), and these were spread across all chromosomes (Fig. S2). Seven
of these allopatric-divergence windows were also repeatedly
differentiated in the benthic-limnetic species pairs, which were located
on chromosomes 4,7, 18 and 19.
To identify candidate regions differentiated along the benthic-limnetic
axis across all 33 solitary populations, we also used genome-wide
association (GWA) mapping. This approach afforded more power by
leveraging the variance of more populations. Using this approach, we
found 1,015 of the 4,985 windows surveyed (20.4%) to be significantly
associated with lake size (our proxy for dietary niche; P< 0.05), although none of these windows remained significant
after very conservative Bonferroni correction. The proportion of
significant windows (20.4%) was much higher than would be expected when
assuming a 5% false positive rate and we were interested in the overlap
between datasets, which would hint at biological significance. Hence, we
proceeded in analyzing the candidate loci classified as significant
before Bonferroni correction. These candidate windows were distributed
across all 21 chromosomes (Fig. 3B), but significantly enriched
(P < 0.05) on chromosomes 11, 16 and 19 (Fig. 3B and
Fig. S4B).
A comparison of our candidates from the lake size GWA mapping in
solitary populations and FST outliers from the sympatric
benthic-limnetic species pairs led to the identification of genomic
regions likely to be important for adaptation to contrasting benthic and
limnetic niches in different geographic settings. There were 55
candidate windows shared between the two datasets (hereafter referred to
as ‘double outliers’), which represents 1.1% of the total 4,985 windows
(Fig. 4). Double outliers were enriched on chromosomes 4, 7 and 19 (Fig.
S4C). A chi-squared test of independence suggested that this overlap
between the two datasets was not more than expected by chance; candidate
windows significantly associated with lake size were not more likely to
be also FST outliers, χ2 (1, N =
4,985) = 0.336, P = 0.562. There was a weak, but significant,
negative correlation between normalized FST values and
association mapping -log10(P ) values across all
4,985 shared windows (Spearman’s rank correlation, ρ = -0.049, P< 0.001). Performing correlation tests for different
categories (non-significant, FST outliers, significantly
associated with lake size, double outliers; Fig. 4) showed a weak
negative correlation for non-significant windows (ρ = -0.079, P< 0.001). Spearman’s rank correlation coefficients were
positive (but not statistically significant) for the significantly
associated with lake size and double outlier categories (Fig. S5).
Characterization of candidate
regions
In total, 13 previously described QTL from different sympatric species
pairs (Arnegard et al., 2014; Conte et al., 2015; Malek et al., 2012)
overlapped with benthic-limnetic FST outlier windows.
These QTL containing FST outliers were mostly located on
chromosomes 1, 4 and 7 and included body shape and trophic traits such
as number of gill rakers (Table S2). An additional 16 QTL included
windows significantly associated with lake size, but not with sympatric
divergence. These QTL included body shape, dorsal spine length, number
of lateral plates and gill rakers and suction feeding index (Table S2).
Three QTL contained double outlier windows and these were associated
with color and shape of males as well as body shape (Table S2). Notably,
most QTL (10 out of 13) were located on chromosomes enriched for
FST outlier windows but no such pattern was observed for
significant lake size candidate windows (1 out of 16). For double
outlier, one out of three QTL mapped to chromosome 7, which was
significantly enriched (Table S2, Fig. 4C).