Genomic data support management of anadromous Arctic Char fisheries in Nunavik by highlighting neutral and putatively adaptive genetic variation

Abstract Distinguishing neutral and adaptive genetic variation is one of the main challenges in investigating processes shaping population structure in the wild, and landscape genomics can help identify signatures of adaptation to contrasting environments. Arctic Char (Salvelinus alpinus) is an anadromous salmonid and the most harvested fish species by Inuit people, including in Nunavik (Québec, Canada), one of the most recently deglaciated regions in the world. Unlike many other anadromous salmonids, Arctic Char occupy coastal habitats near their natal rivers during their short marine phase restricted to the summer ice‐free period. Our main objective was to document putatively neutral and adaptive genomic variation in anadromous Arctic Char populations from Nunavik and bordering regions to inform local fisheries management. We used genotyping by sequencing (GBS) to genotype 18,112 filtered single nucleotide polymorphisms (SNP) in 650 individuals from 23 sampling locations along >2000 km of coastline. Our results reveal a hierarchical genetic structure, whereby neighboring hydrographic systems harbor distinct populations grouped by major oceanographic basins: Hudson Bay, Hudson Strait, Ungava Bay, and Labrador Sea. We found genetic diversity and differentiation to be consistent both with the expected postglacial recolonization history and with patterns of isolation‐by‐distance reflecting contemporary gene flow. Results from three gene–environment association methods supported the hypothesis of local adaptation to both freshwater and marine environments (strongest associations with sea surface and air temperatures during summer and salinity). Our results support a fisheries management strategy at a regional scale, and other implications for hatchery projects and adaptation to climate change are discussed.


| INTRODUC TI ON
Intraspecific diversity is an important part of biodiversity, especially at higher latitudes where species richness is relatively low (Pamilo & Savolainen, 1999). Throughout the Palearctic, the periodic range contractions and expansions brought about by glacial cycles have shaped this diversity through bottlenecks and genetic drift (Hewitt, 2000) and particularly so in fishes (April et al., 2013;Bernatchez & Wilson, 1998). Local adaptation can also arise when species experience different environmental conditions over their geographic ranges (Kawecki & Ebert, 2004;Williams, 1966). Local adaptation has been studied extensively via reciprocal transplant and common-garden field experiments, but these approaches do not provide information on the molecular basis of adaptation (Tiffin & Ross-Ibarra, 2014). New genomic methods are now commonly used to advance our understanding of local adaptation (Grummer et al., 2019;Luikart et al., 2018). Such adaptive genomic variation, as well as contemporary population genetic structure, are of great interest for both conservation and management to ensure actions target biologically significant units Funk et al., 2012).
Salmonids are a diverse family of fishes with high economic and cultural importance. Anadromous salmonids are philopatric, that is, returning to their natal habitat for spawning (Quinn, 1993) and this behavior is known to reduce gene flow amongst populations, thus promoting genetic differentiation and local adaptation at fine spatial scales (reviewed in Fraser et al., 2011). The Arctic Char (Salvelinus alpinus) is a salmonid fish with a circumpolar distribution and is known for its great diversity in life-history characteristics (Klemetsen, 2010). Anadromous individuals spend 3-9 years in cold oligotrophic freshwater at birth (Johnson, 1980), then complete annual migrations between marine habitats for summer foraging and lakes for overwintering. Although straying (i.e., an upstream migration in a non-natal river system) can occur, several studies have shown that Arctic Char maintained philopatric behavior during reproductive years, limiting effective dispersal (Gyselman, 1994;Moore et al., 2013Moore et al., , 2017. Moore et al. (2013) also argued that gene flow could be sufficiently low to allow for local adaptation among populations of eastern Baffin Island in the Canadian Arctic, while Moore et al. (2017) provided some genomic evidence for local adaptation to natal rivers at a fine spatial scale.
Potentially strong demographic bottlenecks during glaciations have greatly reduced intraspecific genetic diversity and lead to divergent glacial lineages that survived in different refugia (Bernatchez & Wilson, 1998;Hewitt, 2000). Secondary contact between lineages has been shown to have commonly occurred among temperate freshwater fishes (e.g., Bradbury et al., 2015;Rougeux et al., 2019;Turgeon & Bernatchez, 2001). In North America, Arctic Char comprises four known mitochondrial DNA (mtDNA) lineages associated with distinct glacial refugia (Brunner et al., 2001;Moore et al., 2015). The Arctic lineage is predominant in the south-eastern Arctic (Brunner et al., 2001;Moore et al., 2015), but populations across northern Labrador are admixed with the Atlantic lineage (Salisbury et al., 2019).
Nunavik, situated in northern Québec (Canada), is one of the last regions in North America to have deglaciated following the last glacial maximum (LGM; Dalton et al., 2020;Dyke, 2004). It is bordered by the Hudson Bay, Hudson Strait, and Ungava Bay (Figure 1). These three marine regions are contrasted in their surface temperature, salinity, productivity, and tidal regimes (Prisenberg, 1984;Savard et al., 2014). In remote coastal communities of Nunavik, Arctic Char subsistence fisheries are key to food security (Laflamme, 2014) and rapid demographic growth has increased harvesting pressures on wild fish populations (Martin, 2011). The current study applies population genomic methods to investigate the relative role of neutral vs.
adaptive processes in shaping contemporary population structure of anadromous Arctic Char over 2000 km of coastline in Nunavik. This will provide a genomic context of relevance for sustainable fisheries management.
Since we expected the contemporary genetic structure of Arctic Char populations in Nunavik to be heavily influenced by their recent postglacial recolonization history, we first sequenced the mitochondrial control region to (i) test for the presence of the Atlantic lineage in Nunavik, placing the present study in the context of previous phylogeographic work in the region that relied heavily on mtDNA. We then used genotyping by sequencing (GBS) to determine whether Arctic Char populations from different regions of Nunavik (ii) are genetically differentiated, (iii) differ in their genetic diversity, and (iv) follow patterns of isolation-bydistance. Additionally, as migration in Arctic Char has been found to be predominantly coastal (Moore, 1975;Moore et al., 2013Moore et al., , 2015Spares et al., 2012), we tested whether (v) the Hudson Strait, a 120-km-wide open water area, acts as a barrier to gene flow between populations on opposite shores. Finally, since we expected broad-scale variation in environments to be the source of divergent selective pressures, we attempt to (vi) detect evidence of local adaptation to both freshwater and marine environments around Nunavik and bordering regions.

| Sampling and DNA extraction
Arctic Chars were sampled in 23 water bodies across Nunavik, southern Baffin Island, and Labrador ( Figure 1). In Nunavik and Baffin Island, adult fish were harvested during their upstream migration using gillnets or counting weirs. In most localities throughout Nunavik, sampling locations were selected in concert with local and regional Inuit wildlife managers, and sampling was done with the assistance of local Inuit guides to prioritize fish populations with an importance for traditional fishing. In two locations (Deception Bay and Hopes Advance Bay), samples were taken in the estuary as well as in two tributary rivers. Juvenile fish were captured by electrofishing in two rivers near Nain, Labrador. Adipose fin clips were collected from each fish and preserved in ethanol 95%. DNA was extracted from fin clips using a modified version of Aljanabi and Martinez (1997). Agarose gel electrophoresis was used to assess DNA quality, and the quantity and quality of DNA were evaluated by NanoDrop spectrophotometer (Thermo Scientific).

| Mitochondrial DNA analysis
To assess whether the Arctic and Atlantic glacial lineages hybridized in Nunavik, we sequenced mtDNA for a subset of individuals from sampling sites in the eastern part of our study area (n = 10-34 per site).
The entire mitochondrial control region was amplified with primers Tpro2 (Brunner et al., 2001) and SalpcrR (Power et al., 2009)

| GBS library preparation, sequencing, and single nucleotide polymorphisms (SNP) calling
Ten (10) µl of DNA was normalized to a concentration of 10 ng/µl using Quant-iT Picogreen dsDNA Assay Kit (Invitrogen) for precision quantification. GBS libraries were prepared with a modified version of the Abed et al. (2019) two-enzyme GBS protocol, using PstI and MspI restriction enzymes. Samples were randomly assigned to libraries to limit batch effects. Sequencing was performed at the Plateforme d'Analyses Génomiques (PAG) at IBIS, Université Laval (http://www.ibis.ulaval.ca/) on Ion Torrent sequencers with p1v3 chips and a median target of 80 million single-end reads per chip.
Each library was sequenced on three separate chips for each batch of 96 individuals, and the volume of DNA from each sample was adjusted after the first chip to reduce the unbalanced representation of individuals in sequences.
We processed the data and filtered the SNP dataset using a RADseq workflow (https://github.com/enorm andea u/stacks_workflow) built around STACKS 2.5 (Rochette et al., 2019). Briefly, the sequences were trimmed at 80 bp and aligned on a Salvelinus sp. reference genome (ASM291031v2; NCBI RefSeq: GCF_002910315.2; Christensen et al., 2021) using bwa mem (-k 19 -c 500 -O 0,0 -E 2,2 -T 0) in BWA-0.717 (r1188; ) and samtools view (-Sb -q 1 -F 4 -F 256 -F 2048) in samtools 1.8 . SNPs were called on polymorphic genotypes with at least 4X coverage F I G U R E 1 Sampling locations in Nunavik (Québec, Canada) and bordering regions. Red extents numbered 1-6 are magnified to show neighboring sampling locations. Samples were either used for genotyping by sequencing (GBS) only (round) or for GBS and mitochondrial DNA sequencing (square). Catchment area for each site is displayed in green per sample (-m), present in at least 60% of samples of each sampling sites (-r) and with the minor allele present in a minimum of three samples. This minor allele sample (MAS) filter is akin to minor allele frequency (MAF), but unlike MAF it is not artificially boosted by frequent RADseq genotyping errors (Linck & Battey, 2019). Samples with more than 20% of missing data or with heterozygosity (F IS ) under −0.2 were removed. For SNP calling and quality control, we combined our samples with 119 other samples from three sampling sites that were not included in our final dataset, as they either had too few samples, had no access to sea, or targeted a population reintroduced from a hatchery brood.
Salmonid fishes have a common ancestor that experienced a whole-genome duplication approximately 60 MYA (Crête-Lafrenière et al., 2012), and many genetic markers identified in our analyses are expected to be situated on paralogous loci of similar sequences.
While these loci may be important for adaptation (Kondrashov, 2012), they were removed due to the fact that these markers do not behave like biallelic SNPs and because genotyping is difficult without very high coverage (>100 reads; Dufresne et al., 2014). SNPs on duplicated loci were categorized and filtered using an adapted HDplot procedure (McKinney et al., 2017), which identifies paralogs by visually comparing the allelic ratio, the proportion of heterozygotes and homozygotes of the rare allele, and the F IS value for each SNP.
We also filtered the markers to avoid physically linked SNPs while keeping a maximum of information: For each pair of SNPs within 50,000 bp, we assessed linkage considering samples without missing data where at least one of the two genotypes contains the rare variant. If the two markers had identical genotypes in more than 50% of these samples, the pair was considered linked and only the first SNP was kept. To mitigate the risk of our inference regarding population structure being affected by an uncharacterized sex-ratio bias (Benestan et al., 2017), we identified sex-related markers by performing a redundancy analysis (RDA) using individual genotypes from samples for which information on sex was available in Hudson Bay (n = 69) and Ungava Bay (n = 85). We classified SNPs as showing a statistically supported association with sex when they loaded with more than 2.5 standard deviation from the mean (p = 0.012) and removed those markers from the dataset. Finally, we measured pairwise relatedness (Yang et al., 2010) in vcftools v0.1.13 (Danecek et al., 2011) and eliminated five related individuals in separate sampling sites, as their high relatedness values are thought to be the result of technical artifacts.

| Identification of putative neutral markers
Markers potentially under selection were identified using two methods: (1) pcadapt (Luu et al., 2017) and (2) BayeScan (Foll & Gaggiotti, 2008). SNPs identified as outliers by at least one method were removed to produce a neutral dataset. The R package pcadapt was used to identify outlier SNPs in relation to population structure according to a principal component analysis (PCA). The first 13 PCs were used, based on visual evaluation of PCA scores and scree plots, and SNPs with minor allele frequencies under 0.05 were excluded from the analysis. We then used BayeScan with prior odds of 1000 and other parameters set to default. For both tests, SNPs with a false detection rate (q-value) under 0.05 were considered putatively under selection.

| Basic statistics and population structure
Population genetic statistics were computed from the neutral dataset. Observed and expected heterozygosity was calculated by population using GenoDive v3.0 (Meirmans & Van Tienderen, 2004). We used vcftools v0.1.13 to measure the proportion of heterozygous SNPs for each individual, and we calculated the number of polymorphic SNPs in each population. Effective population sizes (N e ) and 95% C. I. were estimated with Neestimator v2.01 (Do et al., 2014) using the linkage disequilibrium method on markers with minor allele frequencies over 0.05.
We used a principal coordinate analysis (PCoA) in the R package dartR (Gruber et al., 2018) to document population structure using again the neutral data set. We also estimated ancestry with the maximum-likelihood approach implemented in ADMIXTURE (Alexander et al., 2009) with the number of genetic clusters (K) ranging from 1 to 20. We considered the value of K yielding the lowest cross-validation error to be the most likely number of genetic groups.
Based on this K value, we identified contiguous sampling sites where most individuals shared a common cluster membership and repeated the ADMIXTURE analysis within those sites with K ranging from 1 to 6. Considering that adult Arctic Char is expected to stray to nearby rivers, we regrouped sampling sites in connecting rivers and estuaries for the following analyses on genetic diversity and isolation-bydistance, with the exception of cases where ADMIXTURE showed strong population structure.
Spatial patterns of genetic diversity were compared between a priori regions, using the individual proportion of heterozygous SNPs.
A nested ANOVA was performed with the regions as groups and the sampling sites as subgroups. To account for our unbalanced sampling design, we did a comparison of estimated marginal means (or leastsquare means) on a linear mixed-effect model with a random effect of the sampling site, in the R packages glmmTMB (Magnusson et al., 2017) and emmeans (Lenth, 2021).

| Landscape genomics
Pairwise population F ST were calculated with the R package StAMPP (Pembleton et al., 2013), and 1000 bootstraps were performed to estimate their significance value as the proportion of bootstrapped F ST values under zero. F ST calculations were performed on both the neutral and putatively adaptive datasets. The geographic marine distance separating sampling sites was measured between the coordinates at the mouth of sampled rivers by a least-cost path in the R package marmap (Pante & Simon-Bouhet, 2013) using NOAA bathymetric data at a 4-minute resolution to discriminate land and sea.
We tested for the presence of isolation-by-distance (IBD) with a linear mixed-effect model. Linearized Slatkin, 1995), based on the neutral dataset, was used as the dependent variable and marine distance as a fixed effect. Patterns of IBD between pairs of sampling sites on the same coast versus pairs on either side of Hudson Strait were compared by the inclusion of a categorical fixed effect. To characterize the impact of Hudson Strait as a barrier to gene flow at different spatial scales, analyses were performed on all nonestuarine sites, then only on nonestuarine sites within 250 km of the Hudson Strait. To account for the nonindependence of pairwise distances, the model was run with a maximum-likelihood population effect parameterization (MLPE) (Clarke et al., 2002) with and without restricted maximum likelihood (REML), using the MLPE.lmm function in the R package ResistanceGA (Peterman, 2018). Models without REML were compared with conditional Akaike information criterion (cAIC; Vaida & Blanchard, 2005) in the R package cAIC4

| Gene-environment association
We used the ArcGIS software v10.4 (ESRI, 2011) to extract environmental data from BIO-Oracle v2.0 (Assis et al., 2018), Marspec (Sbrocco & Barber, 2013), and WorldClim v2.0 (Fick & Hijmans, 2017) (Table 1, see Table S5 for values at sites). Tide data were obtained from FES2014, produced by Noveltis, Legos, and CLS and distributed by Aviso+, with support from Cnes (https://www. aviso.altim etry.fr/). Marine variables represent sea-surface values for factors of potential biological importance for Arctic Char and were aggregated in a 20 km radius around each river mouth and within 5 km from the coast, as to best represent the local coastal environment based on existing knowledge of Arctic Char marine habitat use from other geographical regions (Moore et al., 2016;Spares et al., 2015). Freshwater variables comprise the area of the watershed upstream of the sampling site, as well as air temperature and precipitation statistics on these areas. Air temperature is commonly used as a proxy for freshwater temperature in remote areas, as supported by studies linking the growth rate of Lake Trout to air temperature (Black et al., 2013;Torvinen, 2017). To deal with multicolinearity, a PCA was performed on the scaled environmental factors, and the PCs with eigenvalues over 1 were used as explanatory variables in gene-environment association (GEA) analyses.

TA B L E 1 Value range and source of environmental factors considered in gene-environment associations. PCA axis mostly associated with variables is indicated with sign (positive or negative) of the correlation in parentheses
Candidate SNPs associated with environmental variables were identified using the complete genomic dataset with a combination of two redundancy analyses (RDA) and latent factor mixed models (LFMM). Genes within 10,000 base pairs were recorded using bedtools (Quinlan, 2014) (Jombart, 2008). We then used the R package adespatial (Dray et al., 2020) to compute the dbMEMs from the Euclidian distances. This approach was chosen over mere longitude and latitude information to avoid the underestimation of marine migration distances, that is, by preventing considering for instance impossible movements through landmass. Eigenvectors reflecting positive spatial autocorrelation were used as covariables in a partial RDA. As spatial structure could potentially hide evolutionary significant environmental gradients, we performed a second RDA excluding spatial eigenvectors correlated to environmental factors, as suggested by Forester et al. (2018).
Spatial factors were first submitted to a backward model selection procedure with the function ordistep in vegan, and only statistically supported covariables at α = 0.05 were kept. We then repeated this model selection method for environmental factors for both RDAs. To identify SNPs associated with environmental components, Z-scores were obtained for the distribution of individual SNP loadings on RDA axes explaining a significant portion of genetic variation (p < 0.05), and SNPs were defined as outliers if their maximal absolute Z-score was over 3.5 (p < 0.0002).

| LFMM
The LFMM (Frichot et al., 2013) identifies loci-environment associations using a Bayesian mixed model with environmental variables included as fixed effects. Latent factors are derived from a PCA and used as random effects to control for population structure. Missing data in the genetic dataset were imputed based on the most frequent genotype in the sampling site and we build the model using the lfmm_ridge function of the R package lfmm (Caye et al., 2019). The number of latent factors (K value) and regularization parameters (lambda) were respectively set to 9 and 1e −5 , to minimize predictor error estimated by a cross-validation method, as advised in the lfmm manual. p-Values were calibrated using the genomic control method, and false discovery rate (q-value) was calculated following the Benjamini-Hochberg procedure in the qvalue R package (Storey et al., 2020). Associations between a SNP and environmental factor with q-value < 0.01 were considered statistically supported.

| RE SULTS
Based on mtDNA, we only detected secondary contact between glacial lineages in Labrador: Samples from ANA and KAM cor-  (Table S1).
Using GBS, we sequenced a total of 745 samples, of which 95 did not pass our filtering criteria. Retained individuals (n = 650) had 2.217 million reads on average. A total of 30,773 SNPs were called and passed basic filtering in STACK2. Among these, 7244 SNPs were categorized as duplicates and removed ( Figure S1) and 5009 SNPs were pruned during linkage disequilibrium assessment. We then identified and removed 408 SNPs linked to sex, 139 of which mapped on the sex chromosome (15.6% of SNPs on this chromosome) and 269 were elsewhere in the genome (1.5% of SNPs outside the sex chromosome). The final dataset used in subsequent analyses comprised 18,112 SNPs (see Table S2 for exact criteria and number of filtered SNPs), with a global 6.87% missing genotypes.

| Population structure
Population statistics are presented in Table 2.  (Table S3), individual cluster membership for K > 1 was then at most loosely linked to sampling sites ( Figure 3b). However, Labrador sites, where samples were collected as juveniles, are a notable exception.
Individual mean proportion of heterozygous markers varied regionally ( Figure 4). Notably, southern Hudson Bay displayed the TA B L E 2 Summary of Arctic Char sampling and associated descriptive statistics. Samples were collected on adults in rivers and lakes, with exceptions ( † in estuaries, ‡ on juveniles)   (Table 3). However, similar analyses focusing on sampling sites within 250 km of the Hudson Strait supported that the slope and intercept were different for pairs of populations on either side of the Hudson Strait and pairs on the same coast, suggesting that populations separated by the Hudson Strait are more differentiated (Figure 5c).

| Gene-Environment Association
Environmental factors were summarized with three PC axes (Table 1). The first axis (PC1; 38.0% of variation) was associated with primary productivity, dissolved oxygen, precipitation, and winter air temperature. The second axis (PC2; 26.5%) was associated with summer temperature, both in the marine and in freshwater habitats, as well as to turbidity and salinity. Finally, the third axis (PC3; 17.2% of variation) related mostly to tidal amplitude and partly to the size of the catchment area ( Figure S4, S5).
MEM1 was highly correlated with PC3 (r = 0.71) and MEM2 was linked to PC1 (r = 0.74). We conducted a redundancy analysis (hereafter referred to as "RDA") which excluded the spatial covariables Respectively 98, 113, and 84 outliers were found on each of these three axes ( Figure 6). Those outliers were generally most correlated to PC1 (n = 185) or PC2 (n = 90). LFMM identified a total of 173 outlier SNPs (q < 0.01, Figure S10), 13 of which being associated PC1, 102 associated with PC2, and 58 to PC3. A total of six SNPs were significantly associated with more than one environmental component, five of which were associated with both PC2 and PC3.
Of the 585 unique outliers, 57 (9.7%) were also among the loci considered as potentially under selection by either genome scan methods (i.e., pcadapt and Bayescan). This relatively low proportion is not unexpected, given that genome scans are mostly sensitive to largeeffect loci, while RDAs have the power to detect polygenic selection (Forester et al., 2018). A total of 93 (15.9%) outliers were identified as associated with an environmental factor by at least two GEA methods ( Figure 6e, see Figure S11 for spatial distribution of allele frequencies and Table S6 for the number of outliers per linkage group). Among 15 of those 93 candidate SNPs, the LFFM identified PC2 as the environmental factor associated with the candidate SNPs while the factor most correlated to genotypes in the RDA/pRDA was PC1. In those cases, we reported PC2 as it was also statistically correlated to genotypes (r > 0.45, p < 0.05). Thus, 15 candidate SNPs were associated with F I G U R E 2 (a) Population structure assessed by a principal coordinate analysis (PCoA). Individual scores on PCoA axes 1 and 2 are presented as points and colored by a priori geographical region. An ellipse representing a 95% confidence interval was drawn around each sampling site. The percentage of genetic variance explained by each axis is in parentheses PC1 (precipitation, winter air temperature, and productivity), 61 were associated with PC2 (air and sea-surface summer temperature, salinity, and turbidity), and 17 were associated with PC3 (tides). We identified a total of 75 named genes within 10,000 bp of those candidate SNPs, covering a wide range of biological functions, including muscle contraction, development, and circadian rhythm (Table S7).

| D ISCUSS I ON
Given their homing behavior, anadromous salmonids are considered prime models for understanding local adaptation (Fraser et al., 2011;Hendry & Stearns, 2004) and landscape genomics studies have provided key information for the conservation and management of economically important species (Waples et al., 2020). Here, we documented patterns of neutral and adaptive genetic variation in anadromous Arctic Char populations in the Nunavik region using genomic data. By combining fine-and broad-scale sampling, our results revealed a hierarchical genetic structure whereby neighboring hydrographic systems generally harbored distinct populations that could be grouped within oceanographic basins, while differentiation was also influenced by isolation-by-distance. To assess the potential for local adaptation to the contrasted environments in our vast study area, we combined three GEA methods. By doing so, we found genomic signatures that were consistent with local adaptation to either or both freshwater and marine habitats.

| Neutral structure in a postglacial context
By 8000 YBP, most of Hudson Strait and the coast of Labrador was free of ice, while the coasts of Hudson Bay and Ungava Bay deglaciated around 6500-7200 YBP (Dalton et al., 2020). The decreasing genetic diversity we observed along the coast from Hudson Strait to southern Hudson Bay, as well as low F ST values between populations within the Hudson Bay region, is therefore consistent with expectations following a recent range expansion (Eckert et al., 2008;Goodsman et al., 2014). Similar patterns were found in the Coho Salmon (Oncorhynchus kitshush) (Rougemont et al., 2020) Dalton et al. (2020), is represented by shades of gray for 6000-12,000 years before present (YBP) genetic diversity was lower in more recently deglaciated regions (Cauwelier et al., 2018). Following the initial recolonization of rivers in Nunavik by Arctic Char, processes such as isostatic rebound and the formation of proglacial lakes considerably modified the drainage basins (Dubé-Loubert et al., 2018;Jansson, 2003), which could have contributed to gene flow between regions (e.g., Ruzzante et al., 2020).
In Northern Labrador, Salisbury et al. (2019) found extensive admixture of the Arctic mitochondrial lineage, which crossed the Canadian Arctic Archipelago from a glacial refugium in the western Arctic (Brunner et al., 2001;Moore et al., 2015), and the Atlantic lineage, which we expect to have crossed the Atlantic Ocean from the Palearctic during the last deglaciation (Brunner et al., 2001;Wilson et al., 1996). Our results were consistent with these previous studies as we did not observe any Atlantic lineage haplotypes in populations outside Labrador. However, similar to Moore et al.
(2015), we found evidence for introgression from the Atlantic lineage in our nuclear data from southern Baffin Island, and even stronger evidence in Eastern Ungava Bay, where our ADMIXTURE results showed shared ancestry with Labrador samples. Such mito-nuclear discordance is not uncommon in cases of adaptive introgression and could also be linked to sex-biased dispersion (Toews & Brelsford, 2012) as it is suggested that Arctic Char males are more mobile than females (Dempson & Kristofferson, 1987;Moore et al., 2016). Nevertheless, admixture in the eastern part of our study area could explain the higher genetic diversity, as expected during secondary contact of marine fishes (Bay & Caley, 2011;Grant & Bowen, 1998). Pairs of populations on opposite shores of the Hudson Strait did not deviate from IBD, which, contrary to our expectation,  Arctic Char is known for its higher straying rates than most other salmonids, but many studies argue that this dispersal does not necessarily lead to gene flow, as individuals are more prone to straying when overwintering than when spawning (Dempson & Kristofferson, 1987;Moore et al., 2013;Moore et al., 2017;Sévigny et al., 2020).

F I G U R E 4 Individual proportion of heterozygous SNP markers
In this study, spawning and nonspawning individuals were not distinguished upon sampling, which might have resulted in underestimating the real genetic differentiation among rivers. Sites sharing an estuary generally did not exhibit population structure, indeed suggesting high levels of straying between nearby rivers, but our ADMIXTURE results indicate that straying is very rare at the inter-regional scale this study mostly focusses on. Surprisingly, however, the Leaf River (FEU) and the Bérard River (BER) showed strong evidence of genetic differentiation despite sharing an estuary. Interestingly, those two rivers are found on either side of the limit of two major geological provinces (Thériaut & Beauséjour, 2012). As salmonids rely partly on their olfaction to recognize their natal rivers (Keefer & Caudill, 2014), the distinct geologies of the Leaf and Bérard systems could perhaps improve their homing ability and thus limit straying in these systems. Finally, the two sites in Labrador also displayed clearer signs of differentiation than other pairs of neighboring rivers, perhaps because the sampled juveniles did not have the opportunity to stray yet. On the other hand, gene flow between Labrador populations might also be lower than in other regions of the Canadian Arctic, as Li et al. (2020) found population structure in the Lower Northwest Passage, Nunavut, to be similar to what we measured in Nunavik.

| Genomic evidence consistent with local adaptation
We explored how genetic variation in Arctic Char was linked to a range of climatic and physical environmental predictors while considering both freshwater and marine habitats, something that has rarely been done in salmonids (but see Bekkevold et al., 2020 (Bourret et al., 2013;Dionne et al., 2008;Perrier et al., 2017;Sylvester et al., 2018). In contrast, selective pressures in marine habitats have been argued to be weaker since observed mortality rates are lower at sea (Garcia de Leaniz et al., 2007;Quinn, 2005).
However, SST near the mouth of spawning river was found to be correlated with timing of migration in species of Pacific salmon (Kovach et al., 2015) and to phenology-related genes in Arctic Char (Madsen et al., 2019 Recent studies of local adaptation in salmonids have focused on tributary-specific variation in freshwater conditions within a single catchment (e.g., the Columbia River; Hand et al., 2016;Hecht et al., 2015;Micheletti et al., 2018). However, the freshwater environmental factors used in our study are catchment-based given our sampling strategy, which prevents us from knowing the precise spawning site or overwintering lake. In marine systems, genomic evidence for local adaptation and isolation-by-environment has been found both at local scales in heterogeneous habitats (e.g., Lehnert et al., 2019;Miller et al., 2019) and over large geographic distances (e.g., Clucas et al., 2019). Arctic Char is expected to use preferred habitats based on temperature, salinity, and prey availability Spares et al., 2012Spares et al., , 2015. As we averaged near-shore marine conditions around river mouths, this study is limited to broad-scale environmental heterogeneity. This is in line with the suggestion by Fraser et al. (2011) that adaptation of anadromous salmonids to the marine environment should occur at a larger spatial scale than in fresh water.
Regardless of the geographic scale being studied, an ideal sampling design for detecting local adaptation should maximize the environmental variation while minimizing its collinearity with neutral genetic patterns. For example, Lotterhos and Whitlock (2015) suggested sampling pairs of populations with similar ancestry and contrasting environments. Such a design is suitable when studying variables that may change drastically over short distances, such as catchment area and upstream migration distances. However, climatic and physicochemical conditions experienced by geographically close populations of Arctic Char are more likely to be similar than in distant ones. Similar considerations were discussed in Nadeau et al. (2016), where patterns of IBD, isolation-by-environment, and isolation-bycolonization were difficult to disentangle for two pine species in a recently recolonized range. In our study, some genotypes for GEA candidates ( Figure S12) have varying allele frequencies that follow spatial patterns reminiscent of the neutral structure we described earlier, especially in PC1, where environmental variation followed a longitudinal gradient. In some of those cases, it is possible that a GEA was detected even though neutral processes could better explain the distribution of the observed allele frequencies. For example, if colonization of Hudson Bay did occur by rapid demographic expansion as discussed earlier, allele-surfing events, where a mutation can reach high frequency by chance alone (Edmonds et al., 2004), could explain differential allele frequencies (Rougemont et al., 2020).
Alternatively, the introgression of the Atlantic glacial lineage in the eastern part of the study area could also have led to divergence due to drift during the LGM to be falsely identified as linked to environmental variation. On that note, the pRDA and LFMM methods all accounted in some way for neutral or spatial structure. However, it is noteworthy that even the RDA, while excluding spatial correction, only identified a few outliers whose variation was driven by the Labrador sampling sites, that is, from the Atlantic glacial lineage.
As we found a significant part of genetic variation to be explained jointly by environmental variation and spatial patterns, we advocate that including GEA methods that do not account for spatial structure offer a way to acknowledge that local adaptation can also contribute to genetic structure across contrasting environments, for instance by selecting against maladapted migrants (Dionne et al., 2008;Wang & Bardburd, 2014). In that sense, local adaptation in the system studied here could reinforce the neutral structure discussed earlier.

| Implications for conservation and management
Genomic data have the potential to improve the definition of management and conservation units by accounting for neutral genetic structure as well as local adaptation Funk et al., 2012). Arctic Char populations in Nunavik support important small-scale subsistence fisheries, and stocks are managed on a riverby-river basis on the premise that each river contains a single and distinct population (Johnson, 1980). However, the weak genetic differentiation in rivers sharing an estuary and our inability to identify substructure at this level suggest frequent straying of adults, in line with the evidence for the prevalence of mixed stocks of Arctic Char in adjacent rivers (Boguski et al., 2016;Moore et al., 2014Moore et al., , 2017. If genomic tools were to be used for stock assignment of adult fish (e.g., Meek et al., 2016;Moore et al., 2017), further local sampling including juveniles and/or reproducing adults would be useful to better understand fine-scale population structure. The evidence we found for local adaptation to environments is in concordance with the neutral genetic structure at a broader scale, as major oceanographic basins around Nunavik are contrasted both in their environments and in the ancestry of their Arctic Char populations. We therefore suggest that these regional differences between Hudson Bay, Hudson Strait, and Ungava Bay, with distinction of the western and eastern coast for the latter, should form the basis of management actions at a regional level.
Also, there is growing interest in Arctic Char hatchery projects in Nunavik, both for supplementation and for reintroduction of Arctic Char in traditional fishing locations (George, 2007;Rogers, 2015).
The genetic information gathered here could be useful for this purpose, and the adaptative variation explored in this study highlights the need for careful selection of source populations for brood stocks, as maladapted domesticated individuals can waste efforts and resources, in addition of likely being detrimental to wild populations (Fraser et al., 2010;Tymchuk et al., 2006).
As the Arctic warms at a greater pace than any other regions on earth (Cohen et al., 2014), there might be concerns about the response of Arctic Char populations to their changing environment. Traits that are currently optimally adaptive in the present environment could eventually become maladaptive. Species will thus likely need to shift their distributions poleward and/or will need to adapt to persist in their current distribution if there is presence of appropriate genetic diversity/phenotypic plasticity. Layton et al. (2021) highlighted the genomic vulnerability to the future climate in Arctic Char from Labrador, especially at lower latitudes where higher temperatures are expected to be observed.
A temporal study recently showed that Arctic Char populations in Greenland have exhibited stable genetic structure over the last 60 years in face of rapid climate change, and argued that gene flow, although low, could allow for a modest level of evolutionary rescue in the short term (Christensen et al., 2018). Our study shows potential for local adaptation of Arctic Char populations to both their marine and freshwater habitats. As changes in climate might operate at a different pace, scale, and stability in marine and terrestrial ecosystems (Burrows et al., 2011), there is a need for continued research about the interaction of selective pressures over the lifespan of anadromous organisms.

ACK N OWLED G M ENTS
This study was funded as part of the BriGHT project (Bridging the Global change, Inuit Health, and the Transforming Ocean) by Sentinel North through the Canada First Research Excellence