Genome‐phenotype‐environment associations identify signatures of selection in a panmictic population of threespine stickleback

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Genome‐phenotype‐environment associations identify signatures of selection in a panmictic population of threespine stickleback © 2023 the Authors Published version Strickland, K.; Räsänen, K.; Kristjánsson, B. K.; Phillips, J. S.; Einarsson, A.; Snorradóttir, R. G.; Bartrons, M.; Jónsson, Z. O. Strickland, K., Räsänen, K., Kristjánsson, B. K., Phillips, J. S., Einarsson, A., Snorradóttir, R. G., Bartrons, M., & Jónsson, Z. O. (2023). Genome‐phenotype‐environment associations identify signatures of selection in a panmictic population of threespine stickleback. Molecular Ecology, 32(7), 1708-1725. https://doi.org/10.1111/mec.16845 2023 1708 | Molecular Ecology. 2023;32:1708–1725.wileyonlinelibrary.com/journal/mec Received: 14 July 2022 | Revised: 1 December 2022 | Accepted: 13 December 2022 DOI: 10.1111/mec.16845 ORIGINAL ARTICLE Genomephenotypeenvironment associations identify signatures of selection in a panmictic population of threespine stickleback Kasha Strickland1,2 | Katja Räsänen3,4 | Bjarni Kristofer Kristjánsson2 | Joseph S. Phillips2,5 | Arni Einarsson6 | Ragna G. Snorradóttir2 | Mireia Bartrons7 | Zophonías Oddur Jónsson8 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2023 The Authors. Molecular Ecology published by John Wiley & Sons Ltd. 1Institute of Ecology and Evolution, School of Biological Sciences, University of Edinburgh, Edinburgh, UK 2Department of Aquaculture and Fish Biology, Hólar University, Sauðárkrókur, Iceland 3Department of Aquatic Ecology, EAWAG and Institute of Integrative Biology, ETH, Zurich, Switzerland 4Department of Biological and Environmental Science, University of Jyväskylä, Jyväskylä, Finland 5Department of Biology, Creighton University, Omaha, Nebraska, USA 6Mývatn Research Station, Mývatn, Iceland 7Aquatic Ecology Group, University of Vic (UVicUCC), Catalonia, Spain 8Faculty of Life and Environmental Sciences, University of Iceland, Reykjavík, Iceland Correspondence Kasha Strickland, Institute of Ecology and Evolution, School of Biological Sciences, University of Edinburgh, Edinburgh, UK. Email: [email protected] Funding information Icelandic Centre for Research, Grant/ Award Number: 195571052 Handling Editor: Sean Rogers Abstract Adaptive genetic divergence occurs when selection imposed by the environment causes the genomic component of the phenotype to differentiate. However, genomic signatures of natural selection are usually identified without information on which trait is responding to selection by which selective agent(s). Here, we integrate wholegenome sequencing with phenomics and measures of putative selective agents to assess the extent of adaptive divergence in threespine stickleback occupying the highly heterogeneous lake Mývatn, NE Iceland. We find negligible genome wide divergence, yet multiple traits (body size, gill raker structure and defence traits) were divergent along known ecological gradients (temperature, predatory bird densities and water depth). SNP based heritability of all measured traits was high (h2 = 0.42– 0.65), indicating adaptive potential for all traits. Environmentassociation analyses further identified thousands of loci putatively involved in selection, related to genes linked to, for instance, neuron development and protein phosphorylation. Finally, we found that loci linked to water depth were concurrently associated with pelvic spine length variation - supporting the conclusion that divergence in pelvic spine length occurred in the face of gene flow. Our results suggest that whilst there is substantial genetic variation in the traits measured, phenotypic divergence of Mývatn stickleback is mostly weakly associated with environmental gradients, potentially as a result of substantial gene flow. Our study illustrates the value of integrative studies that combine genomic assays of multivariate trait variation with landscape genomics. KEYWORDS adaptive divergence, environmental gradients, Gasterosteus aculeatus, gene flow, genome scans, landscape genomics | 1709 STRICKLAND et al. 1 | INTRODUCTION Elucidating the genetic basis of adaptive divergence in natural populations is an enduring goal of evolutionary biology (Stinchcombe & Hoekstra, 2008). Doing so can provide insight into evolutionary processes occurring in the wild, including the mechanisms associated with adaptive divergence, and the extent to which divergence takes place in the face of gene flow (Räsänen & Hendry, 2008; Rudman et al., 2018). Genetically, adaptive divergence is expected to manifest as blocks of differentiation across the genome, at regions containing genes that contribute to adaptation to divergent local environments (Nosil et al., 2009). Genome scan studies that test these expectations have identified genomic regions associated with adaptation to divergent ecological niches in numerous species (e.g., CampbellStaton et al., 2021; Marrano et al., 2018; Slate et al., 2002). This has been termed a “reverse ecology” approach, whereby loci associated with adaptation may be identified without measuring the traits themselves (Li et al., 2008). However, genome scan studies on wild populations are seldom able to provide precise information on which aspects of the phenotype selection is acting on, or which environmental factors are imposing selection (MacColl, 2011). A comprehensive view on the genomic mechanisms associated with adaptive divergence requires studies that combine phenotypic, environmental and genomic data. Accordingly, integrative approaches that combine association mapping with landscape genomics or selection scans to map genephenotypeenvironment associations could be a powerful means to infer the genomic basis of adaptation (Jones et al., 2013). Association mapping studies (e.g., genomewideassociations [GWA]; Santure & Garant, 2018) identify specific loci that underlie divergent traits, whereas landscape genomic studies can aid in determining loci associated with adaptive divergence, under the assumption that loci should be correlated with environmental variation that is directly or indirectly causing selection (Coop et al., 2010; Eckert et al., 2010). Combining association mapping with landscape genomics can strengthen the identification of genomic signatures of selection by allowing inference on whether causal variants of phenotypic variation are concurrently associated with environmental variation. This would be especially true in cases where correlations between phenotype and environment are mirrored in genetic polymorphisms, where at some quantitative trait loci, allele frequencies differ between groups that inhabit different environments. In the absence of dispersal barriers, many populations remain connected by gene flow during the process of adaptive divergence, often along environmental clines (Feder et al., 2012; Räsänen & Hendry, 2008). Gene flow is expected to constrain divergence, swamping locally adapted alleles and breaking up favourable allele combinations through recombination (Yeaman & Whitlock, 2011). Whilst in cases of substantial gene flow there may be little genomewide divergence, responses to natural selection may be present at specific genomic regions (islands of divergence; Wolf & Ellegren, 2017). Identifying genomic divergence in the presence of gene flow is a major challenge because most genome scan approaches require grouping individuals, which is not usually possible when individuals remain connected (de Villemereuil et al., 2014; Narum & Hess, 2011). Our perspective on adaptive divergence may therefore be biased towards studies where physical barriers to gene flow have facilitated divergence. Although such studies have provided great insight into evolutionary processes, studying processes of divergence in populations connected by gene flow can greatly improve our understanding of the relative roles of natural selection and gene flow in adaptive divergence (Richardson et al., 2014). Here, we employ GWA and landscape genomic approaches to map genephenotypeenvironment associations in threespine stickleback (Gasterosteus aculeatus) that inhabit Mývatn, a highly environmentally heterogeneous lake in NE Iceland. Threespine stickleback is a wellestablished model system in evolutionary biology (Hendry et al., 2013; Reid et al., 2021). Within freshwater systems, there is evidence for repeated adaptive divergence at both phenotypic and genomic levels (Härer et al., 2021; Hendry et al., 2009; Hudson et al., 2021), most commonly across the benthiclimnetic axis (e.g., Härer et al., 2021) but also across a range of other selective agents (e.g., predation; Reimchen, 1992, 2000; Reimchen & Nosil, 2002). However, most of the studies focus on simple environmental contrasts (e.g., benthic vs. limnetic or lake vs. stream), and only few studies have aimed to test intralacustrine divergence across environmental gradients. Mývatn is a large (37 km2) and geologically young lake, formed after a volcanic eruption ca. 2300 years ago. The lake is highly heterogenous, with temperature, water depth, invertebrate, and vertebrate (including stickleback) densities varying over space and time (Einarsson et al., 2004; Ives et al., 2008). Stickleback habitats in this lake can crudely be divided to five main types, across which stickleback vary phenotypically (Kotrschal et al., 2012; Kristjánsson et al., 2002; Millet et al., 2013). Previous work found that male stickleback had relatively larger brains in a “lava” (warm) than a “mud” (colder) habitat (Kotrschal et al., 2012), relatively longer spines in the north basin than the south basin (Millet et al., 2013), and divergence in gill raker morphology and diet among some of the habitats (Kristjánsson et al., 2002; Millet et al., 2013). Evidence for population genetic divergence of stickleback across the lake is mixed. Using samples collected between 1999 and 2002 (Ólafsdóttir et al., 2007) found evidence for genetic divergence using a suite of nuclear and mitochondrial markers between stickleback inhabiting the “lava” and “mud” habitats (microsatellites: FST = 0.08; mtDNA: FST = 0.223), suggesting the presence of two contrasting morphs. In contrast, using samples collected in 2009 and 12 nuclear microsatellite loci (seven of which were the same as in Ólafsdóttir et al., 2007; Millet et al., 2013) found little evidence for neutral genetic divergence of stickleback across five habitat types (average pairwise FST = 0.004), suggesting extensive gene flow. Given the known phenotypic divergence in traits typically under selection in stickleback, coupled with spatial variation in possible selective agents, our main goal here was to identify genomic signatures of selection in Mývatn stickleback occupying different environments. When information on fitness is not available, one common method to identify genomic signatures of selection is to identify as 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 1710 | STRICKLAND et al. genomic regions which are disproportionately divergent between groups compared to the rest of the genome (Hoban et al., 2016). We extended this definition to strengthen our identification of signatures of selection: we expected that genomic regions that bear a signature of selection should be both divergent across ecological axes, and contain loci associated with variation in divergent traits. We further measured single nucleotide polymorphism (SNP)- based additive genetic variation of divergent traits to gain insight into the evolutionary potential of traits that are spatially divergent. 2 | MATERIALS AND METHODS 2.1 | Study system and sampling Lake Mývatn is composed of two basins (North and South basin) that are connected by two narrow channels and vary in a range of abiotic and biotic conditions (Einarsson et al., 2004). The lake is spring fed, with geothermal hot springs (up to c. 23°C) feeding the northeast of the lake and coldwater springs (c. 5°C) feeding the southeastern parts. Most part of the lake follows the ambient temperature, which in summer is around 12– 13°C (Millet et al., 2013). The lake is shallow (1– 4 m), but with some deeper areas (up to c. 7 m) caused by historical diatomite mining in some parts of the North basin (Ólafsson, 1979). Productivity, as well as benthic, epibenthic and pelagic invertebrate abundance and community structure, also varies through space (Bartrons et al., 2015). Based on the combination of water temperature and depth, as well as vegetation and substrate, the habitats occupied by stickleback have previously been classified as warm, rocky shore, cladophorales, pondweed and mined (Millet et al., 2013, see below). In addition, longterm monitoring data shows that stickleback population density varies in both space and time (Einarsson et al., 2004; Phillips et al., 2023), with the North basin having higher densities than the South, and with periodically strong dispersal from the North to the South basin (Phillips et al., 2023). The stickleback population of Mývatn has been surveyed each year since 1991 as part of an ongoing longterm monitoring of population demographics (Phillips et al., 2023). This sampling is done during the third week of June and August each year by laying five unbaited minnow traps at predetermined locations over two 12 h periods (see Millet et al., 2013 for details). During monitoring, stickleback are counted to estimate catchperuniteffort (CPUE) and frozen for later analysis. For phenotyping and genotyping, a random subset of individuals (c. N = 100 per site for each of the day and night catches) have been stored since 2009. To study patterns of spatial divergence, we used stickleback from nine sites collected in June of 2012 due to the availability of detailed ecological data for this time point. 2.2 | Ecological data To characterize the environment that stickleback occupy, we focused on a set of ecological variables which represent putative selective agents. First, we used the same five habitat classifications (warm, mined, pondweed, cladophorales, and rocky shore) previously described in Millet et al. (2013; see Data S1 for details). We also collated data on ecological variables likely to reflect selective agents. These were: water temperature, water depth, stickleback CPUE, piscivorous bird density, and zooplankton abundances and community composition. These were chosen because temperature can affect metabolic processes, development, tolerance to parasite infections (Franke et al., 2017; Karvonen et al., 2013), as well as key life history traits (Kim et al., 2017; Mehlis & Bakker, 2014), whilst depth can affect sensory processes (Veen et al., 2017) invertebrate availability, and stickleback visibility to predators (Rypel et al., 2007). Stickleback CPUE was used as a measure of intraspecific competition (Bolnick, 2004), piscivorous bird density as a measure of predation pressure (Vamosi & Schluter, 2004), and invertebrate data as a measure of prey abundance and composition (Bolnick & Ballare, 2020). Note that although Mývatn stickleback are also subject to predation by brown trout (Salmo trutta) and Arctic charr (Salvelinus alpinus), available data on trout or charr abundances from Mývatn were not at the same spatial scale as the stickleback data (Phillips, Guðbergsson, & Ives, 2022). Hence, we did not integrate those data into our current analysis and note that fish predation is still a likely agent of selection. In general, it should be noted that our measures are only proxies for selection imposed by several correlated ecological factors. Temperature and water depth of each site were used as per Millet et al. (2013). Average temperature at each site was measured between 30 June 2011 and 18 August 2011 with a temperature logger (iButton Maxim Integrated Products), placed at middepth and recording at 3h intervals. CPUE for each site was estimated using count data from the longterm monitoring study from June 2012. To measure piscivorous bird density, we used data collected during the waterfowl census conducted each year at Mývatn, during which all waterfowl observed from predetermined vantage points with known survey areas are counted (Gardarsson, 1979). We used count data collected between 15 May and 10 June 2012 on the following species known to predate on stickleback: horned grebe (Podiceps auritus), redbreasted merganser (Mergus serrator), great northern diver (Gavia immer), red throated diver (Gavia stellata) and goosander (Mergus merganser). Note that the Arctic tern (Sterna paradisaea) is abundant at Mývatn, and predates on stickleback, but this species is not counted during the bird census. We calculated the density of piscivorous birds (summed across all taxa) in each surveyed segment of the lake (number/m2; Figure S1). We used invertebrate data from Bartrons et al. (2015), which were collected by conducting surveys of the epibenthic and zooplanktonic community. Crustaceans (including Daphnia, copepods and epibenthic cladocerans) as well as rotifers are important food sources for stickleback in Mývatn (e.g., Kristjánsson et al., 2002). Although chironomid larvae are a main food source for Mývatn stickleback, data on midge larval abundance were not of sufficient spatial resolution to be used (see Bartrons et al., 2015). However, the benthic community is spatially correlated with the epibenthic and zooplankton 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 1711 STRICKLAND et al. community in the South basin (Bartrons et al., 2015), suggesting that measures of pelagic and epibenthic zooplankton may serve as a proxy measure for the benthic community in Mývatn. Briefly, three transects were conducted between June– July 2012, during which integrated vertical tows of the whole water column were made at each of 31 sites, spaced 500– 600 m apart (see Figure S2 for distribution, and Bartrons et al., 2015 for more details on sampling and sample processing). Each pooled sample of 15 L was filtered through 63μm mesh and counted in entirety under a binocular microscope. We used data from the second transect (25 July 2012) as the spatial resolution in this transect was the greatest. We used data collected from the closest site to each stickleback sampling site (distance to closest zooplankton site: min = 290 m, max = 1365 m). All stickleback sites were within 2.55 km of the nearest site used to collect invertebrate data, which was the distance at which zooplankton communities were found to be spatially autocorrelated. We used number per litre (n/L) of each taxon at the sites closest to stickleback sites. To summarize variance between sites for use in downstream analyses, we ran a principal components analysis (PCA) using the native stats package in R version 4.1.2 (R Core Team, 2021). This summarized the invertebrate data in four general axes, described in detail in Table S1. PC1 (zPC1) described the overall abundance of crustaceans and rotifers, and explained 29% of the variation; PC2 (zPC2) described the negative covariance of rotifer sp. and Alona sp. with planktonic and epibenthic crustaceans, which explained 17% of the variation; PC3 (zPC3) described the negative covariance of the rotifer Keratella and the cladocerans Acroperus harpae and Chydorus sphaericus with Daphnia longispina, and explained 15% of the variation; and PC4 (zPC4) described the negative covariance of the cladocerans Eurycercus lamellatus and Macrothrix hirsuticornis with Daphnia longispina, Cyclops abyssorum and Asplancha, explaining 11% of the variation in the data. Overall invertebrate abundance (described in zPC1), was highly negatively correlated with stickleback CPUE. We therefore used only CPUE and not zPC1 for downstream analyses. 2.3 | Phenotypic data We randomly selected individuals from each habitat type for phenotyping and genome sequencing from the frozen subset of the longterm monitoring samples. Individuals were thawed, weighed on an electronic balance (wet mass, nearest mg) and their total length (TL) measured using a ruler (to the nearest mm). The right pectoral fin was then cut and stored in 96% ethanol for DNA isolation. We measured traits typically under selection in stickleback: body size, defence traits (armour plate number and length of spines) and dietary traits (gill raker morphology and gut length; Härer et al., 2021; Hendry et al., 2009; Reid et al., 2021). Specifically, for each individual we measured the following 10 traits: total length (TL), total gut length (gut length), number of lateral armour plates (plate number), length of the first dorsal spine (DS1), length of the second dorsal spine (DS2), length of the pelvic spine (PS), length of the second gill raker on the first gill arch (GRL2), length of third gill raker on the first gill arch (GRL3), gap width between second and third gill rakers (GRW), and number of long gill rakers on the first gill arch (GRN) (see below). Note that we measured the second and third gill rakers, rather than the first (which is usually used in studies of stickleback trophic phenotype), because in some cases gill arches broke during dissection. After measurement of TL, each individual was dissected to remove the stomach and the gut, and any tapeworm (Schistocephalus solidus) parasites. Gut length (from the sphincter at the end of the oesophagus to the end of the digestive tract) was measured (to the nearest mm) using a ruler. To aid morphological measurements, fish were stained with alizarine red using standard protocols (Millet et al., 2013). Fish were bleached using a 1:1 ratio of 3% H2O2 and 1% KOH and then stained in a solution of alizarin red and 1% KOH (Bell, 1982). After staining, digital images were taken of the left side of the fish with a Canon EOS 600D digital camera, with mm paper for scale. From these images, plate number was counted and the length of the spines (DS1, DS2 and PS) measured to the nearest hundredth of a millimetre. After imaging, we dissected the first gill arch and, where necessary, restained before mounting it between two glass plates and photographing using a digital camera (Nikon Coolpix 4500) mounted to a stereomicroscope (Leica MZ12) with mm paper for scale. We used the digital images of gill arches to measure GRL2, GRL3 and GRW (in mm) and counted GRN. All measurements were taken from the digital images were done using segmented tool in ImageJ (Schneider et al., 2012). 2.4 | Whole genome resequencing and SNP detection Genomic DNA was isolated and purified from the ethanol stored fin clips using MachereyNagel nucleomag tissue kit, following the manufacturer's protocol. Paired end, PCRfree 150bp insert libraries were then prepared for whole genome sequencing using the DNBSeqTM platform by BGIHongkong, generating an average of 40 million cleaned reads per individual (min = 33.4 million, max = 41.8 million, equating to an average depth of coverage of 10×). The clean pairedend reads were aligned to the threespine stickleback genome assembly version 5 (Nath et al., 2021) with BOWTIE2 (version 2.4.1) using default parameter settings (Langmead & Salzberg, 2012), and sorted and indexed using SAMTOOLS (version 1.10) (Li et al., 2009). Variants were discovered using the short variant discovery pipeline of GATK (van der Auwera et al., 2013). SNPs and indels were quality filtered according to GATKs best practices guidelines, using the following hard filters: QualByDepth > 2.0, FisherStrand bias < 60.0, RMSMappingQuality < 40.0, Mappin gQualityRankSumTest < −12.5, ReadPosRankSumTest < −8.0, StrandOddsRatio > 3.0, variant quality score < 30.0. For all following analyses, we removed mitochondrial variants, indels and multiallelic variants, as well as variants identified on either of the sex chromosomes (Peichel et al., 2020). We then filtered the remaining autosomal SNPs for genotype depth less than six or greater than 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 1712 | STRICKLAND et al. 100; minor allele counts less than four; missingness <20%. The sex of individuals was confirmed using the proportion of reads with depth greater than eight mapped to the X versus Y chromosome (Peichel et al., 2020). 2.5 | Statistical analyses We conducted a series of analyses to test for (1) phenotypic and genomic divergence in relation to geographic location and ecological variables, (2) genomic architecture of traits, and (3) associations of genomic variation with environment. 2.5.1 | Phenotypic divergence Phenotypic divergence was analysed using multivariate and univariate Bayesian linear mixed models. All traits were standardized to have a mean of zero and standard deviation of one to improve comparability and model convergence. All models described below were fit with a Gaussian error distribution in the MCMCglmm package in R (Hadfield, 2010). Fixed effects were given weakly informative flat priors, random effects were given default inverse Wishart priors, and we fit full unstructured covariance matrices for the multivariate models. Each model was run for a total of 1,020,000 iterations with a burnin of 20,000 iterations and thinning of 1000 iterations, which resulted in low autocorrelation. Convergence of models was assessed by examining traceplots to visualize sampling mixing and by assessing effective sample sizes and autocorrelation. Effects of site and habitat To investigate the extent of phenotypic variation among sites and habitats, we compared three linear mixed effects models with different random effects structures using the deviance information criterion (DIC) to identify the model with the most support, at ΔDIC = 2 (Spiegelhalter et al., 2014): a model with no random effects (i.e., a “null” model), a model that included site as a random effect, and a model that included habitat type as a random effect. These models were designed to identify whether there was phenotypic variance between all habitats or all sites. We grouped traits in the response variable to investigate phenotypic divergence in functionally correlated traits. TL and gut length were each fit as a univariate response; the four defence traits (plate number, DS1, DS2 and PS) and the four gill raker traits (GRL2, GRL3, GRW, GRN) were fit as multivariate response traits, respectively. This resulted in a total of 12 models (three models with different random effects structures for each of the four responses), all fitted with sex as a fixed factor. Models for gut length, defence traits, and gill raker traits included TL as a covariate. Fitting the interaction between sex and TL did not improve model fits, so we present results from models with sex and TL fit as single term effects. Owing to varying levels of replication for each sex at each site (Table 1) we cannot estimate or exclude sex differences in the patterns of divergence. TABLE 1 Overview of sample size, ecological variables, and population genetic parameters for threespine stickleback sampled across 12 sites and five habitat types in Lake Mývatn Site Nind m f Habitat type Temp (°C) Depth (m) CPUE zPC2 zPC3 zPC4 Bird density (n/m2)HoHeF 23 13 67Cladophorales 12.2 3.24 34 0.689 0.507 −0.492 0.0004 0.262 (0.007) 0.268 (0.027) 0.021 (0.027) 27 97 2 Cladophorales 11.9 2.52 68 −1.026 0.231 −0.898 0.0004 0.266 (0.005) 0.268 (0.018) 0.008 (0.018) 44 15 15 0Cladophorales 12.1 2.95 214 1.074 −1.764 −0.053 0.0007 0.319 (0.043) 0.268 (0.061) −0.191 (0.061) 124 37 31 6Mined 12.6 4.75 885 −0.708 −0.27 −0.069 0.0006 0.298 (0.027) 0.268 (0.001) −0.111 (0.001) 128 20 20 0Pondweed 13 1.62 825 −0.559 0.256 0.007 0.0006 0.280 (0.021) 0.268 (0.078) −0.041 (0.078) 135 8 3 5 Pondweed 11.6 214 1.793 −1.324 −1.059 0.0005 0.261 (0.006) 0.268 (0.023) 0.026 (0.023) CS 40 832 Rocky shore 12.7 1.3 177 −0.546 1.189 −1.1 0.0003 0.281 (0.022) 0.268 (0.082) −0.041 (0.082) DN 10 8 2 Pondweed 13 1.18 920 −0.993 −0.078 −0.569 0.0006 0.286 (0.016) 0.268 (0.059) −0.061 (0.059) HS2 33 20 13 Warm 23.4 0.15 1342 −0.559 0.256 0.007 0.0018 0.292 (0.025) 0.268 (0.093) −0.091 (0.093) Note: The columns indicate site, total sample size (Nind), number of males (m) and females (f) sampled per site, habitat type, average temperature (Temp, °C), water depth (depth, metres), catch per unit effort (CPUE), invertebrate principal component scores zPC2, zPC3, and zPC4, density of piscivorous birds (Bird density) and expected heterozygosity (He), observed heterozygosity (Ho) and inbreeding coefficient (F). Population genetic parameters were measured using vcftools and averaged across all individual, with standard deviations presented in parentheses. 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 1713 STRICKLAND et al. The multivariate models fit full variance– covariance matrices for random effects, allowing us to estimate phenotypic covariances at both individual (residual) and habitat/site levels (Chenoweth et al., 2010). Note that because of the different sample sizes for each group of traits (see below), it was not possible to run a single multivariate model with all traits to measure the full covariance matrix. Fish from sites were considered phenotypically divergent if 95% credible intervals of the posterior distributions of predicted trait values did not overlap. Association with ecological parameters We next ran a suite of univariate, linear mixed effects models that investigated the effect of ecological predictors on each phenotype independently. All models were fitted with TL and sex as fixed effects and site as a random effect (except for model on total length, which only had sex as a fixed factor). Because many of our eight ecological variables were highly correlated (Table 1), we ran one model per ecological predictor per phenotype (total = 7 models per phenotypic trait) and compared each to a null model without any ecological variables using DIC (as above). In these models, all predictor variables were standardized to have a mean of 0 and standard deviation of 1 to improve comparability between models. We present the mode and 95% credible intervals of the posterior distribution for the linear coefficients for each ecological predictor, unless ΔDIC to the null model was >2 because this indicated that the null model was not improved by fitting the ecological variable. 2.5.2 | Genomic divergence To identify the extent of genetic divergence among sites, we used two approaches: (1) Principal component analysis (PCA) was used to explore genetic clustering of Mývatn stickleback, and (2) a modelbased admixture analysis was used to determine population genetic structure by calculating the proportion of an individual's genome that originates from different hypothetical ancestral gene pools (i.e., admixture coefficients). PCA was conducted in ADEGENET package (Jombart & Ahmed, 2011), and as suggested by Jombart et al. (2010) 100% of the initial PCs were retained when identifying the number of clusters. Admixture analyses were run using SNMF, which estimates admixture coefficients using nonnegative matrix algorithms (Frichot et al., 2014) and makes no assumptions about drift or Hardy– Weinberg equilibrium (HWE). SNMF was run using the LEA package in R (Frichot & François, 2015) for number of ancestral populations (K) 1– 10. Pairwise genomewide Nei's FSTs were calculated between all sites and habitats using VCFTOOLS. Contemporary effective population size (Ne) was estimated using the LDbased method, as implemented in NeEstimator version 2 (Do et al., 2014). To remove physical linkage (as required in this method), we first thinned the SNPs across 10 kb windows using VCFTOOLS. To provide confidence intervals of estimates of Ne, we calculated Ne across each chromosome independently. 2.5.3 | Genomewideassociation analyses To estimate SNPbased heritability for each trait, we ran Bayesian mixture models for each trait independently using the BAYESR software (Moser et al., 2015). We selected this method as it has been found to be more accurate in estimating additive genetic variance explained by SNPs than alternative methods (e.g., BSLMM or LMM; Moser et al., 2015). All models included TL and sex as covariates (except when modelling TL explicitly, which included only sex as a factor), and all traits were standardized to a mean of 0 and standard deviation of 1. We also aimed to estimate pairwise genetic correlations between traits using bivariate models in BAYESR. However, none of these models converged, probably owing to the large sample sizes usually required to estimate genetic covariances, and hence are not reported further. To identify SNPs underlying phenotypic variation, we performed a GWA using a linear mixed effects model approach with the program gemma (Zhou & Stephens, 2012). This program fits models that control for relatedness among samples and/or population stratification, which reduced false positives even with small sample sizes. gemma was selected because it implements generalized linear mixed models for traits with nonnormal distributions, and therefore robustly handles data on different scales. These models estimate the linear coefficient for the relationship between each SNP in turn and a given trait. When identifying whether a SNP was putatively associated with trait variance, Wald tests were used to assess significance (pvalue cutoff of - log10(P) = 5). This value was chosen to minimize the number of expected false positives given the number of loci in our data set (12 false positives are expected). All models included sex as a fixed factor and length as a covariate, except for the model with total length as a response variable which just fix sex. The number of SNPs used in each model was dependent on the number of individuals available for each trait, because gemma only models SNPs with <5% missingness. We then identified whether there was overlap between regions that were associated with trait variation in our data, and quantitative trait loci (QTLs) previously mapped on the stickleback genome (Peichel & Marques, 2017), using LIFTOVER and custom R scripts. 2.5.4 | Genomeenvironment association analyses We investigated gene– environment associations using latentfactor mixed models (LFMM) in the R package LEA (Frichot & François, 2015), which models loci against environmental variables while controlling for unobserved latent variables (i.e., spatial structure or autocorrelation). Models were run for 100,000 iterations, with 10,000 burnin cycles and five replicates. zscores were combined from five runs, and we used adjusted pvalues using the genomic control method (a recalibration procedure which decreases the false discovery rate (Frichot & François, 2015)). Following the guidance of Frichot and François (2015), the resulting pvalues for the outlier tests were adjusted for multiple testing 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 1714 | STRICKLAND et al. to Qvalues using the falsediscovery rate method (Benjamini & Hochberg, 1995) as implemented in the “qvalue” 2.6.0 package. A Qvalue for SNPs of <0.05 was considered significant. We compared DIC of models that included one ecological variable at the time to try to identify which ecological variable best predicted the genomic data. To identify whether gene– environment associations were linked to phenotypic variation, we then investigated whether any of the candidate SNPs identified in any of the LFMM analyses were linked (within 5 kb) to SNPs identified in GWA analyses. We used 5 kb windows as a proxy for linkage on the genome, which is a common window size used when physical linkage across the genome is unknown (Artemov et al., 2017, Kingman et al., 2021). To identify the extent to which we would expect an overlap between SNPs identified in LFMM and GWA by chance, we ran permutation tests using custom R scripts. These kept the number of SNPs identified in each analysis and the genomic location of the GWA SNPs constant but randomized the genomic location of environmentally associated SNPs. This was done by randomly selecting SNPs from the global set of loci (N = 1,205,604 SNPs). We ran 1000 permutations and the probability of an overlap occurring by chance was calculated as the proportion of times an overlap occurred in the permutations. To explore the molecular function of genomic regions that showed signatures of selection, we first analysed candidate genes for enrichment of molecular functions. To do this, we identified genes that the candidate SNPs were within 5 kb of and compared these candidate genes with the reference set of 20,805 genes across the stickleback genome (“gene universe”). Gene ontology (GO) information was obtained from the stickleback reference genome on ENSEMBL using the R package BIOMART (Durinck et al., 2009), and functional enrichment was investigated using the package TOPGO 2.42 (Alexa & Rahnenfuhrer, 2020) and the Fisher's exact test (at p < .01). To reduce false positives, we pruned the GO hierarchy by requiring that each GO term had at least 10 annotated genes in our reference list (“nodeSize = 10”). Second, physical overlap on the genome between candidate SNPs identified in environmentassociation analyses and GWA SNPs (as identified above) indicate regions of the genome under selection. For regions of the genome containing both environmentally associated candidate SNPs and GWA SNPs, we therefore identified (1) the genes in this region (within 5 kb), (2) whether haplotypes on that gene in our data set are predicted to cause variation in protein translation, and (3) the function of the gene. We used the program “snpEff” (Cingolani et al., 2012) to detect whether SNPs fall on coding regions changing amino acid sequence. 3 | RESULTS Out of 200 individuals originally sent for sequencing, 14 samples did not pass the quality control and were therefore not sequenced, resulting in a total of 186 sequenced individuals. Quality filtering after variant discovery resulted in a data set of 1,205,604 SNPs. This resulted in an average of one SNP per 270 bp across the stickleback genome. During dissections, some samples were broken resulting in slightly different sample sizes for each trait; total length: N = 186, gut length: N = 106, plate number: N = 160, DS1 and DS2: N = 159, PS: N = 158, GRN: N = 133, GRW: N = 161, GRL2 and GRL3: N = 158 and 159. 3.1 | Phenotypic divergence 3.1.1 | Total length Model comparisons showed that TL varied across sites rather than habitats (Table 2). This pattern was predominantly driven by stickleback from HS2 site being shorter than stickleback from other sites (Figure 1). TL was negatively correlated with both temperature (Table 3) and bird density (Table 3), but the model with temperature as a predictor fit the data better (i.e., had a lower DIC; Table 3). All traits, except plate number and GRN, were positively correlated with TL (Table 2), and all results presented hereafter refer to effects on sizecorrected traits. Defence traits There were no sex differences in the relative length of either dorsal or pelvic spines, but males had more armour plates than females (Table 2). Model comparisons suggested that defence traits tended to vary according to habitat rather than site, although there was only very weak statistical support for this effect (Table 2, Figure 1). Whilst none of the ecological variables predicted relative length of DS1 and armour plate number (Table 3), relative length of DS2 increased as the density of piscivorous birds increased (Table 3) and relative length of PS increased in deeper water (Table 3, Figure 2). However, statistical support for these environmental associations was weak (the ΔDIC to the null model was within 2, and the lower credible interval of the posterior distribution of the linear coefficient was only just above zero; Table 3). Individuals with relatively longer DS1 had correspondingly longer DS2 and PS (pairwise phenotypic covariances [CoV], posterior mode and 95%CI: DS1:DS2 = 0.224[0.165, 0.289]; DS1:PS = 0.133[0.077, 0.190]; DS2:PS = 0.122[0.070, 0.190]), indicating that spine traits covaried at the individual level. However, relative length of spines did not seem to covary with armour plate number (see Table S3). There was no evidence for phenotypic covariance across the habitats between any of the defence traits, suggesting that spatial divergence in defence traits was not correlated at the habitat level (see Table S3). However, to increase statistical power for detecting divergence in phenotypic covariances, greater level of replication is likely needed. Trophic traits Males had relatively longer GRL2 and GRL3, and relatively more gill rakers (GRN), than females (Table 2), but there were no sex differences in GRW (Table 2). Model comparisons showed that gill raker traits varied 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 1715 STRICKLAND et al. TABLE 2 Results of linear mixed effects models used to investigate the extent of spatial divergence of threespine stickleback phenotypes across Lake Mývatn Trait ΔDIC βV Null Site Habitat Length Sex M Site Habitat Residual Length 68.4 08.05 – −0.512 (−0.782, 0.279) 0.495 (0.130, 1.023) – 0.562 (0.451, 0.612) Defence traits Plate number 02.21 1.68 0.012(−0.010, 0.031) 0.339(0.009, 0.762) – 3.919(0.094, 12.241) 0.988(0.766, 1.195) DS1 0.089(0.077, 0.100) 0.186(−0.039, 0.387) – 4.634(0.076, 9.501) 0.309(0.249, 0.388) DS2 0.091(0.081, 0.102) 0.076(−0.136, 0.297) – 4.295(0.059, 11.727) 0.299(0.231, 0.369) PS 0.076(0.064, 0.088) −0.046(−0.271, 0.159) – 3.281(0.080, 9.915) 0.337(0.267, 0.410) Trophic traits GRL2 41.26 02.8 0.069(0.050, 0.089) 0.543(0.222, 0.870) 0.519(0.092, 1.333) – 0.551(0.418, 0.700) GRL3 0.074(0.056, 0.091) 0.559(0.259, 0.882) 0.564(0.089, 1.521) – 0.492(0.372, 0.616) GRW 0.075(0.056, 0.095) 0.088(−0.265, 0.424) 0.430(0.063,1.090) – 0.620(0.478, 0.792) GRN 0.011(−0.013, 0.033) 0.449(0.051, 0.843) 0.945(0.127, 2.318) – 0.772(0.592, 0.982) Gut length 0.39 04.42 0.072(0.057, 0.085) −0.263(−0.497, −0.042) 0.274(0.059, 0.655) – 0.210(0.151, 0.271) Note: Models with different random effects structure (no random effect [“Null”], site [“Site”] or habitat type [“Habitat”]) were compared using DIC, and results from the model with the lowest DIC (at ΔDIC = 2) are shown. (For details of analyses and model structure, see material and methods). The results are shown for the following traits: total length (Length), number of lateral plates (Plate number), lengths of two dorsal spines (DS1 and DS2), the pelvic spine (PS) and two gill rakers (GRL2 and GRL3), as well as gill raker gap width (GRW), number of gill rakers (GRN) and gut length. β reflects the linear coefficient estimate; V is variance estimates of random effects. Posterior means for all parameter estimates are presented, with 95% credible intervals in subscript parentheses. The models included sex and, where appropriate, length as a covariate and therefore all results reflect effect on size corrected traits. 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 1722 | STRICKLAND et al. which activates mitosis and regulates the dynamics of the cell cycle. Whilst these findings suggest candidate molecular functions underlying responses to natural selection, explicit follow up studies (e.g., CRISPR or gene expression studies) would be needed to get at the causal relationships with pelvic spine length variation. Many of the regions of the genome that were correlated with environmental variation in our study were not linked to observed trait divergence. This is not surprising given that we only measured a selected subset of traits within specific functional categories. An understanding of what biological functions these regions are associated with can therefore provide hypotheses about targets of divergent natural selection for future studies. In our study, the biological processes implicated in genomic regions included developmental processes, such as development of nervous and sensory systems. Notably, visual response to abiotic stimulus and perception of sound were enriched biological functions that were associated with the community structure of invertebrates. This is particularly interesting for Mývatn stickleback due to potential for influencing ability to respond to prey stimulus, and because the sensory drive hypothesis posits that an organisms' communication system should be especially sensitive to ecological variation (Endler, 1992; Endler et al., 1993). Although adaptation can occur across exceptionally short time frames (Kingman et al., 2022), selection often acts multifariously and environments are rarely stable temporally, resulting in fluctuating selection pressures (Bell, 2010). Mývatn is a highly dynamic ecosystem, where multiple dimensions of its ecology and, importantly, stickleback population size fluctuate substantially through time (Phillips et al., 2023; Örnólfsdóttir & Einarsson, 2004; Ives et al., 2008). Our environmental, phenotypic and genomic data were collected at a single point in time, and may therefore not accurately reflect the selective environment experienced by Mývatn stickleback within or across generations. Whilst using ecological data collected at the same (single) time point as genomic and phenotypic data is common practice in studies that investigate selection in wild populations (Bolnick & Ballare, 2020; Magalhaes et al., 2021), tracking genephenotypeenvironment associations through time would allow inferences on how patterns of phenotypic and genomic variation withstand, or respond to, temporal fluctuations in selective pressures. In context of Mývatn stickleback, this would also facilitate the exploration of the spatiotemporal balance between adaptive divergence and gene flow. 5 | CONCLUSIONS Our results provide evidence for substantial genetic variation underlying functionally relevant traits, and suggest adaptive divergence in face of gene flow in a large, panmictic population inhabiting an environmentally heterogeneous lake. In particular, genephenotypeenvironment association analyses allowed us to identify genomic signatures of selection by testing which phenotypic traits and genomic variants are associated with putative selective agents. Whilst we found evidence for genomephenotypeenvironment correlations for spine length, we also found evidence for phenotypic divergence in body size (total length) and trophic morphology (gut length, gill raker length and gill raker number) without apparent genomic divergence – despite substantial additive genetic variation in these traits. The lack of genomic trait divergence across environments could reflect a combination of phenotypic plasticity and/or habitat choice (Crispo, 2008; Edelaar et al., 2017; Garant et al., 2007; Westneat et al., 2019), both of which can constrain or accelerate adaptive divergence (Crispo, 2008; Levis & Pfennig, 2020; Schoener, 1974; Wund, 2012). Whilst drawing definitive conclusions about selective processes causing genomephenotypeenvironment correlations (as opposed to drift or plasticity) would require data on phenotypefitness associations, our study sets the stage for a holistic understanding of patterns of divergence and the maintenance of genomic and phenomic variation in spatiotemporally varying wild populations. AUTHOR CONTRIBUTIONS Kasha Strickland conceptualized, designed and led the research, collected and analysed data, and led the writing of the manuscript; Katja Räsänen contributed to conceptual development; Bjarni Kristofer Kristjánsson contributed to generation of phenotype data; Joseph S. Phillips contributed to designing phenotypic analyses; Arni Einarsson collected field data; Ragna G. Snorradóttir generated the phenotype data; Mireia Bartrons contributed ecological data; Zophonías Oddur Jónsson contributed to conceptualisation and all molecular work. All authors contributed to editing the manuscript and approved the final version. ACKNO WLE DGE MENTS We thank all field volunteers who have helped to collect the data presented in this manuscript over the years. We especially thank Antoine Millet for sampling and processing stickleback, and collecting temperature data. These surveys were conducted under the auspices of the Mývatn Research Station, which has government approval for collecting fish specimens from the lake. This study was supported by the Icelandic Research Fund, grant of excellence ( 1 9 5 5 7 1 - 0 5 2 ) . CONFLICT OF INTEREST Authors have no conflicts of interest to declare. DATA AVAILABILITY AND BENEFITSHARING STATEMENT All raw reads generated and analysed in this study have been archived at the European Nucleotide Archive (ENA), accession number PRJEB58765. Phenotypic and ecological data needed to replicate analyses presented in this study have been deposited in the Dryad repository at https://doi.org/10.5061/dryad.mkkwh 7147 (Strickland et al., 2023). ORCID Kasha Strickland https://orcid.org/0000-0002-2490-0607 Zophonías Oddur Jónsson https://orcid. org/0000-0001-5798-9647 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 1723 STRICKLAND et al. REFERENCES Ålund, M., Harper, B., Kjærnested, S., Ohl, J. E., Phillips, J. G., Sattler, J., Thompson, J., Varg, J. E., Wargenau, S., Boughman, J. W., & Keagy, J. (2022). Sensory environment affects Icelandic threespine stickleback's antipredator escape behaviour. Proceedings of the Royal Society B: Biological Sciences, 289, 20220044. Ólafsdóttir, G. Á., Snorrasson, S. S., & Ritchie, M. G. (2007). Postglacial intralacustrine divergence of Icelandic threespine stickleback morphs in three neovolcanic lakes. Journal of Evolutionary Biology, 20, 1870– 1881. Ólafsson, J. (1979). Physical characteristics of Lake Mývatn and river Laxá. Oikos, 32, 38– 66. Örnólfsdóttir, E., & Einarsson, Á. (2004). Spatial and temporal variation of benthic Cladocera (Crustacea) studied with activity traps in Lake Myvatn, Iceland. Aquatic Ecology, 38, 239– 257. Albert, A. Y. K., Sawaya, S., Vines, T. H., Knecht, A. K., Miller, C. T., Summers, B. R., Balabhadra, S., Kingsley, D. M., & Schluter, D. (2008). The genetics of adaptive shape shift in stickleback: pleiotropy and effect size. Evolution: International Journal of Organic Evolution, 62(1), 76– 85. Alexa, A., & Rahnenfuhrer, J. (2020). TopGo: Enrichment analysis for gene ontology. R package version 2.42. 0. 2020. Archambeault, S. L., Bärtschi, L. R., Merminod, A. D., & Peichel, C. L. (2020). Adaptation via pleiotropy and linkage: Association mapping reveals a complex genetic architecture within the stickleback EDA locus. Evolution Letters, 4, 282– 301. Artemov, A. V., Mugue, N. S., Rastorguev, S. M., Zhenilo, S., Mazur, A. M., Tsygankova, S. V., Boulygina, E. S., Kaplun, D., Nedoluzhko, A. V., Medvedeva, Y. A., & Prokhortchouk, E. B. (2017). Genomewide DNA methylation profiling reveals epigenetic adaptation of stickleback to marine and freshwater conditions. Molecular Biology and Evolution, 34(9), 2203– 2213. Bal, T. M. P., LlanosGarrido, A., Chaturvedi, A., Verdonck, I., Hellemans, B., & Raeymaekers, J. A. M. (2021). Adaptive divergence under gene flow along an environmental gradient in two coexisting stickleback species. Genes, 12, 435. Bartrons, M., Einarsson, Á., Nobre, R. L. G., Herren, C. M., Webert, K. C., Brucet, S., Ólafsdóttir, S. R., & Ives, A. R. (2015). Spatial patterns reveal strong abiotic and biotic drivers of zooplankton community composition in Lake Mývatn, Iceland. Ecosphere, 6, art105. Bell, G. (2010). Fluctuating selection: The perpetual renewal of adaptation in variable environments. Philosophical Transactions of the Royal Society B: Biological Sciences, 365, 87– 97. Bell, M. A. (1982). Differentiation of adjacent stream populations of threespine sticklebacks. Evolution, 36, 189– 199. Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B (Methodological), 57, 289– 300. Bolnick, D. I. (2004). Can intraspecific competition drive disruptive selection? An experimental test in natural populations of sticklebacks. Evolution, 58, 608– 618. Bolnick, D. I., & Ballare, K. M. (2020). Resource diversity promotes amongindividual diet variation, but not genomic diversity, in lake stickleback. Ecology Letters, 23, 495– 505. Brodie, A., Azaria, J. R., & Ofran, Y. (2016). How far from the SNP may the causative genes be? Nucleic Acids Research, 44, 6046– 6054. CampbellStaton, S. C., Arnold, B. J., Gonçalves, D., Granli, P., Poole, J., Long, R. A., & Pringle, R. M. (2021). Ivory poaching and the rapid evolution of tusklessness in African elephants. Science, 1979(374), 483– 487. Chenoweth, S. F., Rundle, H. D., & Blows, M. W. (2010). The contribution of selection and genetic constraints to phenotypic divergence. The American Naturalist, 175, 186– 196. Cingolani, P., Platts, A., Wang, L. L., Coon, M., Nguyen, T., Wang, L., Land, S. J., Lu, X., & Ruden, D. M. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso2; i s o - 3 . Fly, 6, 80– 92. Coop, G., Witonsky, D., di Rienzo, A., & Pritchard, J. K. (2010). Using environmental correlations to identify loci underlying local adaptation. Genetics, 185, 1411– 1423. Crispo, E. (2008). Modifying effects of phenotypic plasticity on interactions among natural selection, adaptation and gene flow. Journal of Evolutionary Biology, 21, 1460– 1469. de Villemereuil, P., Frichot, É., Bazin, É., François, O., & Gaggiotti, O. E. (2014). Genome scan methods against more complex models: When and how much should we trust them? Molecular Ecology, 23, 2006– 2019. Do, C., Waples, R. S., Peel, D., Macbeth, G. M., Tillett, B. J., & Ovenden, J. R. (2014). NeEstimator v2: reimplementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Molecular Ecology Resources, 14, 209– 214. Durinck, S., Spellman, P. T., Birney, E., & Huber, W. (2009). Mapping identifiers for the integration of genomic datasets with the R/ Bioconductor package biomaRt. Nature Protocols, 4, 1184– 1191. Eckert, A. J., Bower, A. D., GonzalezMartinez, S. C., Wegrzyn, J. L., Coop, G., & Neale, D. B. (2010). Back to nature: Ecological genomics of loblolly pine (Pinus taeda, Pinaceae). Molecular Ecology, 19, 3789– 3805. Edelaar, P., Jovani, R., & GomezMestre, I. (2017). Should i change or should i go? Phenotypic plasticity and matching habitat choice in the adaptation to environmental heterogeneity. The American Naturalist, 190, 506– 520. Einarsson, Á., Stefánsdóttir, G., Jóhannesson, H., Ólafsson, J. S., Már Gíslason, G., Wakana, I., Gudbergsson, G., & Gardarsson, A. (2004). The ecology of Lake Myvatn and the river Laxá: Variation in space and time. Aquatic Ecology, 38, 317– 348. Endler, J. A. (1992). Signals, signal conditions, and the direction of evolution. The American Naturalist, 139, S125– S153. Endler, J. A. (1995). Multipletrait coevolution and environmental gradients in guppies. Trends in Ecology & Evolution, 10, 22– 29. Endler, J. A., Butlin, R. K., Guilford, T., & Krebs, J. R. (1993). Some general comments on the evolution and design of animal communication systems. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, 340, 215– 225. Feder, J. L., Egan, S. P., & Nosil, P. (2012). The genomics of speciationw i t h - g e n e - f l o w . Trends in Genetics, 28, 342– 350. Franke, F., Armitage, S. A. O., Kutzer, M. A. M., Kurtz, J., & Scharsack, J. P. (2017). Environmental temperature variation influences fitness tradeoffs and tolerance in a fishtapeworm association. Parasites & Vectors, 10, 252. Frichot, E., & François, O. (2015). LEA: An R package for landscape and ecological association studies. Methods in Ecology and Evolution, 6, 925– 929. Frichot, E., Mathieu, F., Trouillon, T., Bouchard, G., & François, O. (2014). Fast and efficient estimation of individual ancestry coefficients. Genetics, 196, 973– 983. Garant, D., Forde, S. E., & Hendry, A. P. (2007). The multifarious effects of dispersal and gene flow on contemporary adaptation. Functional Ecology, 21, 434– 443. Gardarsson, A. (1979). Waterfowl populations of Lake Mývatn and recent changes in numbers and food habits. Oikos, 32, 250– 270. Härer, A., Bolnick, D. I., & Rennison, D. J. (2021). The genomic signature of ecological divergence along the benthiclimnetic axis in allopatric and sympatric threespine stickleback. Molecular Ecology, 30, 451– 463. Hadfield, J. D. (2010). MCMC methods for multiresponse generalized linear mixed models: The MCMCglmm R package. Journal of Statistical Software, 33, 1– 22. Hendry, A. P., Bolnick, D. I., Berner, D., & Peichel, C. L. (2009). Along the speciation continuum in sticklebacks. Journal of Fish Biology, 75, 2000– 2036. 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 1724 | STRICKLAND et al. Hendry, A. P., Peichel, C. L., Matthews, B., Boughman, J. W., & Nosil, P. (2013). Stickleback research: The now and the next. Evolutionary Ecology Research, 15, 111– 141. Hoban, S., Kelley, J. L., Lotterhos, K. E., Antolin, M. F., Bradburd, G., Lowry, D. B., Poss, M. L., Reed, L. K., Storfer, A., & Whitlock, M. C. (2016). Finding the genomic basis of local adaptation: Pitfalls, practical solutions, and future directions. The American Naturalist, 188(4), 379– 397. Hudson, C. M., Lucek, K., Marques, D. A., Alexander, T. J., Moosmann, M., Spaak, P., Seehausen, O., & Matthews, B. (2021). Threespine stickleback in Lake Constance: The ecology and genomic substrate of a recent invasion. Frontiers in Ecology and Evolution, 8, 529. Ives, A. R., Einarsson, Á., Jansen, V. A. A., & Gardarsson, A. (2008). Highamplitude fluctuations and alternative dynamical states of midges in Lake Myvatn. Nature, 452, 84– 87. Jombart, T., & Ahmed, I. (2011). Adegenet 1.31: New tools for the analysis of genomewide SNP data. Bioinformatics, 27, 3070– 3071. Jombart, T., Devillard, S., & Balloux, F. (2010). Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genetics, 11, 94. Jones, M. R., Forester, B. R., Teufel, A. I., Adams, R. V., Anstett, D. N., Goodrich, B. A., Landguth, E. L., Joost, S., & Manel, S. (2013). Integrating landscape genomics and spatially explicit approaches to detect loci under selection in clinal populations. Evolution, 67, 3455– 3468. Karvonen, A., Kristjánsson, B. K., Skúlason, S., Lanki, M., Rellstab, C., & Jokela, J. (2013). Water temperature, not fish morph, determines parasite infections of sympatric Icelandic threespine sticklebacks (Gasterosteus aculeatus). Ecology and Evolution, 3, 1507– 1517. Kim, S.- Y., Costa, M. M., EsteveCodina, A., & Velando, A. (2017). Transcriptional mechanisms underlying lifehistory responses to climate change in the threespined stickleback. Evolutionary Applications, 10, 718– 730. Kingman, G. A. R., Lee, D., Jones, F. C., Desmet, D., Bell, M. A., & Kingsley, D. M. (2021). Longer or shorter spines: Reciprocal trait evolution in stickleback via triallelic regulatory changes in Stanniocalcin2a. Proceedings of the National Academy of Sciences of the United States of America, 118(31), e2100694118. Kingman, G. A. R., Vyas, D. N., Jones, F. C., Brady, S. D., Chen, H. I., Kerry, R., Milhaven, M., Bertino, T. S., Aguirre, W. E., Heins, D. C., von Hippel, F. A., Park, P. J., Kirch, M., Absher, D. M., Myers, R. M., Di Palma, F., Bell, M. A., Kingsley, D. M., & Veeramah, K. R. (2022). Predicting future from past: The genomic basis of recurrent and rapid stickleback evolution. Science Advances, 7, eabg5285. Kotrschal, A., Räsänen, K., Kristjánsson, B. K., Senn, M., & Kolm, N. (2012). Extreme sexual brain size dimorphism in sticklebacks: A consequence of the cognitive challenges of sex and parenting? PLoS One, 7, e30055. Kristjánsson, B. K., Skúlason, S., & Noakes, D. L. G. (2002). Rapid divergence in a recently isolated population of threespine stickleback (Gasterosteus aculeatus). Evolutionary Ecology Research, 4, 659– 672. Langmead, B., & Salzberg, S. L. (2012). Fast gappedread alignment with bowtie 2. Nature Methods, 9, 357– 359. Levis, N. A., & Pfennig, D. W. (2020). Plasticityled evolution: A survey of developmental mechanisms and empirical tests. Evolution & Development, 22, 71– 87. Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G., Durbin, R., & 1000 Genome Project Data Processing Subgroup. (2009). The sequence alignment/map format and SAMtools. Bioinformatics, 25, 2078– 2079. Li, Y. F., Costello, J. C., Holloway, A. K., & Hahn, M. W. (2008). “Reverse ecology” and the power of population genomics. Evolution, 62, 2984– 2994. MacColl, A. D. C. (2011). The ecological causes of evolution. Trends in Ecology & Evolution, 26, 514– 522. Magalhaes, I. S., Whiting, J. R., D'Agostino, D., Hohenlohe, P. A., Mahmud, M., Bell, M. A., Skúlason, S., & MacColl, A. (2021). Intercontinental genomic parallelism in multiple threespined stickleback adaptive radiations. Nature Ecology & Evolution, 5, 251– 261. Marrano, A., Micheletti, D., Lorenzi, S., Neale, D., & Grando, M. S. (2018). Genomic signatures of different adaptations to environmental stimuli between wild and cultivated Vitis vinifera L. Horticulture Research, 5, 34. McGee, M. D., Schluter, D., & Wainwright, P. C. (2013). Functional basis of ecological divergence in sympatric stickleback. BMC Evolutionary Biology, 13, 277. Mehlis, M., & Bakker, T. C. M. (2014). The influence of ambient water temperature on sperm performance and fertilization success in threespined sticklebacks (Gasterosteus aculeatus). Evolutionary Ecology, 28, 655– 667. Millet, A., Kristjánsson, B. K., Einarsson, Á., & Räsänen, K. (2013). Spatial phenotypic and genetic structure of threespine stickleback (Gasterosteus aculeatus) in a heterogeneous natural system, Lake Mývatn, Iceland. Ecology and Evolution, 3, 3219– 3232. Moser, G., Lee, S. H., Hayes, B. J., Goddard, M. E., Wray, N. R., & Visscher, P. M. (2015). Simultaneous discovery, estimation and prediction analysis of complex traits using a Bayesian mixture model. PLoS Genetics, 11, e1004969. Narum, S. R., & Hess, J. E. (2011). Comparison of FST outlier tests for SNP loci under selection. Molecular Ecology Resources, 11, 1 8 4 – 1 9 4 . Nath, S., Shaw, D. E., & White, M. A. (2021). Improved contiguity of the threespine stickleback genome using longread sequencing. G3 Genes|Genomes|Genetics, 11, jkab007. Nosil, P., Funk, D. J., & OrtizBarrientos, D. (2009). Divergent selection and heterogeneous genomic divergence. Molecular Ecology, 18, 375– 402. O'Brown, N. M., Summers, B. R., Jones, F. C., Brady, S. D., & Kingsley, D. M. (2015). A recurrent regulatory change underlying altered expression and Wnt response of the stickleback armor plates gene EDA. eLife, 4, e05290. Peichel, C. L., & Marques, D. A. (2017). The genetic and molecular architecture of phenotypic diversity in sticklebacks. Philosophical Transactions of the Royal Society B: Biological Sciences, 372, 20150486. Peichel, C. L., McCann, S. R., Ross, J. A., Naftaly, A. F. S., Urton, J. R., Cech, J. N., Grimwood, J., Schmutz, J., Myers, R. M., Kingsley, D. M., & White, M. A. (2020). Assembly of the threespine stickleback Y chromosome reveals convergent signatures of sex chromosome evolution. Genome Biology, 21, 177. Phillips, J., Einarsson, Á., Strickland, K., Kristjansson, B., Ives, A., & Rasanen, K. (2023). Demographic basis of spatially structured fluctuations in a threespine stickleback metapopulation. The American Naturalist, 201, 3. https://www.journals.uchicago.edu/ doi/10.1086/722741 Phillips, J. S., Guðbergsson, G., & Ives, A. R. (2022). Opposing trends in survival and recruitment slow the recovery of a historically overexploited fishery. Canadian Journal of Fisheries and Aquatic Sciences, 99(999), 1– 7. R Core Team. (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing. Räsänen, K., & Hendry, A. P. (2008). Disentangling interactions between adaptive divergence and gene flow when ecology drives diversification. Ecology Letters, 11, 624– 636. Reid, K., Bell, M. A., & Veeramah, K. R. (2021). Threespine stickleback: A model system for evolutionary genomics. Annual Review of Genomics and Human Genetics, 22, 357– 383. Reimchen, T. E. (1992). Injuries on stickleback from attacks by a toothed predator (Oncorhynchus) and implications for the evolution of lateral plates. Evolution, 46, 1224– 1230. Reimchen, T. E. (2000). Predator handling failures of lateral plate morphs in gasterosteus aculeatus: Functional implications for the ancestral plate condition. Behaviour, 137, 1081– 1096. 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 1725 STRICKLAND et al. Reimchen, T. E., & Nosil, P. (2002). Temporal variation in divergent selection on spine number in threespine stickleback. Evolution, 56, 2472– 2483. Richardson, J. L., Urban, M. C., Bolnick, D. I., & Skelly, D. K. (2014). Microgeographic adaptation and the spatial scale of evolution. Trends in Ecology & Evolution, 29, 165– 176. Rudman, S. M., Barbour, M. A., Csilléry, K., Gienapp, P., Guillaume, F., Hairston, N. G., Jr., Hendry, A. P., Lasky, J. R., Rafajlović, M., Räsänen, K., Schmidt, P. S., Seehausen, O., Therkildsen, N. O., Turcotte, M. M., & Levine, J. M. (2018). What genomic data can reveal about ecoevolutionary dynamics. Nature Ecology & Evolution, 2, 9– 15. Rypel, A. L., Layman, C. A., & Arrington, D. A. (2007). Water depth modifies relative predation risk for a motile fish taxon in Bahamian tidal creeks. Estuaries and Coasts, 30, 518– 525. Santure, A. W., & Garant, D. (2018). Wild GWAS— Association mapping in natural populations. Molecular Ecology Resources, 18, 729– 738. Schluter, D., & McPhail, J. D. (1992). Ecological character displacement and speciation in sticklebacks. The American Naturalist, 140, 85– 108. Schneider, C. A., Rasband, W. S., & Eliceiri, K. W. (2012). NIH image to ImageJ: 25 years of image analysis. Nature Methods, 9, 671– 675. Schoener, T. W. (1974). Resource partitioning in ecological communities. Science, 185, 185, 27– 139. Slate, J., Visscher, P. M., MacGregor, S., Stevens, D., Tate, M. L., & Pemberton, J. M. (2002). A genome scan for quantitative trait loci in a wild population of Red Deer (Cervus elaphus). Genetics, 162, 1863– 1873. Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, A. (2014). The deviance information criterion: 12 years on. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76, 485– 493. Stinchcombe, J. R., & Hoekstra, H. E. (2008). Combining population genomics and quantitative genetics: Finding the genes underlying ecologically important traits. Heredity, 100, 158– 170. Strickland, K., Räsänen, K., Kristjánsson, B. K., Phillips, J. S., Einarsson, A., Snorradóttir, R. G., Bartrons, M., & Jónsson, Z. O. (2023). Genomephenotypeenvironment associations identify signatures of selection in a panmictic population of threespine stickleback. Dryad, Dataset. https://doi.org/10.5061/dryad.mkkwh 7147 Vamosi, S. M., & Schluter, D. (2004). Character shifts in the defensive armor of sympatric sticklebacks. Evolution, 58, 376– 385. van der Auwera, G. A., Carneiro, M. O., Hartl, C., Poplin, R., del Angel, G., LevyMoonshine, A., Jordan, T., Shakir, K., Roazen, D., Thibault, J., Banks, E., Garimella, K. V., Altshuler, D., Gabriel, S., & DePristo, M. A. (2013). From FastQ data to highconfidence variant calls: The genome analysis toolkit best practices pipeline. Current Protocols in Bioinformatics, 43, 10– 11. Veen, T., Brock, C., Rennison, D., & Bolnick, D. (2017). Plasticity contributes to a finescale depth gradient in sticklebacks' visual system. Molecular Ecology, 26, 4339– 4350. Westneat, D. F., Potts, L. J., Sasser, K. L., & Shaffer, J. D. (2019). Causes and consequences of phenotypic plasticity in complex environments. Trends in Ecology & Evolution, 34, 555– 568. Wolf, J. B. W., & Ellegren, H. (2017). Making sense of genomic islands of differentiation in light of speciation. Nature Reviews Genetics, 18, 87– 10 0. Wund, M. A. (2012). Assessing the impacts of phenotypic plasticity on evolution. Integrative and Comparative Biology, 52, 5– 15. Yeaman, S., & Whitlock, M. C. (2011). The genetic architecture of adaptation under migration– selection balance. Evolution, 65, 1897– 1911. Zhou, X., & Stephens, M. (2012). Genomewide efficient mixedmodel analysis for association studies. Nature Genetics, 44, 821– 824. SUPPORTING INFORMATION Additional supporting information can be found online in the Supporting Information section at the end of this article. How to cite this article: Strickland, K., Räsänen, K., Kristjánsson, B. K., Phillips, J. S., Einarsson, A., Snorradóttir, R. G., Bartrons, M., & Jónsson, Z. O. (2023). Genomephenotypeenvironment associations identify signatures of selection in a panmictic population of threespine stickleback. Molecular Ecology, 32, 1708–1725. https://doi.org/10.1111/mec.16845 1365294x, 2023, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.16845 by University Of Jyväskylä Library, Wiley Online Library on [29/03/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License