Neutral and adaptive loci reveal fine‐scale population structure in Eleginops maclovinus from north Patagonia

Abstract Patagonia is an understudied area, especially when it comes to population genomic studies with relevance to fishery management. However, the dynamic and heterogeneous landscape in this area can harbor an important but cryptic genetic population structure. Once such information is revealed, it can be integrated into the management of infrequently investigated species. Eleginops maclovinus is a protandrous hermaphrodite species with economic importance for local communities that are currently managed as a single genetic unit. In this study, we sampled five locations distributed across a salinity cline from Northern Patagonia to investigate the genetic population structure of E. maclovinus. We used restriction site‐associated DNA (RAD) sequencing and outlier tests to obtain neutral and adaptive loci, using FST and GEA approaches. We identified a spatial pattern of structuration with gene flow and spatial selection by environmental association. Neutral and adaptive loci showed two and three genetic groups, respectively. The effective population sizes estimated ranged from 572 (Chepu) to 14,454 (Chaitén) and were influenced more by locality than by salinity cline. We found loci putatively associated with salinity suggesting that salinity may act as a selective driver in E. maclovinus populations. These results suggest a complex interaction between genetic drift, gene flow, and natural selection in this area. Our findings also suggest several evolutionary significant units in this area, and the information should be integrated into the management of this species. We discussed the significance of these results for fishery management and suggest future directions to improve our understanding of how E. maclovinus has adapted to the dynamic waters of Northern Patagonia.


| INTRODUC TI ON
Advances in genome sampling methods have reduced complexity (e.g., Restriction associated DNA sequencing, RADseq) and allowed the collection of an unusual amount of data to analyze the genome concerning conservation and management problems Xuereb et al., 2021). These data provide a great way to solve unanswered questions and have the advantage of allowing the quantification of adaptive variation unlike microsatellites Funk et al., 2012). Currently, researchers can differentiate neutral and adaptive variation across population information that can be incorporated into management and conservation programs to obtain better solutions on them (Xuereb et al., 2021).

Identification of evolutionary significant units (ESU) is important
to guide management and conservation efforts (Funk et al., 2012) and maximize the evolutionary potential for environmental change (Bernatchez, 2016). At a fine scale, a management unit, which is included in an ESU, refers to demographically independent populations and shows significant divergence results with low gene flow (Funk et al., 2012;Moritz, 1994;von der Heyden, 2017). Finding genetic differences at a fine scale is challenging because it depends on the biological characteristics of the organisms studied (e.g., vagility) and the geomorphological conformation and environmental heterogeneity of their geographical distribution (Canales-Aguirre et al., 2016;Jørgensen et al., 2005). For the latter, habitats such as fjords may greatly affect population genetic diversity in marine organisms due to unique environmental characteristics.
Patagonia in Chile includes a vast coastal area (240,000 km 2 ; Pantoja et al., 2011). The northern region extends from latitude 41.5°S (Reloncaví Fjord) to latitude 46.5°S (San Rafael Lagune; Rodrigo, 2008) with high ecosystem productivity and heterogeneous geomorphological and physical-chemical oceanographic conditions (Pérez-Santos et al., 2014;Ríos et al., 2016;Yevenes et al., 2017). For example, Patagonia has a saline cline pattern resulting from freshwater runoff from melting ice from the Andean Mountains (annual average caudal >300 m 3 /s) and an annual average precipitation >1000 mm (Garreaud et al., 2013). Salinity can play a role in the ontogeny resulting in a differentiated vertical distribution of eggs and yolk sac larvae (Petereit et al., 2009), and also, it influences population abundance as a result of freshwater discharges (Ojaveer & Kalejs, 2010). Also, a heterogeneous-salinity environment can promote local adaptation in marine populations that can result in genetic population structure differences (Berg et al., 2015;Limborg et al., 2009;McCairns & Bernatchez, 2008). Such landscape characteristics influence a high diversity of marine organisms and hierarchical levels, from populations to ecosystems (Beuchel et al., 2006;Canales-Aguirre et al., 2010, 2016Kristoffersen & Salvanes, 2009;Olsen et al., 2002), that support and sustain economically important fisheries. Unfortunately, this area has been understudied, and no population genomics studies with relevance to fisheries management have been conducted despite the fact that this unique landscape may promote large genetic differentiation and local adaptation. Genomic studies can reveal discrete genetic groups that can be integrated into future conservation and management programs (e.g., Euclide et al., 2021;Larson, Seeb, Pascal, Templin, & Seeb, 2014;.

Eleginops maclovinus (Cuvier and Valenciennes 1830) is a mono-
typic species and one of the few species of the Notothenioidei with a non-Antarctic distribution (Bargelloni et al., 2000;Matschiner et al., 2015;Near, 2004). Endemic to South America, it is found close to estuaries from ~33°S on both sides of the Pacific and the Atlantic Ocean to the Beagle Channel (~54°S) (Pequeño, 1989), including in the Malvinas/Falkland Islands (Gosztonyi, 1974). This species is economically important for local communities and is caught by artisanal and recreational fishers (Gastaldi et al., 2009), with 182 tons landed in 2019 in Chile (Sernapesca, 2019). In Chilean Patagonia, there is little scientific information about this fish that can be used to identify appropriate management units. Eleginops maclovinus is a partial spawner with a spawning peak during late autumn in estuaries of the Chilean coast (Ruiz, 1993); it is also the most fecund of the Notothenioidei (~550,000 eggs per female, Gosztonyi, 1979). This species is a protandrous hermaphrodite, and males are often smaller in length than females (10-52 cm males and >53 cm females) (Brickle et al., 2005).
Connectivity and dispersal between populations of E. maclovinus seem to be biased by sex. Mechanical tags from females have been found up to 60 nautical miles away from their tagging location (Brickle et al., 2003), while parasites used as biological tags suggest that males are residents that do not move large distances (Brickle & MacKenzie, 2007). This information supports the idea that E. maclovinus has a structured geographic distribution; however, this has not been supported by research based on molecular approaches.
For instance, phylogeographic analysis using mtDNA, which only should reflect female dispersal, showed that E. maclovinus had weak genetic differentiation, shared haplotypes among locations, and recent population expansion that would have occurred as a result of Quaternary Glaciations (Ceballos et al., 2012). In addition, a subsequent study using microsatellites revealed a low but significant level of genetic differentiation between Pacific and Atlantic populations, with Atlantic populations showing a mixed membership from the two main genetic clusters. Some degree of genetic heterogeneity was suggested within the Pacific at a lower hierarchical level (Ceballos et al., 2016). The microsatellites study also suggested that northern populations, in both Atlantic and Pacific Oceans, may harbor more genetic variability as revealed by the number of private alleles (Ceballos et al., 2016). This suggests that Northern Patagonia fjords may harbor a cryptic population genetic structure at a fine scale, as a result of the geomorphological and heterogenous landscape, which was not observed probably due to lack of sampling locations from this region or low resolution of the molecular information used.
In this study, we investigate the population genomics of E. maclovinus from Northern Patagonia, using RADseq to assess both neutral and adaptive variation. We specifically aim to (i) evaluate the extent of neutral and adaptive genetic diversity, population differentiation between locations, and estimation of effective population size (ii) correlate putative adaptive loci to environmental variables, (iii) identify putative functions for candidate loci, and (iv) discuss the implications of our results for conservation and management of the species.

| RADseq library and genotyping
We obtained high-quality DNA using two steps. First, we used a traditional phenol-chloroform DNA extraction protocol to obtain a large amount of genomic DNA; second, we purified the DNA using DNeasy Blood & Tissue Kit (Qiagen®) following the manufacturer's instructions but skipping the lysis step. We quantified the total double-stranded DNA using a Qubit® 3.0 Fluorometer (Invitrogen) and the dsDNA BR Assay Kit, following the manufacturer's instructions. Each DNA sample was dried-down and then normalized with sterile water to 20 ng/μl in a final volume of 10 μl (i.e., 200 ng per individual).
Each DNA sample was digested using the SbfI-HF® restriction enzyme (New England Biolabs) and ligated with an 8 bp unique barcode adaptor. Barcoded individual DNA samples were pooled into master libraries and fragmented to ~300-500 bp with 12-14 30 s cycles in a Q500 sonicator (Qsonica). The fragmented DNA library was bound to Dynabeads™ M-280 Streptavidin magnetic beads (Invitrogen). Subsequently, nontarget fragments were removed by washing with a TLE buffer, and DNA was released from the beads using an incubation step. We conducted three purification steps for each library: (i) nonligated adaptors or small fragments of DNA, and enzymes for each library were removed using AMPure XP beads (Beckman Coulter); (ii) all fragments of DNA that were not barcoded were removed using Dynabeads™ M-280 Streptavidin magnetic beads leaving only the DNA fragments with the restriction site and barcode; (iii) any residuals of small fragments of DNA and all enzymes and other impurities were also discarded using AMPure XP beads. NEBNext® Ultra™ DNA Library Prep Kit for Illumina® was used for ligation of master library barcodes, a 250bp insert size selection, and a 12-cycle PCR enrichment. We confirmed enrichment and size selection by visualizing PCR products on a 2% agarose E-Gel (Invitrogen). A final AMPure XP purification clean-up step was conducted prior to quantifying the DNA with a Qubit® 2.0 Fluorometer. All prepared libraries for paired-end libraries were sent to Novogene for sequencing on the Illumina NovaseqS4 platform.

| SNP discovery
To discover and genotype SNPs from the raw RADseq data, we used a similar pipeline to the one suggested by Rochette and Catchen (2017) using the software STACKS v2.41 (Rochette et al., 2019). We selected a subset of 60 samples for testing and to select an optimal set of parameters as suggested by Paris et al. (2017). Raw sequences were demultiplexed by barcode using PROCESS_RADTAGS program, where intact barcode and cut-site of SbfI restriction enzyme were checked, and whole reads with quality issues were discarded and trimmed to 140 bp (parameter flags: -e SbfI -c -q -r -t 140 --filter_illumina --bestrad). Individuals that showed <950,000 reads retained were excluded from the downstream STACKS pipeline. We kept a total of 112 individuals. Reads were then assembled to build unique stacks to identify putative loci using a maximum likelihood framework (Hohenlohe

| Bioinformatics and genotyping quality filters
We used the iterative filtering process described by McKinney et al. (2020), which applies a soft and then a more stringent threshold to remove poor-quality loci and samples. Filtering order was: (i) minor allele frequencies, (ii) genotype rate for loci, and (iii) genotype rate for sample, after which we recalculated the proportion of missing data, and ran steps (i), (ii), and (iii) with more stringent thresholds. For minor allele frequencies (MAF), we start by removing loci with MAF ≤ 0.05, and then, MAF ≤ 0.1. For the genotype rate for loci, the first threshold was set at 25% and the second at 90%. For the genotype rate for samples, the first threshold was set at 50% and the second at 85%. Loci in Hardy-Weinberg disequilibrium (p < .05) were removed if they deviated in three or more populations. For tags with more than one SNP, we kept the putative SNP with the highest F ST to reduce the influence of linked loci in the results. Additionally, paralog sequence variants were identified using HDPLOT (McKinney, Waples, et al., 2017) and then removed (parameter flags: H < 0.6; |D| < 5). Paralogs result from gene duplication events that have affected the evolution of the notothenioid genome (Chen et al., 2008(Chen et al., , 2019. Paralogs are generally difficult to genotype reliably with RADseq data due to insufficient read depth (McKinney et al., 2018).
All variants that met with these criteria were retained for further analyses ( Table 1).

| Data analyses: Identifying neutral and adaptive loci
We used three methods to detect loci under selection (PCADAPT, FSTHET, and Redundancy Analysis). The first method inferred outliers based on principal component analysis (PCA) and was implemented using the PCADAPT package v4.3.3 (Privé et al., 2020).
This method assumes that markers excessively related to population structure are candidates for local adaptation (Luu et al., 2017).
The PCADAPT method uses the Mahalanobis distance (D) statistic, where a vector of the z-score is derived for regressing each SNP with K principal components Luu et al. (2017). We applied Cattell's rule to choose the K number of the principal components (Cattell, 1966). The p-values were obtained from transforming Mahalanobis distance (D) based on the chi-square distribution (Cattell, 1966). To avoid confounding effects of the population structure, we identified an optimal K-value testing from K = 1 to K = 10, and we checked the proportion of variance explained by each principal component using a scree plot using the pcadapt function and PCA ( Figure S1). K = 3 was retained, and we calculated the false discovery rate of the p-values associated with Mahalanobis distances using the QVALUE package (Dabney et al., 2010). A list of putative adaptive loci was obtained under an expected false discovery rate of α = 0.1. The second method corresponds to FSTHET, which identifies candidate loci by calculating smoothed quantiles between loci with strong differentiation F ST relative to their expected heterozygosity (Flanagan & Jones, 2017). This approach does not require any assumptions about the underlying population structure and is therefore more broadly applicable than other outlier detection methods. We calculated the empirical F ST based on Wright (1943) and expected heterozygosity. Loci were binned based on their expected heterozygosity values, sorted by F ST value, and quantiles were calculated. Loci that showed departures from a 95% confidence interval were considered under positive or balancing selection whether they surpassed superior or inferior confidence intervals, respectively ( Figure S2). Finally, we used a Genotype-Environment Association approach (GEA) conducting a Redundancy analysis (RDA) to identify signatures of local adaptation to environmental variables. RDA is a multivariate ordination method that combines PCs from allele frequency and multivariate environmental distance matrices to produce canonical axes predicting relationships between environments and particular loci (Forester et al., 2018;Rellstab et al., 2015). We used the Bio-Oracle dataset (Assis et al., 2018;Tyberghein et al., 2012) to obtain environmental variables for salinity, temperature, pH, oxygen, silicate, current velocity, primary production, phosphate, phytoplankton, iron, nitrate, chlorophyll a, and calcite (https:// www.bio-oracle.org/). For all these variables, we obtained values for the maximum, minimum, mean, and range of each variable when possible. RDA was performed using the VEGAN v. 2.3.4 R package (Oksanen et al., 2015). Variance inflation factor (VIF; vif. cca function of VEGAN) was used to ascertain the lack of multicollinearity among variables (Hair et al., 1995;Zuur et al., 2010) and excluded variables with a VIF ≥ 10 (Hair et al., 1995). Outliers were identified on each of the first three ordination axes as SNPs with a "locus score" that was ±3 SD from the mean score for that axis RDA, as suggested by Forester et al. (2018) to minimize false-positive and false-negative results. We then determined the correlation between each candidate SNP and one or more environmental variables.
The VennDiagram package (Chen & Boutros, 2011) was used to identify all putative adaptive loci identified and shared among the three software (PCADAPT, FSTHET, and RDA). Three datasets were built after identifying putative loci under selection: (i) neutral, (ii) adaptive loci merged, and (iii) adaptive loci shared. The neutral dataset included all loci that were not included in the adaptive loci merged (12,026 SNPs). The adaptive loci merged dataset included all loci that were identified as an outlier in each software (356 SNPs).
The adaptive loci shared dataset included all loci that were shared among three software (13 SNPs). These three datasets were used for further population genomic analyses.

| Summary statistics, population divergence, and effective population size
The summary statistics of genetic diversity expected heterozygosity (H E ) and observed heterozygosity (H O ) for each subset were calculated by location and conducted using the HIERFSTAT v0.04-10 package (Goudet, 2005). The number of polymorphic loci and the effective population size (N e ) of each location were obtained in NEESTIMATOR v2 (Do et al., 2014). The N e was estimated only for the neutral dataset using the LD method (Waples, 2006) updated for missing data and following Peel et al. (2013). Values of N e within corresponding 95% confidence intervals for each population were estimated using the following parameters: a minimum allele frequency cutoff of 0.01 and a random mating model.
We estimated the individual ancestry coefficients based on sparse non-negative matrix factorization algorithms (sNMF) using the package "LEA" in R (Frichot & François, 2015). In this package, we tested each dataset to reveal the population genetic structure. We identified the best number of genetic clusters (K) based on cross-validation and on an information theoretic measure, the cross-entropy criterion (Alexander & Lange, 2011;Frichot & François, 2015). We iteratively tested from K = 1 to K = 10, with 10 replicates, and with 10,000 permutations per K using the function obj.snmf in LEA. We conducted the statistical procedure PCA to reduce the multivariate SNP multilocus data into two orthogonal axes using the ADEGENET v2.0 package (Jombart, 2008;Jombart & Ahmed, 2011). We used the PCA approach in all subsets obtained. We used the PCA to seek a summary of the genetic diversity among the sampled individuals ignoring the assumptions of the Hardy-Weinberg equilibrium and linkage disequilibrium, which are often required in other individual-based models. Finally, we calculated pairwise F ST values for populations and performed significance tests for pairwise using 10,000 permutations in the STAMPP package (Pembleton et al., 2013). The F ST estimation was following the Wright (1949) method but was corrected by the unequal population size as updated by Weir and Cockerham (1984) (see Pembleton et al., 2013). This analysis was completed for each dataset. To understand the process that led to the population structure, we calculated directional migration rates among locations using the divMigrate function (Sundqvist et al., 2016) included in the diveRsity package (Keenan et al., 2013). Using all datasets, we explored the migration rates using the effective number of migrants index (Nm; Alcala et al., 2014), as a statistic to calculate relative migration. We used 1000 bootstraps to test the statistical significance of directional migration, and no filter threshold was applied to see all migration rates estimated. The conception behind the divMigrate function is that for each pair of populations, a hypothetical pool of migrants is created using the allelic frequencies inferred from the two populations compared. Then, a measure of genetic differentiation is estimated for the hypothetical pool and between each pair of population. This directional genetic differentiation obtained is then used to calculate the relative migration between the two populations. Additionally, a The neutral dataset included all loci that were not included in the adaptive loci merged. c The adaptive loci merged dataset included all loci that were identified as an outlier in each software. d The adaptive loci shared dataset included all loci that were shared among three software (PCADAPT, FSTHET, and RDA).
one-way anova was performed to evaluate whether the Nm was different for the three different datasets. This analysis was conducted using the R package rstatix v0.6.

| Putative function from Blast2Go
We conducted loci annotation to identify a putative function for candidate variants underlying positive selection obtained from PCADAPT, FSTHET, and RDA. We used the software Blast2Go included in OmicsBox following the annotation pipeline described by Götz et al. (2008). Briefly, we compared our candidate loci against the NCBI genomic database translating the sequences from nucleotide to protein using BLASTX. Then, we mapped homologous sequences to Gene Ontology (GO) terms. Finally, sequences were annotated by applying the Blast2GO annotation rule (see Götz et al., 2008). We tabulated the Tag_SNP from our read; a method that identified the outlier, environmental variable associated, Gene name, Gene Ontology ID, and their respective GO names. to current velocity ( Figure 2b). Finally, a total of 356 adaptive loci were identified among PCADAPT, FSTHET, and RDA, while13 loci were shared among the three software. Four loci were unique for PCADAPT, 267 for FSTHET, and 18 for RDA ( Figure S3).

| Summary statistics, population divergence, and effective population size
Summary statistics of genetic diversity revealed similar values for each location within each dataset (   Table 2). Only the Chaitén confidence interval included an infinite value.
The lowest value of the cross-entropy suggested K = 1 for neutral, K = 3 for adaptive merged, and K = 6 for adaptive shared dataset ( Figures S4-S6). We surveyed the admixture result from K = 2 to K = 5 and decided to keep K = 2 for neutral, K = 3 for adaptive merged, and K = 2 for the adaptive shared dataset (Figure 3a- (Table 3). Overall, relative migration rates using neutral dataset (Figure 5a p < .0001), where neutral loci showed a higher relative migration rate comparing with adaptive loci (Figure 5d).

| Putative function from Blast2GO
From 356 putative candidate loci (FSTHET, PCADAPT, and RDA), 124 loci were blasted in BlastX, 111 loci mapped were of homologue sequences to GO terms, and finally, 98 loci were annotated to GO terms. Based on these results, we found a variety of candidate genes whose functions involved mitotic cytokinesis, epithelial cell differentiation, embryonic morphogenesis, chondrocyte differentiation, and kidney development, among others (Table S1).

| DISCUSS ION
Under the assumption that heterogeneous landscape, low vagil-

ity, and biological conditions can result in population divergence in
Eleginops maclovinus, we aimed to (i) disentangle the differences in neutral and adaptive genetic variation, (ii) correlate putative adaptive loci to environmental variables, and (iii) identify putative functions for candidate loci. PCA and membership analyses revealed two (neutral loci) and three (adaptive loci) clusters, none of them previously described for this species. Neutral loci suggest a spatial pattern of structuration with gene flow, while adaptive loci suggest spatial selection by environmental association. We identify candidate loci for divergent selection mainly associated with biological processes (e.g., DNA repair, sodium ion transmembrane transport, and metabolism).
Contrasting N e estimations among populations were found, show-

| Neutral genetic variation
We found fine-scale genetic spatial pattern of structuration with gene flow based on neutral genomic data. We identified two ge- Pacific and Atlantic distribution. They suggested current and historical connectivity between populations from their expansion from the middle Pleistocene. Later, using microsatellites, low but significant regional differentiation between the Pacific Ocean and Atlantic Ocean locations was found (Ceballos et al., 2016). Similar studies using microsatellites did not find population structure in the Patagonian area from the Pacific Ocean (Canales-Aguirre et al., 2010. Our findings differ from at scale level with previous results and reinforce that genotyping tools such as RADseq increase the power to resolve shallow population structure in fish with gene flow (Hollenbeck et al., 2019;Larson, Seeb, Everett, et al., 2014;Luikart et al., 2003) when other lower resolution genomic approaches cannot.
Dispersal and reproductive attributes could contribute to the observed population pattern. Theoretically, hermaphroditic species tend to have more structured populations than gonochoristic species (Chopelet et al., 2009;Coscia et al., 2016), but information supporting this hypothesis is scarce. For example, Chopelet et al. (2009) conducted a meta-analysis testing this hypothesis and found no supporting evidence. They suggested that dispersal capacities and environmental barriers can play an underlying role in the variance of the genetic structuring of marine fish populations. Dispersal behavior in E. maclovinus has been recorded using biological (parasites) and mechanical tags (Brickle et al., 2003;Brickle & MacKenzie, 2007). Both studies indicate that juveniles (mainly males) tend to be residents and larger fish (mainly female) can migrate, comparatively, large distances (up to 60 nautical miles, i.e., 111 km) (Brickle et al., 2003;Brickle & MacKenzie, 2007). In our case, the maximum geographic distance between sampling locations is approximately 230 km (REL-CHE). It seems large enough to avoid connectivity but that is in a scenario where sampling locations are unique populations in the sampling area where E. maclovinus inhabits, which is not true. In is hard to suggest how drifting of early stages can connect populations through gene flow as has been suggested for several marine species (Benestan et al., 2021). Overall, this result indicates that the life-history trait is not enough to promote population structure in E.
maclovinus, but this is maintained by a complex interaction between migration and selection (see Adaptive variation section).
The geomorphology of the Northern Patagonia and the pattern of sea currents provide further evidence for the differences observed. The presence of channels, estuaries, close sounds, and fjords can be an efficient barrier for dispersal at different lifehistory stages. For example, it has been suggested that genetic isolation is related to the shallow sill depth and life-history behavior in Benthosema glaciale (Kristoffersen & Salvanes, 2009). The Reloncaví Estuary is 55 km long, with geomorphological and environmental characteristics that make this a unique area (Castillo et al., 2016(Castillo et al., , 2017. For instance, superficial currents resulting from wind, flowsdown to the Reloncaví Sound in winter and flows-up in spring and summer (Castillo et al., 2016(Castillo et al., , 2017. The latter could promote a safe environment for the early stages of marine organisms. Our samples were collected in the head of the estuary, a sheltered area, less influenced by sea tidal change and more by freshwater runoff. The sea current in Reloncaví Estuary mainly goes from the head to the mouth of the estuary (Castillo et al., 2016(Castillo et al., , 2017; however, it seems that some individuals are more similar to Hornopirén and Chaitén suggesting dispersal based on PCA. This is supported by the admixture of some individuals in the membership analyses and the cluster distribution by location. This singular geomorphology could promote a similar pattern of genetic diversity for other species. Unfortunately, to date, there are no genetic studies in fish or other marine taxa using genomic analyses for this specific microgeographical area.

| Adaptive variation
We identified two and three groups of samples in the shared and merged adaptive datasets, respectively. The 13 loci for the shared and 356 loci for the merged adaptive datasets give some clues about their effect on structuration patterns. We suspected that a small amount of adaptive SNPs could make a difference in the structuration pattern. These SNPs showed only a small fraction of the genome, leaving genomic regions not explored that may have likely loci under selective pressures. These 13 loci only reflect adaptive differences found for REL but not for CHE such as merged adaptive dataset shown. The underlying idea to highlight the results of both, shared and merged adaptive datasets was to be aware of outcomes minimizing the type I error and type II error, respectively; thus, do not miss loci under weak selection (Hoban et al., 2016;Lotterhos & Whitlock, 2015). Different statistical approaches to find outlier loci can identify different loci under selection (Lotterhos & Whitlock, 2015). What approach is better will depend on the demographic history of the population surveyed TA B L E 3 Pairwise Note: Pairwise Fst index was calculated using 12,026 SNPs for neutral loci, 356 SNPs for merged adaptive loci, and 13 SNPs for shared adaptive loci. F ST estimation was using Wright (1949) methods but corrected by Weir and Cockerham (1984) for uneven population size (see Pembleton et al., 2013).
Abbreviations: CHE, Chepu; CHA, Chaitén; HOR, Hornopirén; MAN, Manao; REL, Reloncaví Estuary. (Lotterhos & Whitlock, 2015). In our case, we do not know the real demographic history of E. maclovinus, therefore using all approaches and including all putative adaptive loci give us a broad view of how selective pressure is acting. Both adaptive datasets reinforce the idea that REL is different from the rest of the locations, but only the merged adaptive loci permit us to separate also CHE. This latter pattern fit well with the different seascape in the studied area (i.e., REL, influenced by runoff of fresh water; MAN-HOR-CHA, influenced by a mix of runoff fresh and marine waters; CHE, influenced by marine water).
The GEA analysis provided evidence for local adaptation to current velocity, primary productivity, oxygen, and salinity; regardless of the database used. Those environmental variables seem to be exerting selective pressures on E. maclovinus. A major effect of these variables is expected in the early development stages, where survival to recruitment is crucial. For instance, recruitment variability is explained by environmental stability and food availability (Houde, 2008). A low sea current velocity provides a stable environment, while primary productivity is an indirect proxy for food availability (Gove et al., 2016;Fox et al., 2018). Oxygen plays a role in energetic metabolisms and neuromuscular processes, which are important during embryonic development and environmental stress (Epelboin et al., 2016;Moreira et al., 2018). In addition, salinity modifies the egg and larvae distribution (Petereit et al., 2009) -Chacoff et al., 2015) and with food availability conditions (Vanella et al., 2017). Oxygen consumption as a proxy of metabolisms showed minor changes comparing gill and kidney tissues with liver tissue, suggesting that E. maclovinus can adapt to different salinities (Vargas-Chacoff et al., 2014). Although we did not find the association with temperature, F I G U R E 5 Directional relative migration network including all relative migration values among locations by (a) neutral loci, (b) adaptive merged loci, and (c) adaptive shared loci. Box plots (d) for one-way anova comparing differences among Nm distance and the three different datasets. *p < .05; ****p < .0001. REL, Reloncaví Estuary; MAN, Manao; HOR, Hornopirén; CHE, Chepu; CHA, Chaitén.
studies have shown that E. maclovinus could reorganize components for the intermediary metabolism to respond to climate change (Oyarzún et al., 2017). Additionally, it has been proposed a trade-off between growth and swimming activity to allocate energy, where juveniles allocate more energy to swimming in low temperatures than physiological functions like growth (Vanella et al., 2012).
All genetic groups identified with the adaptive datasets are distributed in an environmental cline and some locations are more influenced by marine conditions than others. Divergence in adaptive loci reinforces the idea that such markers may experience selection by environmental pressures (Hollenbeck et al., 2019). where samples were collected for this study. They found that factors such as salinity, water discharge by rivers, glacial melt, and precipitation cause differences in the habitat between the Reloncaví Estuary and the exposed area in Chiloe (Araneda et al., 2016).
Two main genes identified in our study, SNX19 and DYNC2H are candidates for salinity adaptation. The SNX19 participates in regulating vesicle trafficking, which can confer salt tolerance and the DYNC2H1, which has been associated with "renal water homeostasis," "vasopressin-regulated water reabsorption," and "urea transport" (Zhou et al., 2018).

| Management implications
Our results provide the first report of the fine-scale spatial pattern of structuration with gene flow and spatial selection by the found that males and females mature early in locations with intensified recreational or commercial fisheries (Hamilton et al., 2007).

| Future directions
We showed that using relatively dense genomic information allows us to refine the population structure pattern in E. maclovinus in a particular area in its distributional range, whereas previous studies indicated only weak genetic differentiation in a large geographical area (Ceballos et al., 2012(Ceballos et al., , 2016. Future research should include the whole distribution of the species; making it possible to assess more populations along and out of the heterogeneous landscape in Patagonia and identify whether environmental variables are associated with the current diversity and population structure. Additionally, other types of genomic variation should be examined such as copy number of variations or structural variants (e.g., Barth et al., 2019;Cayuela et al., 2021), which are also informative for identifying structuration patterns (Mérot et al., 2020). Knowing the polygenic architecture associated with adaptation to environmental variables -combining approaches of population genomics and quantitative genetics (Gagnaire & Gaggiotti, 2016)  Investigación" awarded by Universidad de los Lagos. We thank two anonymous reviewers for providing valuable comments of the manuscript and Susan Cleveland for first checking and improving the language or the English style.

FU N D I N G I N FO R M ATI O N
This research and the APC was funded by ANID/CONICYT FONDECYT Iniciación no. 11180897.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
All raw data analyzed for this manuscript have been uploaded to NCBI's SRA database (BioProject ID PRJNA883977). Genotypes are archived on DRYAD (https://doi.org/10.5061/dryad.mkkwh711v).