Parallel evolution and adaptation to environmental factors in a marine flatfish: Implications for fisheries and aquaculture management of the turbot (Scophthalmus maximus)

Abstract Unraveling adaptive genetic variation represents, in addition to the estimate of population demographic parameters, a cornerstone for the management of aquatic natural living resources, which, in turn, represent the raw material for breeding programs. The turbot (Scophthalmus maximus) is a marine flatfish of high commercial value living on the European continental shelf. While wild populations are declining, aquaculture is flourishing in southern Europe. We evaluated the genetic structure of turbot throughout its natural distribution range (672 individuals; 20 populations) by analyzing allele frequency data from 755 single nucleotide polymorphism discovered and genotyped by double‐digest RAD sequencing. The species was structured into four main regions: Baltic Sea, Atlantic Ocean, Adriatic Sea, and Black Sea, with subtle differentiation apparent at the distribution margins of the Atlantic region. Genetic diversity and effective population size estimates were highest in the Atlantic populations, the area of greatest occurrence, while turbot from other regions showed lower levels, reflecting geographical isolation and reduced abundance. Divergent selection was detected within and between the Atlantic Ocean and Baltic Sea regions, and also when comparing these two regions with the Black Sea. Evidence of parallel evolution was detected between the two low salinity regions, the Baltic and Black seas. Correlation between genetic and environmental variation indicated that temperature and salinity were probably the main environmental drivers of selection. Mining around the four genomic regions consistently inferred to be under selection identified candidate genes related to osmoregulation, growth, and resistance to diseases. The new insights are useful for the management of turbot fisheries and aquaculture by providing the baseline for evaluating the consequences of turbot releases from restocking and farming.


| INTRODUC TI ON
The detection of genetic structure in marine species represents a challenge due to generally high effective population sizes and high gene flow facilitated by the absence of physical barriers, which lead to genomic homogenization across populations (Danancher & Garcia-Vazquez, 2011;Vandamme et al., 2014;Vilas et al., 2015).
However, various factors can bring about genetic differentiation, such as habitat shifts (ecotones) and oceanic currents (Blanco-Gonzalez, Knutsen, & Jorde, 2016;Galarza et al., 2009;Nielsen, Nielsen, Meldrup, & Hansen, 2004;Vera et al., 2016a), and natural selection in response to environmental variation (Milano et al., 2014;Vandamme et al., 2014;Vilas, Bouza, Vera, Millán, & Martínez, 2010;Vilas et al., 2015). Distinguishing between neutral and adaptive genetic variation has become a central issue in evolutionary biology, allowing for understanding of population structure in both historical/demographic and adaptive terms (Bernatchez, 2016;Nielsen, Hemmer-Hansen, Larsen, & Bekkevold, 2009), thereby providing essential information for the conservation and management of wild populations. Genetic diversity in the wild represents, in turn, the raw material for the foundation of aquaculture broodstock and consequently, a reference to identify selection signatures for targeted traits in farmed populations through genome scanning (Liu et al., 2017).
The turbot, Scophthalmus maximus L., is a marine flatfish of the family Scophthalmidae, Order Pleuronectiformes, which lives in the Northeast Atlantic Ocean (from Morocco to the Arctic Circle) and in the Mediterranean Sea as well as in the Black Sea (Froese & Pauly, 2016). It has a demersal lifestyle and inhabits sandy coastal habitats (Bouza et al., 2014). Postlarval stages are relatively sedentary, with generally short migration distances (<10 km) being reported for both juvenile and adults of the species (Florin & Höglund, 2007;Nielsen et al., 2004). In contrast, the high dispersal potential of pelagic larvae (until ~30 days posthatching) mediated by oceanic currents coupled with the high fecundity provides potential for gene flow on larger spatial scales (Johannesson & André, 2006). Turbot is currently classified as a vulnerable species according to the IUCN European Red List assessment (Nieto et al., 2015). In some Atlantic regions, turbot fisheries are close to depletion and its main fisheries are located in the North Sea and in the Baltic Sea (ICES 2017a(ICES , 2017b, where turbot is exploited as a by-catch species. An analysis of historical survey data in the North Sea suggests that turbot biomass has importantly declined, which might be associated with important biological changes in growth rate and reproduction (Bouza et al., 2014). In the Black Sea, the turbot is one of the most valuable commercial species and it is subjected to intensive fishing which has led to be characterized as exploited unsustainably and at risk of collapse (Nikolov et al., 2015). This scenario has led to restocking depleted areas with hatchery turbot with unknown outcomes, as its monitoring has not been accomplished (Bouza et al., 2014). On the other hand, as turbot breeding programs are at its beginning , it is still feasible to introduce genetic variation from natural resources, especially in new geographical areas with particular environmental conditions. Although most turbot farms are located in the Atlantic area of Spain, France, and Portugal, its high commercial value is promoting new facilities in the Adriatic Sea (Croatia) and in the Black Sea (Turkey; FEAP, 2013).
Turbot experience a diverse physical and biological environment across its range. The Atlantic Ocean has a subtle salinity gradient running roughly from north to south, while sharp differences are found between the Northern Atlantic Ocean (≈35 PSU-practical salinity units) and the Baltic Sea (up to ≈2 PSU in the northern area; Environmental Marine Information System (EMIS) database; http:// mcc.jrc.ec.europa.eu/emis/). Within the Baltic Sea, excluding the transition area with the North Sea, salinity can also reach higher values in the south (≈13 PSU). In the Mediterranean Sea, the salinity is even higher than in the Atlantic Ocean (≈38 PSU), but drops abruptly in the transition to the Black Sea, whose salinity levels resemble the Baltic Sea (≈11 PSU). Contrasting patterns of surface temperature also occur across latitude and between seasons. A north-south cline exists in the Atlantic area (annual average from ≈7°C in Norway up to ≈16°C off the Spanish coast), which increases further in the transition to the Mediterranean Sea (≈21°C), especially during summer. A sharp winter versus summer variation is also found within the inner and also when comparing these two regions with the Black Sea. Evidence of parallel evolution was detected between the two low salinity regions, the Baltic and Black seas. Correlation between genetic and environmental variation indicated that temperature and salinity were probably the main environmental drivers of selection.
Mining around the four genomic regions consistently inferred to be under selection identified candidate genes related to osmoregulation, growth, and resistance to diseases. The new insights are useful for the management of turbot fisheries and aquaculture by providing the baseline for evaluating the consequences of turbot releases from restocking and farming.
The uneven distribution of turbot has been associated with phylogeographic events related to rapid adaptive radiation following the last glaciation and to heterogeneous environmental conditions across its distribution range (Bouza, Presa, Castro, Sánchez, & Martínez, 2002;Vandamme et al., 2014). Some life-history traits, such as growth, survival, reproduction, and fecundity, have been shown to be influenced by temperature (Felix, Vinagre, & Cabral, 2011). Previous population genetic studies have shown low or no genetic population structure over large geographical areas, such as in the Northeast Atlantic Ocean. This has been attributed to the advection of larvae and, in some cases, also to the active migration of adults (Bouza, Sánchez, & Martínez, 1997;Bouza et al., 2002Bouza et al., , 2014Coughlan et al., 1998). Genetic divergence has been documented to be mainly associated with isolated areas in the Western and Eastern Mediterranean Sea (Atanassov et al., 2011). Low but significant differentiation between the Atlantic Ocean and Baltic Sea turbot has also been reported (Florin & Höglund, 2007;Nielsen et al., 2004;Vilas et al., 2010).
Despite the overall high genetic homogeneity recorded for turbot, strong site fidelity to the spawning sites and limited dispersal of adults during the reproduction period have been documented in the Baltic Sea (Florin & Franzén, 2010), suggesting that geographical segregation, even within continuous areas, might occur. Additionally, evidence suggestive of adaptation to temperature and salinity at a regional level has been reported Vilas et al., 2015). Reproductive success and growth differences in turbot between the Atlantic Ocean and the Baltic Sea have been associated with salinity (Nissling, Johansson, & Jacobsson, 2006) and have been explained either by phenotypic plasticity (Florin & Höglund, 2007) or divergent selection (Vilas et al., 2010(Vilas et al., , 2015. Significant divergence at specific markers (SNPs or microsatellites) has been reported between turbot sampled from the Irish shelf and the North Sea, the English Channel and the Biscay Gulf, between southern and northern North Sea, and between Baltic and North Sea Vilas et al., 2010Vilas et al., , 2015. An earlier study identified different hemoglobin genotypes, which suggested that turbot populations in the Northern Atlantic Ocean might not be entirely homogeneous (Imsland, Scanu, & Naevdal, 2003). These data highlight the need for more detailed studies using larger genomic coverage to clearly elucidate both neutral and adaptive genetic differentiation. Moreover, despite turbot being well-studied, its population structure has not been explored across its full distribution range; knowledge on wild populations is mostly limited to the Atlantic and Baltic areas and genomic features associated with environmental variables have only been recently investigated in a limited number of populations (Vilas et al., 2010(Vilas et al., , 2015 and markers . With the turbot genome recently sequenced (Figueras et al., 2016)  2017), it is now practical to consider conducting a genomewide scan of turbot populations, in order to better describe and understand its population genetic structure.
In this study, a panel of 755 SNPs evenly distributed over the turbot genome were genotyped using the double-digest RAD sequencing (ddRAD-seq) technology (Peterson, Weber, Kay, Fisher, & Hoekstra, 2012) and used to screen populations at a large geographical scale to (i) evaluate the level of population differentiation at small and large scales, (ii) test whether similar environmental conditions led to parallel evolution/adaptation in geographically distant/independent populations, and (iii) test the discrimination power of neutral versus outlier markers to define turbot populations for later applications in traceability studies. The information gathered is useful for a sustainable management of genetic resources in the wild and for guiding selection of genetic raw material for the growing turbot aquaculture.

| Sampling
A total of 697 individuals were collected at 20 sampling sites, mostly exceeding 20 individuals per sample, and often above 30 ( Figure 1,

| SNP calling and genotyping
A ddRAD genotyping-by-sequencing (GBS) approach was carried out following the procedure first described by Peterson et al. (2012) with slight modifications, detailed in full elsewhere (Brown et al., 2016  Effective population size was estimated by NeESTIMATOR v2.01 (Do et al., 2014) using the linkage disequilibrium (LD) method and loci with a MAF >0.02.

| Genetic diversity and structure
Pairwise F ST values between sampling sites and between geographical regions were estimated with ARLEQUIN v3.5 and tested for significance using 10,000 permutations. To further investigate population structure, a Bayesian individual clustering approach was applied using the program STRUCTURE v2.3.4 (Pritchard, Stephens, & Donnelly, 2000) based on the admixture ancestry model and correlated allele frequencies. A burn-in of 10,000 iterations was used, followed by a MCMC procedure (Markov chain Monte Carlo) of 50,000 repeats, and tested K (number of genetic clusters or units) from 1 to 10. Five independent runs were used for each K estimate when the full marker dataset was used, while 10 runs per K were tested for the neutral and outliers marker subsets. The locprior model, which specifies the population of origin of each individual, was also used, as previously suggested for data showing weak structure (Hubisz, Falush, Stephens, & Pritchard, 2009;Pritchard et al., 2000;Vandamme et al., 2014).
Results from STRUCTURE were processed with the program STRUCTURE HARVESTER v0.3 (Earl & vonHoldt, 2012) to estimate the best-fitted number of clusters K based on the ΔK method described by Evanno, Regnaut, and Goudet (2005). However, as hierarchical analysis assumed by STRUCTURE may overshadow other K-values and affect the detection of substructure, HARVESTER was also run excluding K = 1. Under this scenario and for all tests, several K assignments were explored to better resolve subtle structuring, as previously suggested (Pritchard et al., 2000;Vandamme et al., 2014).
Additionally, because STRUCTURE inference may be affected by uneven sampling (Puechmaille, 2016), as in our case due to the ex- For all statistical tests, significance level (p < .05) was corrected for multiple comparisons using the sequential Bonferroni method (Rice, 1989).

| Detection of outlier loci and mining of the turbot genome
To identify candidate loci subjected to selection, three different algorithms and related software packages were applied as previously suggested (Narum & Hess, 2011;Shimada, Shikano, & Merilä, 2011;Souche et al., 2015). BAYESCAN v2.01 (Foll & Gaggiotti, 2008 were considered as outlier loci depending on the analysis. Taking into account the low structuring of turbot across its distribution range (Bouza et al., 2014), the identification of outliers for divergent selection would be expected more feasible and confident than those for stabilizing selection, especially in the Atlantic Ocean and Baltic Sea, the most representative regions of the species.
Accordingly, all loci detected by any of the aforementioned methods were considered the initial set of candidates for divergent selection. A more stringent approach was used to identify divergent outliers with the highest confidence. Among the programs applied, BAYESCAN follows the most conservative approach and identifies a smaller number of markers than LOSITAN and ARLEQUIN, which usually return a high proportion of false positives (Narum & Hess, 2011;Shimada et al., 2011). Accordingly, all outlier loci detected with BAYESCAN (p < .050) were retained, and additionally, only those that were identified both with LOSITAN and ARLEQUIN at the highest confidence level (CI < 99 and FDR = 0.010 for LOSITAN and p < .010 for ARLEQUIN). Outlier loci for balancing selection are more difficult to detect, and a high rate of false positives is returned with most programs (Narum & Hess, 2011); this fact could be even stressed considering turbot structure. Thus, only markers detected with at least two of the three approaches were considered significant in our study, additionally considering the environmental scenario of the area analyzed (see Results).
Outliers were investigated using all samples or subgroups of samples according to the observed regional structure (see Results (2). The AD sample was not included in this analysis due to its limited representativeness of the Mediterranean area.

Sequences (90 bases) including the most consistent set of outlier
SNPs were mapped to the turbot genome to obtain their genomic position using local BLASTn (e-value <1e-20). These outliers were subsequently anchored to the turbot genetic map (linkage group (LG) and predicted position in cM) to investigate their putative relationship with previously reported outliers (Vilas et al., 2010(Vilas et al., , 2015 or QTLs . Gene mining around identified outlier loci was accomplished using the turbot genome browser (http://denovo.cnag.cat/genomes/turbot/; Figueras et al., 2016).
Gene lists were compiled along an interval of ±250 kb flanking the outlier locus position, this considered appropriate due to the low population linkage disequilibrium observed in turbot (Saura et al., 2017). Suggestive candidate genes involved in immunity and growth, or putatively related to environmental variables, along with those in the vicinity of previously described QTLs were retrieved from gene lists. Functional enrichment of the gene lists against the turbot genome was undertaken with BLAST2GO (Conesa et al., 2005).

| Seascape analyses
The environmental variables "sea surface temperature" and "sea bottom temperature" (SST and SBT, respectively, in °C), "sea surface salinity" and "sea bottom salinity" (SSS and SBS, respectively, in PSU), and "primary production" (PP in g C m −2 day −1 ) were retrieved from the Environmental Marine Information System (EMIS) database Genetic differentiation explained by spatial variables (i.e., latitude and longitude) and the registered environmental variables (SST, SBT, SSS, SBS, and PP) was studied by a canonical redundancy analysis (RDA) using the R platform VEGAN software (Oksanen, 2015). For each sample, allele frequencies were estimated for all analyzed loci with the ADEGENET package on the R platform using the "makefreq" option applied on the ADEGENET "genpop" file. Only two loci showed missing values to estimate allele frequencies, and consequently, they were removed for analyses as VEGAN does not allow missing information in their input files. ANOVA tests were performed to check the significance of the variance associated with the different environmental variables using 1,000 random permutations with VEGAN.
Variance inflation factor (VIF) was estimated to explore collinearity (correlation) among environmental variables in our dataset. VIF values <5 show no collinearity problems, values from 5 to 10 represent moderate problems, and values >10 indicate important collinearity problems (Marquardt, 1970). Different models were adjusted by automatic forward selection based on significant variable criteria. These analyses were performed using the PACKFOR package in R (Dray, Legendre, & Blanchet, 2009).
Forward selection corrects for highly inflated type I errors and overestimated amounts of explained variation . Thus, the reduced panel of explanatory variables was used to recalculate the total proportion of genetic variation in the variance partitioning. The weight of the different loci on the significant environmental vectors was obtained using VEGAN. All these analyses were performed separately for the full and the neutral SNP datasets.

| Sequencing
A total of 783,829,288 read-pairs were retrieved from the sequencing platform, with ~1 million read-pairs on average per sample. The distribution of reads per sample was rather unbalanced, ranging from 125,563 to 5,842,110, mostly attributed to the low DNA quality of some archived samples, along with the specific methodology employed; that is, samples were combined immediately following adapter ligation (Brown et al., 2016) and not amplified, size selected and re-quantified individually before pooling into a library as described for the original protocol (Peterson et al., 2012). After barcode and restriction enzyme site identification, an average of 85% of read-pairs passed the first quality-filtering step, leaving 659,968,515 sequences to feed into the STACKs pipeline. We developed a cus-

| Genetic structure
Moderate but highly significant genetic differentiation (F ST = 0.090, p < .001) was detected across the whole turbot distribution using the 755 SNP panel, mostly attributable to divergence between the more isolated areas (Black Sea and the Adriatic Sea) and the Atlantic Ocean and the Baltic Sea (Table 3 and Table S1). The most likely number of genetic units (K) was four as revealed by STRUCTURE (i.e., ATL, BAS, AD, BLS; Figure 2 and Figure S1) and DAPC ( Figure   S2a), corresponding with the four main geographical sampled areas ( Figure 1). The divergence was moderate-high between the southern (AD and BLS) and the Atlantic region (F ST = 0.137 and 0.134, respectively; p < .001), and low, but significant, between the Baltic and the Atlantic region (F ST = 0.005; p < .001; Table 3). All pairwise F ST (Table S1)
These outliers were mostly detected in BAS & BLS (7), two distant regions which share a low salinity environment, and an additional one was detected when comparing BAS populations (Table 4).
LG9: Outlier 1056_25 was located within an intron of MTMR7 (myotubularin-related protein 7), a gene related to lipid and protein metabolism (Safran et al., 2010). This gene is within the confidence interval of an A. salmonicida resistance QTL and tightly linked to a previously reported outlier (~ 63 kb; Sma-E117; Vilas et al., 2015).
Data mining revealed the presence of NFIL3 (nuclear factor, interleukin 3 regulated), a candidate gene for VHSV resistance (Figueras  LG, linkage group; Selection: type of selection and geographical region; Gene/Function: official gene symbol and function (IG: immunoglobulin); data mining: relevant genes identified related to growth (G), VHSV resistance (VH), Aeromonas salmonicida resistance (AS), reproduction (R) and osmoregulation (OS); QTLs and previous outliers from 1  and PATL1 (PAT1 homolog 1, processing body mRNA decay factor; Grady et al., 2014). The outlier 1916_69 was located in an intergenic region, relatively close to a previously reported SNP outlier (<3 cM; SmaSNP247; Vilas et al., 2015). Functional enrichment around both outlier loci showed genes related to metabolism and immunity (Table 5).
LG10: Outlier 5848_28 was located within an intron of HDAC7 resistance (Figueras et al., 2016), and TWIST2 (twist family BHLH transcription factor 2), a gene related to osmoregulation (Grady et al., 2014). Functional enrichment displayed GO terms related to primary metabolism and immunity (Table 5).

| Seascape analysis
All environmental variables showed VIF values below 5, suggesting no collinearity issues between them. Redundancy analyses including spatial and environmental variables (Table 1)

| Genetic structure based on neutral and outlier markers
Population structure was reanalyzed using the most consistent neutral (513 loci) and outlier (25 loci) SNP datasets to compare patterns of genetic differentiation related to demographic parameters and selection (Table 4). A third subpanel was also used to analyze differentiation in the ATL and BAS using the six outlier loci detected in these comparisons (16278_38, 1916_69, 1056_25, 16775_23, 6850_51, and 7550_55; Table 4).

Levels of genetic differentiation estimated for the neutral panel
were lower than estimated from the full SNP dataset (Table 3

| Genetic diversity and population structure
In this study, the population genetic structure of turbot was analyzed for the first time using a genomewide SNP panel (>1 SNP per Mb) developed with a genotyping-by-sequencing (GBS) strategy, namely ddRAD (Peterson et al., 2012). This allowed us to screen populations from the entire distribution range avoiding the ascertainment bias that can affect analyses with fixed SNP panels (e.g., SNP chips) developed from only a few populations. The combination of balanced multiplexing of individuals per library along with the progressive lowering cost of high-throughput sequencing made the approach cost-effective. A set of 755 highly confident SNPs (0.5% genotyping error) were obtained covering homogeneously the turbot genome at >1 SNP/Mb from the initial 5,564 SNPs. The proportion of filtered SNPs was in line with those previously reported in marine fish species (Palaiokostas et al., 2015) and the genotyping error was even lower (Mastretta-Yanes et al., 2015), supporting the robustness of our custom filtering pipeline.
Our analysis represents the largest effort to date to analyze the genetic diversity and structure of turbot across its entire distribution range including the Atlantic Ocean and the Mediterranean Sea, as well as the inner Baltic and Black seas. Our sample collection was rather uneven due to the overrepresentation of the Atlantic area (14 of 20 samples), which suggests some caution for statistical bias when analyzing the whole dataset. However, this sampling reflects quite well turbot distribution, which is common in the Atlantic Ocean and very scarce in the Mediterranean Sea, and therefore, the analysis should take this fact into consideration. Anyway, we split our analysis also by region and made all meaningful comparisons between regions to obtain the most comprehensive view of turbot structure.
Genetic diversity in the Atlantic Ocean (H E : 0.095) was much lower than previously reported (Vilas et al., 2015; H E : ~0.300), mainly due to the high number of samples used to identify polymorphic loci in our study (697) and the low MAF cutoff used (0.002). This resulted in the detection of a large fraction of loci with very rare alleles compared to previous studies (sample size ~30; Vera et al., 2013;Vilas et al., 2015). In fact, genetic diversity estimated using the fraction of loci at P 95 almost doubled (H E : 0.184) approaching to previous results (Vilas et al., 2015). Anyway, our data suggest that genetic di- Given the general lack of physical barriers in the sea, marine fish, such as turbot, with wide distribution ranges, high fecundity, pelagic eggs, and larvae, are expected to be at best weakly structured across large areas. The only study exploring genetic structure in turbot covering a roughly similar distribution range was based on allozyme marker data (Blanquer et al., 1992 Very low genetic differentiation between samples was observed within the Atlantic region, supporting previous findings (Bouza et al., 2002;Coughlan et al., 1998;Florin & Höglund, 2007;Nielsen et al., 2004;Vandamme et al., 2014;Vilas et al., 2015), which indicates relatively high levels of gene flow. This is the case notwithstanding the presence of different well-known current fronts inside the large marine ecosystems (LME) surveyed in this study, such as the Iberian Coastal, Irish and Biscay Shelf, and North Sea (Belkin & Cornillon, 2007;Belkin, Cornillon, & Sherman, 2009;Vandamme et al., 2014).
Within this relatively homogeneous genetic scenario, subtle, but rather consistent differences along with low N e estimates were detected at both geographical extremes, that is, Norwegian and Spanish Coast samples. Vilas et al. (2015) also suggested that Iberian turbot should be considered genetically distinct from elsewhere in the Atlantic region. Restricted gene flow to more northern regions due to oceanic fronts off the Galician coast has been suggested for other species with passive larval dispersal such as the flat oyster Ostrea edulis (Vera et al., 2016a).
Low but significant genetic differentiation was detected between the Atlantic and Baltic Sea turbot (F ST = 0.005; p < .001), a substructure previously reported based on other SNP and microsatellite genotypes (Nielsen et al., 2004;Vandamme et al., 2014;Vilas et al., 2010Vilas et al., , 2015 and, in part, attributed to the biogeographic barrier occurring between the North Sea and the Baltic Sea (Johannesson & André, 2006). Marginal or geographically isolated populations are often more prone to the effects of genetic drift and show higher genetic divergence and lower diversity than those closer to the center of the species distribution (Kawecki, 2008;Lira-Noriega & Manthey, 2014). Other marine fish distributed in the Baltic Sea also show lower genetic diversity than conspecific Atlantic populations due to varying processes of isolation over 4,000-8,000 years since colonization after the last glaciation (Littorina period; Johannesson & André, 2006). The more pronounced deviations from HWE observed in the Baltic samples may suggest a Wahlund effect, the Baltic Sea being a more heterogeneous environment and having a relatively recent history of colonization from the Atlantic, constituting a partial transitional environment (Nielsen et al., 2004).

| Local adaptation in turbot
The detection of signals for divergent selection in genomes is favored in scenarios of relatively large N e , as in turbot, because the genetic signal left by the demographic history will be easier to erode and the ability to detect high differentiation outliers is favored by a low baseline level of neutral genetic differentiation between populations. Although turbot populations generally exhibited low to moderate genetic structure, populations from the Black Sea and the Adriatic Sea showed evidence of geographical isolation. The discrimination between neutral and selective variation in the Baltic-Atlantic transition zone may be further complicated as correlations between genetic and environmental variation might be due to reasons other than natural selection (Bierne, Roze, & Welch, 2013;Bierne et al., 2013). Therefore, outlier loci were carefully assessed in our study and were only considered reliable when they showed strong statistical support and previous phenotypic information consistent with predefined hypotheses. Genomic co-localization with previously reported outlier loci or association with growth or disease-resistance-related QTL evaluated so far in turbot  was also considered, as they might point out to genomic islands of divergence, as reported in other fish (Bradbury et al., 2013).
Evidence of divergent selection in turbot was mainly detected in the comparisons between the Atlantic region and those regions with low salinity (BAS and BLS). Differentiation due to selection (and genetic drift) may be favored by limited gene flow related to differences in salinity tolerance. It is known that Atlantic turbot eggs do not survive at the lower salinities of the Baltic Sea (Florin & Höglund, 2007), and, in addition, because turbot eggs are not buoyant at salinities below 20 PSU, eggs from the Baltic are demersal rather than pelagic (Nissling et al., 2006). Three of the five markers that were statistically significant in the ATL-BAS comparison (1916_69 at LG9, 6850_51 and 7550_55 at LG2) were also significant in the comparison between ATL and BLS, strongly supporting that their divergence might be related to adaptation to differences in salinity. Furthermore, the marker 1916_69 is closely linked to a previously reported outlier (SmaSNP247; Vilas et al., 2015), and both of these markers show a pattern consistent with divergent selection and are located relatively close (3 cM distant; Figueras et al., 2016).
Although 1056_25, also at LG9, was only divergent in the comparison with BAS, it was located within an important functional region including genes related to osmoregulation and is tightly linked to a previously reported divergent outlier (SmaE117) between the Atlantic Ocean and the Baltic Sea (Vilas et al., 2015). These results suggest that several loci within a narrow region at LG9 show a spatial pattern of structuring consistent with adaptation to environments that differ markedly in salinity.
Marker 7574_88 selectively diverged in the ATL-BLS comparison with all three analytical tests employed, but not in the ATL-BAS comparison. Most important, this outlier is located in a genomic region related to growth (Figueras et al., 2016;Rodríguez-Ramilo et al., 2014), and further, it showed a gradual cline from Baltic through to Black Sea samples. A similar pattern was detected for 2372_27, a linked marker at LG1 (~400 kb distant from 7574_88), which was significant for divergent selection with LOSITAN and ARLEQUIN.
Interestingly, the outlier 5397_68, located very close to 2372_27 at LG1 (<100 kb distant), showed signals of stabilizing selection in the BAS-BLS comparison. The marker SNP 5397_68 is located in the vicinity (~100-400 kb distant) of genes related to growth (e.g., GPAA1, EXOSC4, PRKACB), and further, Norman, Ferguson, and Danzmann (2014) detected several QTLs related to salinity tolerance in an orthologous region in Arctic charr (Salvelinus alpinus, LG32), which contains three genes related to osmoregulation (EFID, PPM1L and UCK2). This genomic region at LG1, which includes several outlier loci, growth and VHS resistance QTLs, and growth and salinity tolerance candidate genes, seems relevant to explain the genetic structure of turbot.
Although still a controversial topic, balancing selection (or parallel evolution) has been identified as an important factor in the evolution of some marine species, for example, related to coral reef fishes (Gaither et al., 2015), European sea bass (Lemaire et al., 2000), and three-spined stickleback (Feulner et al., 2013;Guo, DeFaveri, Sotelo, Nair, & Merilä, 2015). Moreover, it has also been suggested to play an important role in invasive processes of aquatic organisms (Vera, Díez-del-Molino, & García-Marín, 2016b). Adaptive variants maintained by balancing selection have been reported in different species, such as those in immune-related genes in three-spined stickleback (Feulner et al., 2013;Guo et al., 2015) and in ribosomal structure and regulation genes in the albacore tuna (Laconcha et al., 2015). Seven of the eight outlier loci identified as being influenced by balancing selection in the current study were detected in the BAS-BLS comparison, a result consistent with parallel adaptation to low salinities.
In summary, our data not only confirm previous genetic diversity and population structure data based on different markers, but reveals crucial novel information on population adaptation and connectivity using a combination of neutral and adaptive genetic variation. Low but significant differentiation was detected between the Atlantic and Baltic regions, while high differentiation was observed between the Atlantic and the southeastern-most regions (Adriatic and Black seas), indicating that demographic and historical factors have contributed to shaping population structure of turbot across its natural distribution. The information reported here also provides for the first time new insights on turbot adaptation, especially in the Black Sea, and suggested parallel evolution between areas with similar environmental conditions. Both strong neutral evolutionary forces and adaptive selection appear to be acting simultaneously on geographically isolated populations in the natural distribution of turbot. However, subtle neutral differentiation and local adaptations might also be occurring within regions. Candidate outlier loci, mostly anchored to the turbot genome, and especially at specific regions of LG1 and LG9, showed a positive correlation with environmental variables related to salinity and temperature.

| Management implications
Our results represent useful information for the management of wild stocks, and they can be valuable for breeding programs of farmed turbot. An improved definition of management units considering both demography and adaptation to environmental variation along the whole distribution range can now be delineated, allowing the future definition of adaptive management units (AMU, Bernatchez et al., 2017). Four main operational units can be defined related to the four main genetic regions identified along the distribution range, but further refinement should be considered within the Atlantic area, where both Norway and Spanish samples showed slight differentiation from the Atlantic core both using all data and outlier loci. Our study did not detect significant subdivision in the British Isles and the North Sea as previously suggested by Vandamme et al. (2014). Further, although only two samples were analyzed, the North and South Baltic Sea showed significant differences, both when considering the full dataset, and outlier and neutral loci separately. In addition, our data represent the baseline to monitor restocking performed in the Atlantic area with unknown consequences, and when the genetic composition of the broodstock of the main turbot farms be at hand, to evaluate the impact and introgression from farm escapes to evolve toward a sustainable aquaculture. Finally, providing breeders with information of natural resources regarding environmental variables will be highly useful to boost breeding programs fitting them to market demands, to found broodstock at farms in new geographical regions, and to face the new challenges of climatic change.

ACK N OWLED G EM ENTS
The authors thank all Aquatrace partners for sampling support and Data available from the Dryad Digital Repository: https://doi. org/10.5061/dryad.pq3376q.

CO N FLI C T O F I NTE R E S T
None declared.