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).