Sampling took place in Upper Lake Constance, Germany, in spring and summer 2019 using three different methods: trawling in the pelagic zone, gillnetting in the littoral zone and electrofishing in three tributaries of Lake Constance (Fig. 1; Table 1). Details on the respective fishing methods can be found in Gugele et al. (2020). It is important to highlight that the trawling in the pelagic zone was performed along transects covering the entire area of Upper Lake Constance (Gugele et al. 2020) and that individuals for sequencing were randomly selected from individuals caught across the entire transect, thus representing the lake-wide pelagic diversity. A total of 95 sticklebacks were sampled. Fish from trawling and electrofishing were euthanised with clove oil (1 ml/l water, Euro OTC Pharma, Bönen, Germany) and total length was recorded. Each individual was photographed laterally using a digital camera (Pentax K3 II, 18–135 mm lens, fixed focal length). For genetic analysis, a piece of the caudal fin (approx. 0.5 cm²) was fixed in pure ethanol. The sex was determined by dissection and the fish were stored at -80 °C until further processing. Fish sampling was carried out according to local regulations (“Landesfischereiverordnung Baden-Württemberg”, LFischVO).

Figure 1.

Overview map showing the position of the individual sampling locations in Lake Constance and locations of outgroup populations (insert) used in the genetic analysis. Population information and geographic coordinates for all samples are provided in Suppl. material 2. Pelagic zone samples were sampled by trawling in transects across Upper Lake Constance and the location in the map is not highlighting the exact location (see Methods). Outgroup samples from the United Kingdom are freshwater and marine samples from the River Tyne, whereas outgroup samples from Norway and Germany are only freshwater samples (lake and stream).

To determine the number of lateral plates, formalin-fixed sticklebacks were stained with alizarin red according to a protocol modified from Peichel et al. (2001) (see Suppl. material 1: Methods for details). The mean number of lateral plates for both sides of each individual was used for further analysis (Fig. 2A). Individuals were classified into “fully plated” (> 30 lateral plates), “partially plated” (< 34 and > 10 lateral plates, with a gap) and “low plated” (< 10 lateral plates) according to Bell (2001). Differences in lateral plate number between sampling sites were tested using pairwise Steel-Dwass tests and differences in total length were tested using an ANOVA with a Tukey-Kramer HSD post hoc test in JMP 16.0.0 (SAS Institute, Cary, USA).

Figure 2.

Determination of morphometric and meristic traits in sticklebacks and results of the principal component analysis (PCA) using landmarks A (Top) Stickleback with stained bony structures (for more details, please refer to the text), allowing the determination of the lateral plate number. (Bottom) Position of 18 unique reference points (“landmarks”) on the body for shape analysis. (Right) Description of the individual locations of the landmarks B combined boxplot and violin plot to illustrate variation in size (left) and lateral plate number (right) of sticklebacks from different sampling sites. Lower case letters indicate statistically significant differences (size: ANOVA + Tukey-Kramer HSD post hoc test, lateral plate number: pairwise Steel-Dwass test). Box plots defined in the insert panel on the right C scatterplot showing the first two principal components (PCs), which explain most of the variance of the data (see axis labelling in percent). Sticklebacks were grouped according to sampling site (lit_zone = littoral zone, pel_zone = pelagic zone, Bru_inflow = Brunnisach inflow, Non_inflow = Nonnenbach inflow, Alt_inflow = Alter Rhein inflow) and sex (m = male, f = female). Wireframe graphs of the shape changes along the first two PCs in the PCA are shown on the right.

For the morphometric analysis, 18 unique reference points (‘landmarks’) were placed on digital images using TPSDig v.2.31 (Rohlf, Stony Brook University, New York, USA) (Fig. 2A). Statistical analysis was performed using “geomorph” v.4.0.3 (Baken et al. 2021). A Generalised Procrustes Analysis (GPA) was performed (function: “gpagen”) to remove differences in size, position and orientation. Possible errors in landmark placement were identified using the “plotOutliers” function and affected individuals were excluded from further analysis as appropriate. We tested for differences in shape between sites and sexes using a model that included fish size (i.e. allometric effects) as a covariate, sampling site and sex as a fixed factors, as well as the interaction of both fixed factors: coords = log(size) * site + sex + site:sex. A permutation-based Procrustes ANOVA using residual randomisation (function: “procD.lm” and “anova”, permutations = 9999; estimation method: ordinary least squares) was used to examine which factors have a statistically significant effect on shape. A pairwise post hoc test of sampling sites was performed with R package “RRPP” V.1.2.3 (Collyer and Adams 2018, 2021) (“pairwise” function), based on the previously developed model and a null model considering only size and sex (permutations = 9999, test type = distance between vectors (“dist”), confidence = 0.95). P-values were Bonferroni-Holm adjusted (Holm 1979). Procrustes shape coordinates for each individual were used for principal coordinate analysis (“gm.prcomp” function; (see Suppl. material 1: Methods for details).

Due to problems with storage of caudal fin tissue samples, brain tissue of the fish stored at -80 °C was dissected and used for DNA extraction. DNA for RAD sequencing was extracted from frozen brain tissue using the PureLink Genomic DNA Mini Kit (Invitrogen, Carlsbad, USA) according to the manufacturer’s instructions. Tissue was homogenised using 2.4 mm diameter metal beads (3 s, max. speed; Bead Ruptor 4, OMNI international, Kennesaw, USA), the homogenate digested for three hours at 55 °C (150 RPM) and RNA was removed using a 1h RNAse treatment. DNA was eluted in 50 μl elution buffer (NanoDrop 2000c, Thermo Fisher Scientific, Waltham, USA), sample concentrations checked using fluorometry (Qubit 3, Life Technologies, USA) and DNA quality determined using agarose gel electrophoresis. Non-degraded samples were shipped to Floragenexs (Eugene, Oregon, USA) for RAD library preparation and sequencing. Restriction-site Associated DNA Sequencing (RAD-seq) libraries were generated as detailed in Etter et al. (2011) using the restriction enzyme SbfI (New England Biolabs).

FASTQ data files were demultiplexed using GBSX v.1.3 (Herten et al. 2015) and reads were mapped to the G. aculeatus reference genome (release 92; Ensembl) using BWA v.0.7.17 (Li and Durbin 2010). Reads with a Phred quality score < 20 were removed before hard-calling genotypes using the Stacks v.2.6 gstacks module (Rochette et al. 2019). The Stacks output was subsequently filtered using VCFtools v.0.1.16 (Danecek et al. 2011) to remove indels (--remove-indels) and SNPs with minor allele frequencies below 0.05 (--maf 0.05), more than 25% missing data (--max-missing 0.75) and genotype quality below 10 (--minGQ 10). Only loci with a minimum read depth of 5 (--minDP 5), a minimum mean depth of 5 (--min-meanDP 5), a maximum depth of 50 (--maxDP 50) and a maximum mean depth of 50 (--min-meanDP 50) were kept.

Geographic structure across all samples was explored using PCA on filtered SNP data (n = 11,184 unlinked SNPs) in ade4 v.1.7-16 (Dray et al. 2007). SNP data (n = 28,194 SNPs) were pruned for linkage disequilibrium (LD) in Plink v.1.9 (Purcell et al. 2007) in 25 SNP windows with a five SNP window shift and an r² of 0.5. Chromosome 19 containing sex-determining regions was removed to avoid any sex-based bias (Peichel et al. 2004). Individual ancestry was assessed using Admixture v.1.3.0 (Alexander et al. 2009), using the LD-pruned SNP dataset, assessing between 1 and 7 clusters (K) to determine the optimal number of K using the tenfold cross-validation.

Genetic differentiation between sampling sites was investigated using haplotype-based relative allelic differentiation (F_ST') and absolute divergence (D_xy') between population pairs (populations module in Stacks) for all loci containing filtered SNPs. Haplotype-based estimates have the advantage of accommodating loci with more than two alleles contrary to SNP-based statistics (Bassham et al. 2018). Gene diversity (Hs), a haplotype-based equivalent to nucleotide diversity which is corrected for sampling bias originating from sampling small sample sizes, was calculated for each sampling site.

To identify genomic regions potentially under selection between littoral and pelagic sticklebacks, z-transformed F_ST' and D_xy' estimates for the pelagic-littoral population comparison were computed. Loci with z-transformed F_ST' ≥ 3 were classified as differentiated outlier loci. Loci with increased absolute divergence, zD_xy' ≥ 3, potentially highlight the differential sorting of ancient alleles between habitats. A z-transformed value ≥ 3 approximately corresponds to a p-value below 0.01. Outlier loci were further tested for signs of selection by comparing their interpopulation gene diversity differences (∆Hs = Hs_littoral – Hs_pelagic) to the genome-wide background. The expectation was that ∆Hs would be higher in loci under selection in the pelagic populations compared to the genomic background, driven by reduced Hs in the pelagic population. We compared median ∆Hs values for outliers (∆Hs_outlier) to the genomic background (∆Hs_bg) using a non-parametric two-sided Wilcoxon rank sum test and further compared the distributions of values using a Kolmogorov-Smirnov test. Genetic differentiation amongst the remaining populations and outlier loci overlaps amongst all pairwise population comparisons were also estimated.

Furthermore, we performed phylogenetic analyses of Lake Constance together with whole-genome data from outgroup populations from across Europe to confirm that Lake Constance stickleback cluster with the Baltic lineage (see Suppl. material 1: Methods for details) (Marques et al. 2019; Berner 2021).

Genetic association mapping using Genome-wide Efficient Mixed Model Association (GEMMA v.2.1; Zhou and Stephens (2012); Zhou et al. (2013)) was conducted to identify genetic association for body shape and mean lateral plate number. Associations with total length were not performed, as fish were sampled at different times of the year and, therefore, differed in age and size. We fitted Bayesian Sparse Linear Mixed Models (BSLMM) for each phenotype, providing the same genotype and relatedness matrix input files. As the BSLMM does not allow covariate files to correct for sex, we corrected sex-biased traits (body shape) using a linear model in R and used the residual body shape PC scores as input for the BSLMM analyses. We fitted five separate BSLMMs for each phenotype, averaged the results across chains and subsequently across runs. The BSLMM also estimates hyperparameters describing the genomic architecture of a trait, such as the proportion of variance in phenotypes explained by all SNPs (PVE), the proportion of variance explained by sparse effect loci (PGE) and the number of variants with major effects (n gamma). We estimated the means, median and 95^th confidence interval (CI) for these parameters. Furthermore, we identified SNPs as those with an average posterior inclusion probability (PIP) above 0.01 as “associated” (Comeault et al. 2014) and those with PIP above 0.1 as “strongly associated” with the studied phenotype (Chaves et al. 2016).

We tested if phenotype-associated SNPs were also significant outlier loci or showed increased genetic differentiation between littoral and pelagic sticklebacks, which would suggest potential selection acting on these phenotypes. To test if phenotype-associated SNPs showed increased genetic differentiation and divergence compared to a random genomic background, we performed random resampling of the same number of SNPs from the entire SNP dataset, estimated the mean F_ST' and D_xy' for the corresponding haplotype and repeated this 10,000 times to create a null distribution. Subsequently, we compared the means F_ST' and D_xy' of the phenotype-associated SNPs and the null distribution using a Wilcoxon test. We did this for each phenotype and for the sex-chromosome and autosomes separately.

Total length of sticklebacks differed between sampling sites (ANOVA: F_4,90 = 23.1534, p < 0.001; Fig. 2B), with pelagic individuals being smaller than all other groups (pelagic = 5.4 ± 0.3 mm, littoral = 6.4 ± 0.5 mm; inflow: Brunnisach = 6.3 ± 0.4 mm, Nonnenbach = 6.0 ± 0.5 mm, Alten Rhein = 5.9 ± 0.3 mm). Lateral plate numbers did not differ between sampling sites (Steel-Dwass test: p > 0.05, Fig. 2B), with 77.9% of sticklebacks being fully plated, 20.0% partially plated and 2.1% low plated.

The morphometric analysis, based on landmarks, showed a significant effect of size on body shape (ANOVA: p < 0.001; Suppl. material 1: table S1), indicating allometric effects. As these effects were unique to all sites (ANOVA: p = 0.019; Suppl. material 1: table S1), a size correction was not possible. Furthermore, both sex and sampling site had significant effects on shape (ANOVA: sex: p < 0.001, sampling site: p < 0.001; Suppl. material 1: table S1). Pairwise comparison of sampling sites revealed that fish from the littoral zone differed statistically significantly in shape from all other sites, regardless of size and sex (Table 2). The principal component analysis (PCA) generally visualised the results of the ANOVA, showing fish from the littoral zone being isolated from other sites in both sexes (Fig. 2C, Suppl. material 1: table S2). Changes in body shape along the first two principal components were most evident in the head region, the positioning of the pectoral fin and general body contour (Fig. 2C). The utilised broken stick model indicated that the first two principal components are statistically “meaningful”.

The admixture analysis suggested that sticklebacks in Lake Constance were highly admixed, as the genetic structure was best explained by 1 cluster (K = 1; Suppl. material 1: fig. S1). Low population structure was supported by the PCA (Fig. 3). PC1 and PC2 only explained negligible variation in genetic structure amongst individuals (1.75% and 1.52%, respectively). Furthermore, genetic differentiation was low between all sampling sites, both in relative divergence (mean F_ST' = 0.005 ± 0.002) and absolute divergence (mean D_xy' = 0.001 ± 2.61 × 10^-5) (Suppl. material 1: table S3). Differentiation between the pelagic and littoral three-spined sticklebacks was the lowest amongst all pairwise comparisons (F_ST' _LIT-PEL = 0.002). We did not detect any structure for the pelagic samples, which were sampled across Upper Lake Constance.

Figure 3.

Population structure. Principal Component Analysis representing individual structuring across Lake Constance populations. Analysis was performed on pruned data excluding the sex chromosome (11,184 SNPs). Colours indicate different sampling sites, while shapes represent habitats – inflow (circle), littoral (square), pelagic (triangle). Smaller, lighter data points show individual variation, while the larger shapes with a black centre indicate population Principal Component centroids, which were calculated as the mean of both the 1^st and 2^nd axes.

Furthermore, phylogenetic analyses showed that stickleback from Lake Constance clustered with the Baltic Lineage (Eastern European) stickleback (Suppl. material 1: figs S2, S3, table S4), supporting earlier findings that sticklebacks in Lake Constance are derived from the Baltic lineage (Berner 2021).

We detected 333 loci with zF_ST' ≥ 3, distributed across the entire genome (Fig. 4). Outlier loci showed on average increased absolute divergence (D_xy') compared to the genomic background, both on autosomes (Wilcoxon: W = 2993041, p < 0.001) and the sex chromosome (Wilcoxon: W = 6203, p = 0.012) (Fig. 4B). However, only 16 outlier loci showed strongly increased D_xy' values (zD_xy' ≥ 3). Comparisons between littoral and pelagic populations with inflow populations showed a similar picture (Suppl. material 1: figs S4, S5).

Figure 4.

Signatures of selection A Z-transformed haplotype-based F_ST' (zFst') estimates for loci (dots) across all chromosomes (noted on the x-axis). Outlier loci with zF_ST' ≥ 3 are shown in orange B absolute divergence (Dxy') between outlier loci (orange) and the genomic background on autosomes and the sex chromosome. Individual dots denote genomic loci and the distribution of values is shown by density plots. The sex chromosome was analysed separately due to lower recombination rates compared to autosomes and, therefore, potentially higher absolute divergence C comparison of delta gene diversity (∆Hs) between the genomic background (zFst' < 3; grey), outlier loci (orange) and outlier loci showing increased absolute divergence (zDxy' ≥ 3; red). Delta gene diversity was estimated by subtracting gene diversity in pelagic individuals from gene diversity in littoral individuals. Positive ∆Hs values are indicative of reduced gene diversity in littoral individuals and vice versa. Individual dots denote genomic loci. Box plots defined in Fig. 2B.

We further tested for signals of divergent selection by comparing gene diversity (Suppl. material 1: fig. S6) between littoral and pelagic populations. While the mean between-population difference in gene diversity of outlier loci (∆Hs_outlier) was not lower than the genomic background (∆Hs_bg), outlier loci showed more extreme values than the genomic background (Kolmogorov-Smirnov test: D = 0.1694, p < 0.001) (Fig. 4C). Outlier loci with increased absolute divergence had on average lower gene diversity in littoral sticklebacks (positive ∆Hs) (Fig. 4C), suggesting that these highly divergent loci are likely under selection in the littoral population.

Genome-wide association analyses for lateral plate number identified 41 associated SNPs (mean PIP > 0.01), with 7 SNPs (17.1%) showing very strong associations (mean PIP > 0.1). These were mainly located on chromosome 4, with one strongly-associated SNP on chromosome 2 (Fig. 4A). For body shape, the BSLMM detected 104 associated SNPs with PIP > 0.01 and one strongly-associated SNP on chromosome 21 with PIP > 0.1. Although body shape values were corrected for sex, a large proportion of associated SNPs (n = 34; 32.7%) were located on sex chromosome 19.

The proportion of variance explained by all loci was similar for lateral plate number (PVE_PN = 88.4%) and body shape (PVE_BS = 85.3%), but the proportion of PVE explained by sparse effect loci (PGEBS = 61.5%, PGE_PN = 77.8%) and the estimated number of sparse effect loci (mean n gamma_BS = 28; mean n gamma_PN = 6) were smaller for body shape compared to lateral plate number (Fig. 5B).

Figure 5.

Genome-wide association analyses (GWAS) for phenotypic traits A posterior Inclusion Probabilities (PIP) from BSLMMs for all SNPs (dots) across the genome are shown, with outliers SNPs passing significance threshold (PIP > 0.01) shown in red. Manhattan plots are shown for GWAS results with mean lateral plate number, body shape PC1-6 and total length B hyperparameters from BLSMMs are plotted, as the mean (large dot) and 95% confidence intervals (grey lines), for body shape (BS: yellow) and lateral plate number (PN: blue) C the distribution of genetic differentiation (F_ST‘) values for the permuted null distribution is shown as a histogram and mean differentiation for phenotype-associated loci is indicated as a red line. Results are shown for SNPs associated with later plate number (blue) and body shape (yellow), for autosomes and the sex chromosome separately. No trait-associated SNPs were detected for later plate number on the sex chromosome.

Body shape-associated loci were not strongly differentiated (i.e. outlier loci), but autosomal loci associated with body shape showed increased genetic differentiation between littoral and pelagic sticklebacks compared to the genomic background (Fig. 5C). Loci associated with lateral plate number and body shape-associated loci on the sex chromosome, did not show increased genetic differentiation compared to the genomic background (Fig. 5C).

Sticklebacks exhibit a high degree of phenotypic diversity between habitats, with some morphological traits having evolved in parallel during their postglacial dispersal into new freshwater habitats (McPhail 1993; Bell and Foster 1994). In particular, the reduction of lateral plate armour in freshwater populations is regularly observed (Hagen and Gilbertson 1972; Bell 2001). Although this process can happen very quickly (Bell et al. 2004), our findings show that sticklebacks in Lake Constance still exhibit largely full armouring. The littoral zone, where sticklebacks were found before the mass abundance in the open waters, is inhabited by several predatory fish species known to prey on sticklebacks (Donadi et al. 2017; Jacobson et al. 2019). Thus, full plating may have remained relevant as effective predator protection even for the pelagic population during the littoral breeding season (Reimchen 1994; Kitano et al. 2008; Rennison et al. 2019).

Littoral and pelagic sticklebacks differed slightly in snout length and body depth, with longer snouts and deeper bodies in littoral fish (Schluter and McPhail 1992; Arnegard et al. 2014). However, body shape was also clearly affected by size and sex of the fish. A pronounced sexual dimorphism is known for sticklebacks, with males having larger heads and mouths (Kitano et al. 2007). The underlying cause of the differences in body shape between habitats is unclear. A recent study was able to show that littoral and pelagic sticklebacks do not strongly differ in trophic position (Gugele et al. 2023). Habitat-specific morphological variation is common in sticklebacks (Gow et al. 2008; Webster et al. 2011), also in sympatric occurring limnetic and benthic populations (Nagel and Schluter 1998). Stable isotope ratios in muscle did not reveal differences in trophic position between sticklebacks from pelagic and littoral habitats in Lake Constance, but rather a mere preference to forage in the pelagic zone (Gugele et al. 2023). However, differences in δ13C levels in the liver, which were attributed to habitat shifts (Gugele et al. 2023), suggest that body shape differences are potentially related to differences in habitat usage rather than trophic eco-morphology.

Observed size differences between sticklebacks from the pelagic and littoral zone are likely related to differences in sampling time rather than growth rate, with fish from the pelagic zone having been captured in late March, while fish from the littoral zone were captured in June and July. It can, therefore, be assumed that this is a time-dependent increase in size over the course of the year. Future common garden experiments, temporal sampling throughout the year and more detailed phenotypic and trophic analyses could shed light on the eco-morphological basis of the rapid pelagic invasion of Lake Constance stickleback.

To date, it has been unclear whether pelagic stickleback in Lake Constance, which have increased rapidly in abundance since 2012 (Eckmann and Engesser 2019), originated within Lake Constance or are the result of a separate introduction. Our genetic results suggest a recent intralacustrine origin of pelagic stickleback in Lake Constance. Genome-wide patterns of genetic differentiation were overall weak between habitats, suggesting either a very recent expansion into the pelagic zone by sticklebacks and/or ongoing gene flow. A very recent origin and ongoing gene flow are supported by annual hydroacoustic surveys (Eckmann and Engesser 2019) and spatio-temporal sampling (Gugele et al. 2020), respectively. In Lake Constance, analyses of spatio-temporal movement of sticklebacks suggested temporal migration of sticklebacks from the pelagic zone to tributaries (‘inflow populations’) and back (Gugele et al. 2020), yet genetic analyses did not suggest a closer relationship between inflow and pelagic stickleback. The slightly weaker genetic divergence between littoral and pelagic sticklebacks, compared to inflow populations, suggests that pelagic sticklebacks might have originated from the littoral population and/or that gene flow is higher between the littoral and pelagic zone compared to tributaries (Gugele et al. 2020). These genetic differences are very subtle, though and not sufficient to confirm putative spawning locations of pelagic sticklebacks. Furthermore, the lack of population structure within the pelagic samples, which were sampled from across Upper Lake Constance, suggests that the pelagic population is genetically and spatially homogeneous, which might be expected under rapid expansion into the pelagic zone.

Stickleback from pelagic and littoral habitats were sampled during slightly different times in our study, potentially biasing estimates of genome-wide genetic differentiation between habitats. However, differences in sampling time would likely result in even lower estimates of genetic differentiation (compared to the ‘true’ differentiation), if pelagic sticklebacks, which were sampled in the spring, move into the littoral zone to spawn in the summer and are caught together with littoral sticklebacks. Thus, we believe that, overall, there is no strong genome-wide differentiation between habitats, in line with a recent expansion under ongoing gene flow.

The phylogenetic clustering of sticklebacks from Lake Constance as sister to Baltic stickleback from northern Germany (Suppl. material 1: fig. S2) is further in line with the theory that Lake Constance was historically colonised by individuals which shared a substantial proportion of ancestry with marine-like sticklebacks that were repeatedly introduced to streams and ponds in the Lake Constance system from catchments south to the Baltic Sea (Muckle 1972; Marques et al. 2019; Berner 2021; Hudson et al. 2021).

Despite the likely recent colonisation of the pelagic zone and minimal genome-wide differentiation between habitats, or lack thereof, we detected a polygenic signal of divergence with hundreds of outlier SNPs across the genome showing increased genetic differentiation between individuals from littoral and pelagic habitats. Such polygenic patterns of divergence between benthic and limnetic sticklebacks were also observed in Canadian populations (Härer et al. 2021). While many studies have identified individual large-effect loci associated with rapid local adaptation in different systems (Barrett et al. 2019; Thompson et al. 2020; Schluter et al. 2021), polygenic responses to selection, as observed in our and other studies (Laurentino et al. 2020; Salmón et al. 2021), can also lead to rapid local adaptation despite ongoing gene flow (Jain and Stephan 2017). Polygenic differentiation under gene flow and recent divergence could have two non-mutually exclusive explanations: a) rapid adaptation via sharing and re-assembly of adaptive alleles through selection on standing genetic variation (e.g. Terekhanova et al. 2014; Roesti et al. 2015; Fang et al. 2020); and b) strong divergent selection pressures arising from the differing environments in the pelagic versus littoral zone (Moser et al. 2016), which acts on adaptive phenotypes, such as body shape.

Genetic differentiation can occur without divergent selection, for example, through linked selection in low recombination regions or genetic drift due to population bottlenecks, yet these are unlikely explanations in this system. Firstly, linked selection is less likely to lead to increased differentiation over such short evolutionary timescales (Burri 2017) and patterns of genetic diversity between habitats do not indicate the presence of genetic bottlenecks, which is supported by large observed stickleback populations (Eckmann and Engesser 2019; Gugele et al. 2020). Furthermore, we observed increased absolute divergence in outlier loci, with strong divergence for a small subset of loci, suggesting a contribution of divergent selection in the genetic differentiation between habitats. Absolute divergence takes longer to build up than genetic differentiation, thus suggesting that increased divergence between pelagic and littoral sticklebacks is potentially due to the sorting of ancient adaptive alleles between habitats, which has been shown to play a role in freshwater adaptation in sticklebacks (Nelson and Cresko 2018). Differences in gene diversity at outlier loci between pelagic and littoral populations compared to the genomic background further support that genetic differentiation is potentially driven by divergent selection rather than variation in genetic diversity across the genome, which would be expected to lead to reduced diversity in both populations (Cruickshank and Hahn 2014; Burri 2017).

Overall, the fact that we observed genetic differentiation of many loci across the genome, despite low levels of genome-wide differentiation, indicates that habitat preferences might be, at least partially, genetically determined and not purely plastic, although a plastic component cannot be excluded.

Increased differentiation of body shape-associated autosomal loci between pelagic and littoral sticklebacks suggests that body shape, a trait that seems to differ between populations in these habitats, is under divergent selection between habitats. The observed signal is likely not due to chance, as loci associated with lateral plate number, a trait that does not differ between pelagic and littoral sticklebacks, do not show increased differentiation. We also recovered a well-studied lateral plate number associated genomic region on chromosome 4, further suggesting that we had sufficient power to detect large-effect loci. Hence, the genetic differentiation of body shape-associated loci suggests that the observed divergence between littoral and pelagic stickleback is not purely due to phenotypic plasticity, but is at least partly genetically determined. Variation in morphology could ultimately lead to assortative mating and divergence into distinct ecotypes over time (Garduno-Paz et al. 2020).

We did not test for genotype-association with variation in body size in our dataset, as individuals were sampled at slightly different time-points throughout the year. Whilst body size divergence between Lake Constance and stream sticklebacks has been demonstrated to be plastic and driven by differences in food availability (Moser et al. 2015), divergence in body size is a substantial driver of reproductive isolation in sticklebacks as it affects reproductive behaviour such as mate choice (Moser et al. 2015) and, thus, may facilitate speciation in future generations (Berner et al. 2017).

The authors have declared that no competing interests exist.

No ethical statement was reported.

This study was supported by the grant “SeeWandel: Life in Lake Constance - the past, present and future” within the framework of the Interreg V programme “Alpenrhein-Bodensee-Hochrhein (Germany/Austria/Switzerland/Liechtenstein)” whose funds are provided by the European Regional Development Fund as well as the Swiss Confederation and cantons. Arne Jacobs was supported by a Leverhulme Trust Early Career Fellowship (ECF-2020-509) and a University of Glasgow Lord Kelvin/Adam Smith Leadership Fellowship. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

A.B. conceived the project. Sampling was performed by A.R. and phenotypic analysis was performed by S.R. Genomic analysis was performed by C.D and A.J, with input from K.R.E. C.D. and S.R. wrote the first version of the manuscript with input from A.J. All authors contributed to the final version.

Carolin Dahms https://orcid.org/0000-0002-3283-7820

Samuel Roch https://orcid.org/0000-0003-3441-9830

Kathryn R. Elmer https://orcid.org/0000-0002-9219-7001

Albert Ros https://orcid.org/0000-0003-3241-9722

Alexander Brinker https://orcid.org/0000-0002-2433-56520

Arne Jacobs https://orcid.org/0000-0001-7635-5447

All newly-generated sequence data for this study are available in the NCBI SRA under the BioProject PRJNA1090479 with the following run accessions: SRR28409948–SRR28410042.