Drivers of neutral and adaptive differentiation in pike (Esox lucius) populations from contrasting environments

Understanding how eco‐evolutionary processes and environmental factors drive population differentiation and adaptation are key challenges in evolutionary biology of relevance for biodiversity protection. Differentiation requires at least partial reproductive separation, which may result from different modes of isolation such as geographic isolation (allopatry) or isolation by distance (IBD), resistance (IBR), and environment (IBE). Despite that multiple modes might jointly influence differentiation, studies that compare the relative contributions are scarce. Using RADseq, we analyse neutral and adaptive genetic diversity and structure in 11 pike (Esox lucius) populations from contrasting environments along a latitudinal gradient (54.9–63.6°N), to investigate the relative effects of IBD, IBE and IBR, and to assess whether the effects differ between neutral and adaptive variation, or across structural levels. Patterns of neutral and adaptive variation differed, probably reflecting that they have been differently affected by stochastic and deterministic processes. The importance of the different modes of isolation differed between neutral and adaptive diversity, yet were consistent across structural levels. Neutral variation was influenced by interactions among all three modes of isolation, with IBR (seascape features) playing a central role, wheares adaptive variation was mainly influenced by IBE (environmental conditions). Taken together, this and previous studies suggest that it is common that multiple modes of isolation interactively shape patterns of genetic variation, and that their relative contributions differ among systems. To enable identification of general patterns and understand how various factors influence the relative contributions, it is important that several modes are simultaneously investigated in additional populations, species and environmental settings.


| INTRODUC TI ON
Understanding how eco-evolutionary processes and environmental factors drive population differentiation and adaptation remain key challenges in evolutionary biology. In addition to informing about the origin and dynamics of biodiversity, knowledge about drivers can promote successful management and protection of populations and species. For differentiation to occur populations must be partially reproductively separated. That geographic isolation imposes reproductive separation is self-evident. However, reproductive isolation may also arise through other mechanisms. For example, the isolation by distance model (IBD; Wright, 1943) proposes that there is a negative association between geographic distance and gene flow, which should translate into clinal population structures with higher degrees of differentiation on longer distances. Landscape or seascape characteristics (e.g., roads, mountains, open waters, and currents) can further influence the feasibility of dispersal (resistance), and the isolation by resistance model (IBR;McRae, 2006) proposes that in heterogeneous environments there should be a negative relationship between resistance and gene flow. In addition to geography and landscape/seascape properties, ecological factors can influence the degree of gene flow (Rundle & Nosil, 2005). Differences in ecology and environmental preferences (e.g., different ecotypes) can instigate population differentiation via isolation by environment (IBE; Wang & Bradburd, 2014). There is also potential for genetic structuring to evolve as a result of isolation by time (IBT), whereby reproductive separation occurs between groups of individuals due to differences in reproductive timing (Hendry & Day, 2005).
Differences in ecology and environmental preferences are generally accompanied by divergent selection and the evolution of different local adaptations, which can cause a feed-back loop that reinforces the reproductive separation and speeds up the differentiation process (a case of IBE referred to as isolation by adaptation [IBA]; Nosil et al., 2008). The reinforcing effect is probably of particular importance when divergent selection acts upon traits directly associated with reproductive isolation and reproductive success, for example, habitat preferences (Hendry et al., 2007), spawning segregation (Nosil, 2012), or early life history traits (Momigliano et al., 2017), and might lead to "ecological speciation" (Schluter & Rambaut, 1996).
There is a growing awareness that multiple modes of isolation can interactively influence differentiation. In a review of 70 studies, Sexton et al. (2014) found evidence of IBD in 20.0%, IBE in 37.1%, and both patterns in 37.1%. Despite that more than a third of the examples found evidence of both IBD and IBE, attempts to investigate the relative importance of different mechanisms are scarce.
Northern temperate freshwater fish species present good opportunities to study differentiation and eco-evolutionary dynamics (Hume et al., 2018), and to compare the relative effects of different drivers within study systems. It was not until after the retreat of the glaciers (15,000-10,000 years ago) that many areas were colonized via range expansions (Petit et al., 2003;Rowe et al., 2004). The geographical isolation in glacial refugia and subsequent range expansions have allowed for founder events, bottlenecks, divergent selection, and adaptions to local conditions to influence genetic structure and diversity. In aquatic systems, gene flow is generally expected to be high due to the lack of apparent dispersal boundaries (Gagnaire et al., 2015;Puebla et al., 2012). Despite this, genetic structuring has been reported within open and connected waterbodies (Bergek & Björklund, 2007;Momigliano et al., 2017;Nordahl et al., 2019;Willing et al., 2010;Yıldırım, Anderson, et al., 2018), even on small spatial scales. One example is the Baltic Sea, one of the largest brackish-water ecosystems in the world, that exhibits north-south gradients in salinity and temperature (Bendtsen et al., 2007). This makes the Baltic Sea an excellent system for studying differentiation and local adaptation (e.g., in Guo et al., 2015Guo et al., , 2016. One of the most common, coastal, predatory fish species in the Baltic Sea is pike (Esox lucius). It is of freshwater origin but has colonised brackish waters with salinities up to approximately 15 ppt (Craig, 2008). This range expansion has been accompanied by the evolution of three different ecotypes that differ with regard to migratory behaviour. The original (freshwater) ecotype spends the entire life in freshwater, whilst the other ecotypes (anadromous and resident) spend part or their entire life in saline (brackish) waters. The two latter ecotypes coexist in the coastal waters of the Baltic Sea during most of the year (Westin & Limburg, 2002), and separate only for a short period during spawning when anadromous individuals migrate to freshwater localities whilst the resident ecotype stays in the brackish coastal waters (Engstedt et al., 2010;Larsson et al., 2015;Muller, 1986). An important feature of Baltic Sea pike, which affects the population structure of the species, is their homing behaviour (Engstedt et al., 2014;Larsson et al., 2015;, with adults returning to their natal spawning ground for reproduction (Engstedt et al., 2010(Engstedt et al., , 2014Jacobsen et al., 2017;Muller, 1986). Another important aspect that has been suggested to influence population structuring in pike is that open pelagic environments act as dispersal barriers (Nordahl et al., 2019), which probably reflect that pike are associated with macrophytes and structured habitats, and negatively select open pelagic habitats.
While a previous study of Baltic Sea pike has indicated genetic differentiation between the resident and anadromous ecotypes (Nordahl et al., 2019), studies of genetic structuring within the ecotypes show conflicting results. Some studies suggest weak structuring (Laikre et al., 2005;Wennerström et al., 2016), whilst others report fine-scaled genetic structuring among anadromous populations, and point to a role of isolation by distance (Bekkevold et al., 2015;Möller et al., 2020;Nordahl et al., 2019;Sunde et al., 2020a).
There is also evidence to suggest that environmental differences among spawning locations have resulted in adaptive phenotypic differentiation in salinity tolerance (Sunde et al., 2018), temperature tolerance , vertebral number , body size and growth rate , early life history traits, and reproductive investment . Pike in the Baltic Sea therefore offers possibilities to study differentiation at different structural levels (between allopatric and sympatric populations, and within and among ecotypes), to evaluate the relative contributions to genetic structure of different mode of isolation (such as IBD, IBR, IBE, IBA, and IBT), and investigate the genetic underpinnings of local adaptations.
Previous population assignments of Baltic Sea pike have mainly been based on neutral microsatellite markers (Bekkevold et al., 2015;Eschbach et al., 2021;Möller et al., 2020;Nordahl et al., 2019;Wennerström et al., 2016). However, knowledge about how different evolutionary processes shape genetic structure and diversity requires that population genetic studies use molecular markers that capture variation also in coding regions (Andrews et al., 2016;Candy et al., 2015;Ford et al., 2015;Holderegger et al., 2006;Sunde et al., 2020a). Restriction site associated DNA sequencing (RADseq) yields thousands of single nucleotide polymorphisms (SNPs) residing in both coding (functional) and non-coding (mainly selectively neutral) genomic regions (Andrews et al., 2016). Comparisons of genetic structure and diversity obtained using separate analyses of neutral and adaptive SNPs can inform about connectivity and gene flow, provide insights about the roles of both stochastic and deterministic processes, allow for identification of candidate genes influenced by selection by environmental factors (Andrews et al., 2016;de Villemereuil & Gaggiotti, 2015;, and therefore have potential to provide a comprehensive understanding of the eco-evolutionary processes that jointly influence genetic diversity and shape genetic structure of natural populations. Here, we report on the results of a population genetic study utilizing the RADseq method to investigate genetic structure and decipher the roles of different ecological and evolutionary processes for differentiation and adaptation in pike. For this we used 11 populations representing the three ecotypes and spanning 8.7 degrees in latitude, from Denmark in the south to Umeå in northern Sweden ( Figure 1, Table 1), that experience different environmental conditions. We investigated genetic structure using both neutral and outlier SNP datasets, and evaluated the relative contributions of IBE, IBD and IBR at different structural levels (between allopatric and sympatric populations, and within and among ecotypes). We also performed outlier analyses to identify loci putatively under selection and to pinpoint specific environmental factors contributing to evolutionary divergence.

| Study species
Pike is a large, long-lived, iteroparous predatory fish that plays an important role in many aquatic ecosystems by regulating abundances of species in lower trophic levels (Donadi et al., 2017). It is economically important for recreational and commercial fishing (Lehtonen et al., 2009;Pierce et al., 1995), and an established model F I G U R E 1 Map of study area, and genetic structuring among the pike, Esox lucius, populations. The map shows the locations and spawning ecotypes of the study populations (blue, freshwater; black, anadromous; and red, resident). It was created in Adobe Photoshop CC, v. 2015.0.1, and is a combined and modified version of two base maps (one of Scandinavia, and one of Sweden) that are available from Wikimedia Commons under the nonrestrictive creative commons license. The distruct plots show the genetic structuring among the populations for the most likely number of genetic clusters (K) suggested by fastSTRUCTURE for the full (neutral) and outlier datasets for 234 individuals [Colour figure can be viewed at wileyonlinelibrary.com] species in studies of ecology and evolution .
Increased knowledge is required to aid successful management.

| Study populations and sampling procedure
A total of 234 pike from 11 populations (N = 8-27) were sampled for this study. The locations of the populations ranged from 63.6-54.9°N, and included the three spawning ecotypes (freshwater, anadromous, and resident; for details see Figure 1 and Table 1). Our initial plan was to include both anadromous and resident populations throughout the latitudinal range. However, due to the aforementioned decreases of pike in the Baltic Sea, many of the resident populations are no longer found in their former spawning locations, and we were therefore only able to include one resident population. Individuals were captured using fyke nets, rod-and-reel fishing, or electrofishing, and a nonlethal DNA sample (fin clip) was collected for each individual before they were released back into the water. The fin clips were placed in separate 1.5 ml Eppendorf tubes with 70% ethanol, which were stored in a freezer (-20°C) until the molecular work was conducted.

| Molecular workflow
DNA was extracted from fin tissue with DNeasy blood and tissue kit (Qiagen), and digested with HF EcoRI (New England Biolabs).  (Sunde et al., 2020a) and retrieved from the NCBI Sequence Read Archive (Sunde et al., 2020b). In total this yielded a dataset of ~3200 M raw reads (mean 13.6 M per sample).
Quality filtering of the raw reads was conducted using Trimmomatic (Bolger et al., 2014) and process_radtags in the Stacks pipeline (Catchen et al., 2011(Catchen et al., , 2013. The ~3000 M reads (95%) that passed the quality filtering, were further processed with the integrated approach of the Stacks pipeline (using version 2.2 of Stacks) described in (Paris et al., 2017). The integrated approach is a combination of the de novo and the ref_map (reference) pipelines in which consensus sequences obtained from the de novo runs are aligned to the reference genome to integrate the alignment information into the de novo output (Paris et al., 2017). This method was chosen as it has been shown to be more successful than the ref_map pipeline (Paris et al., 2017). To determine parameter settings, we performed an initial parameter optimization (Figures S1 and S2) before the integrated approach was carried out with the parameter settings for m, M, and n as suggested by the optimization (4 for all), and otherwise default settings. In short, in the integrated approach ustacks, cstacks, sstacks, and gstacks were used to create consensus sequences and detect SNPs. The consensus sequences were aligned to the refer-

| Analysis of neutral genetic diversity and structure
The Populations software in the Stacks pipeline was used to assess whether the populations were in Hardy-Weinberg equilibrium (HWE). Only a small proportion of the loci deviated from HWE (100-283 loci for each population, P < 0.05) and the F is distributions were unimodal and peaked at zero for all populations ( Figure   S3), suggesting that the majority of loci were not affected by null alleles and that the populations were in HWE. All loci were therefore retained for the subsequent analyses. The same software was used to obtain estimates of number of private alleles, observed het-

erozygosity (H O ), expected heterozygosity (H E ), and fixation index
F is for each population, and to assess pairwise population differentiation based on fixation index F ST by Weir and Cockerham (1984).
To test whether genetic diversity differed among populations a distance matrix of pairwise genetic similarities (proportion of alleles shared) was analysed with the PERMDISP analysis (Anderson, 2006) in the PERMANOVA+software in PRIMER (Anderson et al., 2008) to test whether the average dispersion from the centroid (estimate of within-population variation; Parker et al., 2012;Yıldırım, Tinnert, et al., 2018) differed among populations. Within-population diversity between freshwater and anadromous populations was compared using a t-test.
To assess genetic structuring and determine the most likely number of genetic clusters (K), the full dataset was analysed using fast-STRUCTURE (Raj et al., 2014). Because fine-scaled differentiation might be concealed by differentiation on higher levels (Evanno et al., 2005), we also ran the analysis with a subset of only the anadromous populations to test if further differentiation became evident. For details see "Analysis of genetic diversity and population structure" in the Supporting Information.
To assess whether the full dataset was representative of neutrally evolving diversity, the results were compared to those generated by a subset "neutral" dataset. This "neutral" dataset was created by excluding the 5% tails of the F ST distribution from the full dataset, which should exclude loci under strong selection (both divergent and balancing). Patterns of genetic structure indicated by such a neutral dataset should therefore be reflective of stochastic processes. The comparison revealed that the results obtained for the full and the "neutral" datasets were qualitatively similar (not shown), indicating that the full dataset was representative of neutral diversity. We therefore proceeded with the full dataset for the comparisons with adaptive genetic variation and structure (see below).

| Phylogenetic analysis
The evolutionary relationship among the samples was investigated with a maximum likelihood based phylogenetic inference in RAxML v8.2.10 (Stamatakis, 2014). This was done by using the full dataset following the steps described in "The Hard & Slow Way" in the RAxML manual (Stamatakis, 2016), which includes optimizing the parameter settings for number of rate categories and initial rearrangement. After the optimization, the final analysis with 200 inferences and 1000 bootstraps was run using 10 as initial rear-

| Testing for effects of different modes of isolation
To investigate the relative roles of different modes of isolation for neutral genetic structure, the full dataset was analysed using distance-based redundancy analysis (db-RDA; Legendre & Anderson, 1999) and maximum likelihood population effect (MLPE) mixed modelling (Clarke et al., 2002) in RStudio 2 (RStudio Team, 2015) with R (R Core Team, 2012). For this, we created a genetic distance matrix based on pairwise genetic similarity between individuals (proportion of alleles shared), and geographic distance, resistance, and environmental dissimilarity were estimated to test for effects of IBD, IBR and IBE, respectively.
For geographic distance, latitude was used as a proxy (   fastSTRUCTURE suggested that the most likely number of genetic clusters (K) was either 6 (for model complexity that maximizes marginal likelihood) or 7 (for model used to explain structure in the data) ( Figure 1). The results for K = 6 and K = 7 were generally similar, and only differed with regard to K = 7 assigning Snäckstavik to a population-specific cluster, whilst it was assigned to a shared cluster for K = 6. All three freshwater lake populations were mainly assigned to population-specific clusters, which

| Loci putatively under selection
All outlier analysis approaches identified loci putatively under selection. Locus-specific effects were found for 28 loci with BayeScan (q-values < 0.05), 231 loci with Fdist (P < 0.01), and 635 loci with LOSITAN (P < 0.01). Of all these loci, 17 were identified by all three software ( Figure S7). The LD window approach showed that seven of the outliers SNPs were located within previously annotated genes, and that 26 genes in total resided within 20 kb from the genomic positions of the outliers (and were thus considered as putatively under selection) (for details see Table S1).
When utilizing BayeScEnv to search for associations between allele frequencies and the environmental variables (midrange salinity and latitude), 13 loci were identified as outliers (q-value < 0.05; Figure S7). Eight of the loci were located within previously annotated genes and 23 candidate genes in total were covered by the LD decay windows (for details see Table S2). Of these 23 genes, five were also identified in the test of locus-specific effects.

| Neutral versus adaptive diversity
We found that the patterns of neutral and adaptive genetic diversity differed. More specifically, the neutral variation varied among populations and was generally higher in the freshwater populations, whereas the level of adaptive variation was similar for most of the probably reflects that pike is of freshwater origin (Craig, 2008), and that younger populations tend to have lower genetic diversity due to founder effects and genetic drift (García-Verdugo et al., 2009;Haag et al., 2005). The finding that intrapopulation adaptive diversity was not associated with neutral diversity across populations is consistent with the findings in previous studies (Leinonen et al., 2008;Reed & Frankham, 2001;Whitlock et al., 2013;Yıldırım, Tinnert, et al., 2018), and indicates that different processes have shaped the patterns of neutral and adaptive variation, as expected based on evolutionary theory. A key challenge for conservation is to design actions that maintain functional genetic and phenotypic diversity both within and among populations (Hutchinson, 2008;Stephenson, 1999;Tamario et al., 2019). The realization that estimates of neutral genetic diversity cannot reliably inform about evolvability and viability of populations is of utmost importance for successful management.
That the adaptive variation was higher than the neutral variation indicates that multiple alleles for the outlier loci (associated with candidate genes) are present in the populations. This is in contrast to what would be expected by the traditional selective sweep model, which suggests that favourable alleles rapidly should reach fixation (Smith & Haigh, 1974). A possible explanation for the higher adaptive variation is that the traits associated with the outlier loci might be polygenic, that is, controlled by multiple genes (Roff, 1997). The selection pressure on each locus will therefore not be that pronounced, and the adaptive responses associated with such traits commonly manifest as smaller allele frequency shifts at many (hundreds to thousands) loci (Pritchard & Di Rienzo, 2010;. The higher adaptive diversity could also reflect temporal environmental heterogeneity, with fluctuating selection having favored different genotypes across time and thus retained a variety of genotypes within the populations (Bell, 2010;Hedrick, 2006;Roff, 1997). This fits well with the finding that the Lervik population had higher adaptive variation than Harfjärden, as Sunde et al. (2018) showed that Lervik experience a temporally unstable salinity regime whilst Harfjärden is temporally stable. However, information about temporal heterogeneity in the other study locations is lacking, thus preventing formal evaluation of this hypothesis. The finding that the populations have a variety of alleles in the outlier loci nevertheless points to an adaptive potential which could be vital for their capability to cope with future environmental changes associated with anthropogenic impacts such as eutrophication and climate change.
Genetic and phenotypic variation is required for populations to respond to selection and adapt to changing and novel environmental conditions (Charlesworth & Charlesworth, 2017;Forsman, 2014;Roff, 1997;Wennersten & Forsman, 2012). There is also potential for the consequences of genetic variation to go beyond the level of the species, as it can influence community structure and ecosystem functioning (Des Roches et al., 2018;Hughes et al., 2008). Being important predators, competitors, and prey to other species, there are many ways by which pike and other species of fish can affect the functioning of lakes, rivers, coastal ecosystems, and open oceans (Donadi et al., 2017;Nilsson et al., 2019;Post et al., 2008;Tamario et al., 2019). To maintain ecosystem function, evolutionary potential should hence be safe-guarded within and across species.

| Genetic differentiation
All populations were significantly differentiated, and the differentiation ranged from low to high (F ST : 0.06-0.25; Figure 3), which is consistent with previous studies (Bekkevold et al., 2015;Larsson et al., 2015;Möller et al., 2020;Nordahl et al., 2019;Sunde et al., 2020a;Wennerström et al., 2016). The F ST values in the present study were high compared to those reported in other large scale studies, and this despite that Weir and Cockerham F ST tends to underestimate differentiation when based on large SNP datasets (Bhatia et al., 2013). Notably, studies on smaller spatial scales report higher levels of differentiation (

| Gene flow
Despite that all populations were significantly neutrally differentiated, the anadromous Swedish mainland populations clustered together (Figure 1). This pattern is consistent with what has been reported in a previous large-scale study of pike (Laikre et al., 2005), and a similar pattern has also been shown for pikeperch (Sander lucioperca) in the northern part of the Baltic Sea (Säisä et al., 2010).
Clustering of populations is indicative of gene flow and/or recent divergence. Given that all anadromous populations were differentiated from each other (Figure 3), that the genetic clustering was not as apparent when only the anadromous populations where analysed ( Figure S6), and that previous studies have reported low levels of gene flow among anadromous pike populations (Nordahl et al., 2019;Sunde et al., 2020a;Tibblin et al., 2015), it seems unlikely that gene flow would be sufficient to explain the observed clustering, and instead points to a role of more recent divergence.

| Post-glacial expansions
Previous studies based on mitochondrial DNA from pike across

| Interactive effects of isolation by distance, environment and resistance
We found that all three modes of isolation (IBD, IBE and IBR) affected neutral genetic distance. The significant interaction effects for both IBR x IBD and IBR x IBE ( Table 2) further indicated that a complex interplay between distance, environment and seascape have shaped the neutral genetic structuring in our study system. However, the relative importance of the modes differed, and resistance (seascape features) appears to be the main driver of the neutral genetic structuring by preventing gene flow. This was evidenced by the much larger proportion of variation explained by IBR (35.2%-36.7%) as compared to both IBD (1.5%) and IBE (3.0%-3.2%) ( Table 2), and was also corroborated by that the highest levels of neutral differentiation were found in the comparisons including populations characterised by high resistance (the anadromous population Harfjärden or at least one of the three freshwater populations) (Figure 3).
When all populations were considered, the effect of IBR was mainly attributable to the geographic isolation of the freshwater populations ( Figure 5a). It has been suggested that pike experienced drastic population declines or bottlenecks following postglacial recolonization across the Northern Hemisphere (Jacobsen et al., 2004), and that the succeeding isolation of local river and lake systems might have resulted in genetic drift and differentiation among freshwater populations (Bekkevold et al., 2015). Given this, and that the three freshwater pop- Some studies suggest that resistance is the most important mechanism affecting neutral genetic variation (Camurugi et al., 2021;Oliveira et al., 2019;Xuereb et al., 2018). Other studies point to a main effect of the environment (Jiang et al., 2019;Quéméré et al., 2016;Rödin-Mörch et al., 2019), and yet others indicate that distance is the main determinant of genetic structure (Ferreira et al., 2020;Vidal et al., 2019).
The relative importance of the different modes of isolation could depend on environmental conditions, ecology of the species, landscape/ seascape properties, sampling scale and the specific set of samples (individuals, populations, and/or lineages) studied (Ferreira et al., 2020;Huang et al., 2016;Ruiz-Gonzalez et al., 2015). The connectivity among populations, environmental heterogeneity, dispersal capacity, longevity, environmental requirements, and whether generations are discrete or overlapping that influence the evolutionary dynamics of populations (Roff, 1997) might also modulate the effects of the different drivers.
Unfortunately, the few studies evaluating the contributions of different modes of isolation within study systems differ in several important aspects (e.g., study species and the inherent differences in their ecology, spatial scale, and modes studied) which complicates systematic evaluations. Firm conclusions are therefore not yet possible.

| Drivers of adaptive genetic differentiation
As for genetic diversity, we found that the structuring and main modes of isolation differed between neutral and adaptive variation. The adaptive genetic structuring was not as clear as the neutral structuring ( Figure 1). This may partly reflect the much lower number of loci (17 in the outlier dataset compared to 5993 in the full [neutral] dataset). It could also, in theory, reflect that the outliers were not truly adaptive loci, but rather represented outlier loci arisen from random processes such as genetic drift. However, that the adaptive intrapopulation variation was higher than the neutral variation indicates that the outliers did not represent loci with low variability due to drift. That two of the strongly neutrally differentiated freshwater populations (Kosta and Hamnaryd) shared the same adaptive genetic cluster (Figure 1) strengthens the inference that the outliers represent loci associated with common adaptations to similar environments.

| Adaptive isolation by environment mediated by latitude and salinity
Despite that the adaptive structuring was weaker, there were associations with latitude and midrange salinity (Figure 1). Together with the finding that the only single mode of isolation that was identified by both approaches (MLPE and db-RDA) as significantly affecting the adaptive structuring was IBE ( Even though IBE had the highest relative importance, it accounted for no more than 6% of the variation (Table 2), meaning that a large proportion of the adaptive structure is not explained by any of the environmental variables considered in this study. This might result from an incomplete representation of the population-specific selective environments. We used midrange salinity as a proxy for the environment. Previous studies have shown that salinity is an important environmental factor for pike (Jørgensen et al., 2010;Sunde et al., 2018), and our present outlier analyses identified candidate genes associated with salinity tolerance. However, selection imposed by additional environmental variables (e.g., temperature, This population experiences divergent environmental conditions during spawning, which has resulted in the evolution of local adaptations during early fry development (temperature, Sunde et al., 2019;salinity, Sunde et al., 2018). The island population generally spawns earlier than the mainland populations . It is therefore possible that in addition to IBE, IBT also contributes to the high level of adaptive differentiation of Harfjärden. To the extent that the timing of spawning migration is heritable , differences in timing of reproduction among populations may impair the success of individuals that attempt to spawn in a population different from where they were born, and thereby reduce gene flow. However, further studies would be required to thoroughly evaluate potential contributions of IBT.

| Added value of analysing both neutral and adaptive variation
Our study revealed contrasting patterns and identified different isolating mechanisms as important drivers of neutral and adaptive genetic variation. Specifically, neutral genetic clustering of the anadromous Swedish mainland populations was observed, whereas adaptive structuring showed a pattern seemingly associated with midrange salinity. Our results add to the existing body of evidence that estimates of neutral genetic diversity cannot be used as substitutes for adaptive or functional genetic diversity to infer evolutionary potential and the ability of populations and species to cope with environmental change (Leinonen et al., 2008;Reed & Frankham, 2001;Whitlock et al., 2013;Yıldırım, Tinnert, et al., 2018). This also exemplifies how the utilization of neutral and functional markers together can provide a more comprehensive understanding of the ecoevolutionary processes that jointly influence genetic diversity and structure of natural populations (Andrews et al., 2016;Holderegger et al., 2006;Rödin-Mörch et al., 2019). In particular, our outlier analyses identified some specific genes associated with adaptative structure, the functions of which also point to certain selective pressures, as discussed below.

| Candidate genes putatively under selection
Some of the identified candidate genes putatively under selection have been found to be associated with salinity tolerance. One of these genes was potassium voltage-gated channel subfamily A member 10 (KCNA10), which encodes a voltage-dependent potassium-selective channel, and which has been reported to be associated with salinity stress in blue mussels (Mytilus spp.) (Lockwood & Somero, 2011).
Another gene was vesicle-associated membrane protein 7 (VAMP7), which is crucial for calcium regulated lysosomal exocytosis, and has been found to be involved in salt-tolerance in Arabidopsis thaliana (Leshem et al., 2006). These genes therefore seem particularly interesting as candidate genes involved in adaptation to salinity. Previous studies also report on evolution associated with salinity tolerance in other fish species in the Baltic Sea, including three-spined stickleback (Gasterosteus aculeatus) (Guo et al., 2015;Hasan et al., 2017) and European flounder (Platichthys flesus) (Momigliano et al., 2017), emphasizing the general importance of salinity adaptation.
Another gene that seems particularly interesting is zinc-finger protein 436-like (ZNF436-like), identified as putatively under selection also in a previous study of three of our study populations (Sunde et al., 2020a). The exact function of ZNF436 in fish is not known, but it is a transcription factor that regulates early cardiac development in humans (Fu et al., 2018). Other transcription factors have been found to be important in responses to heat-stress in sea urchin (Strongylocentrotus intermedius) (Zhan et al., 2019), and some ZNF genes are associated with acclimation to low salinity in the euryhaline fish half-smooth tongue sole (Cynoglossus semilaevis) (Si et al., 2018) and migratory behaviour in brown trout (Salmo trutta) (Lemopoulos et al., 2018). It is also worth to note that in pike the exact function of several of the candidate genes are not known, and that the SNPs for which no candidate genes were found might be associated with genes that have not yet been identified in the pike genome. It is therefore likely that additional environmental variables not identified here also impose selection that contributes to the adaptive population genetic structure.

| CON CLUS I ON S AND FUTURE DIREC TIONS
We utilized RADseq to study neutral and adaptive genetic variation and structure in 11 populations of E. lucius pike from contrasting environments. The results show that in addition to distance, both the environment and landscape/seascape can play important roles in shaping patterns of genetic variation, and exemplify how a combination of analyses of neutral and adaptive variation can help disentangle the complex interplay of different stochastic and deterministic contributing processes. Not only were the patterns of neutral and adaptive genetic variation different, the importance of different modes of isolation also differed between the two, yet were consistent across structural levels (between allopatric and sympatric populations, and within and among ecotypes). In our pike study system, neutral genetic variation was influenced by an interplay of all three modes of isolation (IBD, IBR and IBE), with seascape features (complete and partial dispersal barriers) being the single mode of highest importance, whereas adaptive variation was mainly influenced by the environment. The gene-environment associations and outlier loci suggest that the adaptive structure was largely driven by selection targeting genes associated with salinity tolerance. However, the exact functions of several of the outlier loci are yet not known, and selection imposed by other factors (e.g., temperature) may also be important. Finally, (iii) we note that studies comparing the contributions of different modes of isolation are still scarce. As more studies become available, systematic evaluations and meta-analysis approaches may be used to assess the generatlity of findings and determine whether and how the roles of IBD, IBE, IBR and IBT are modified by species characteristics, spatial scales and environmental settings.

ACK N OWLED G EM ENTS
We would like to thank Anssi Laurila, Dennis Hasselquist, Ane Timenes

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

E TH I C A L A PPROVA L
The study was carried out in accordance with all relevant applicable national guidelines for the care and use of animals. The study was

B EN EFIT-S H A R I N G S TATEM ENT
Benefits from this research accrue from the sharing of our data and results on public databases as described above.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at the Dryad digital repository: https:// doi.org/10.5061/dryad.sf7m0 cg79 (Sunde et al., 2021a), and the NCBI Sequence Read Archive (SRA): BioProjects PRJNA586770 (submission ID: SUB6466329; Sunde et al., 2020b) and PRJNA639676 (BioID: SAMN 15246073 -15246242; Sunde et al., 2021b).

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw sequence reads and metadata have been made available in the Sequence Read Archive (SRA) in NCBI (BioProject: PRJNA639676 and BioID: SAMN 15246073 -15246242; Sunde et al., 2021b), and the vcf genotype file is available in the Dryad Digital Repository (Sunde et al., 2021a).