Population genomics and history of speciation reveal fishery management gaps in two related redfish species (Sebastes mentella and Sebastes fasciatus)

Abstract Understanding the processes shaping population structure and reproductive isolation of marine organisms can improve their management and conservation. Using genomic markers combined with estimation of individual ancestries, assignment tests, spatial ecology, and demographic modeling, we (i) characterized the contemporary population structure, (ii) assessed the influence of space, fishing depth, and sampling years on contemporary distribution, and (iii) reconstructed the speciation history of two cryptic redfish species, Sebastes mentella and S. fasciatus. We genotyped 860 individuals in the Northwest Atlantic Ocean using 24,603 filtered single nucleotide polymorphisms (SNPs). Our results confirmed the clear genetic distinctiveness of the two species and identified three ecotypes within S. mentella and five populations in S. fasciatus. Multivariate analyses highlighted the influence of spatial distribution and depth on the overall genomic variation, while demographic modeling revealed that secondary contact models best explained inter‐ and intragenomic divergence. These species, ecotypes, and populations can be considered as a rare and wide continuum of genomic divergence in the marine environment. This acquired knowledge pertaining to the evolutionary processes driving population divergence and reproductive isolation will help optimizing the assessment of demographic units and possibly to refine fishery management units.

the drivers that may explain genomic distance between individuals from different areas, by using multivariate methods to describe the association between explanatory (i.e., environmental variables) and response (i.e., genomic distances) variables (summarized in Riginos et al., 2016). Both population genomics and spatial ecology have their limitations, but combining these approaches allows these potential weaknesses to be better addressed and reinforces the interpretation of the observed marine connectivity scheme.
Furthermore, coupling population genomics and spatial ecology approaches with information about the demographic history of species and populations can deliver a broader portrait by delineating the moment and the geographic conditions (i.e., allopatry vs sympatry) under which the genetic clusters identified in the population genomics approach formed.
Genetic clusters can be formed following habitat fragmentation and climatic/environmental change, for instance as it occurred during the Quaternary period (Hewitt, 2000;Rougeux et al., 2017). This period characterized by change in important environmental driver has shaped the connectivity of populations and the dynamic distribution of marine species in the North Atlantic (Maggs et al., 2008). Demographic modeling is fundamental toward determining whether species diverge in the face of gene flow under the action of divergent ecological selection alone (i.e., ecological speciation, Rundle and Nosil, 2005) or if an initial isolation period precedes subsequent divergence (Endler, 1977;Roux et al., 2016). Methods, either using coalescent simulations Roux et al., 2016) or solving the site frequency spectrum through a diffusion approach (Gutenkunst et al., 2009), are suitable to test for these scenarios and also to model local effects in the genome, which can improve demographic inferences (Cruickshank and Hahn, 2014;Ewing and Jensen, 2016;Rougeux et al., 2017;Roux et al., 2014;Schrider et al., 2016; more details provided in the methods).
With more than 110 closely related species worldwide (Hyde and Vetter, 2007), the marine genus Sebastes has attracted the attention of evolutionary biologists that aim to study their speciation history (Johns and Avise, 1998). Redfish are ovoviviparous (i.e., internal fertilization), resulting in lecithotrophic larvae feeding exclusively on the yolk of the egg (Johns and Avise, 1998). In the North Atlantic Ocean, Sebastes mentella Travin 1951 and S. fasciatus Storer 1954 are two well-suited species to study speciation that have overlapping (i.e., sympatric) and nonoverlapping (i.e., allopatric) distributions. Both species co-occur in the Northwest Atlantic and inhabit at different depth intervals (Atkinson, 1987). The barrier to gene flow among Sebastes spp. in the North Atlantic Ocean is partially porous. In the Gulf of St. Lawrence (GSL) where S. mentella and S. fasciatus live in sympatry, S. mentella releases its larvae approximately three to four weeks earlier than S. fasciatus (Sévigny et al., 2000). Yet, evidence of introgression between S. mentella and S. fasciatus is present across their overlapping distribution but also between S. mentella and S. norvegicus or S. viviparous in the Northeast Arctic (Roques et al., 2002;Saha et al., 2017). Both S. mentella and S. fasciatus also show some intraspecific population structure.
For S. mentella, two ecotypes have been characterized across the North Atlantic Ocean, namely S. mentella shallow and S. mentella deep, which inhabit specific depths between 300 and 500 m and greater than 500 m, respectively (Cadrin et al., 2010;Saha et al., 2017;Valentin et al., 2014). Environmental gradients related to depth were suggested to play a role in the diversification of S. mentella ecotypes (Shum et al., 2014). S. fasciatus inhabits depths between 100 and 400 m and was previously characterized with weak population structure in the Northwest Atlantic using 13 microsatellite loci (Valentin et al., 2014). Along the Canadian coast, four populations were possibly present, namely (1)  however, the abundance of Sebastes spp. juveniles located in the northern GSL was reported to be 230 times and 4 times higher in 2017 than their respective average abundances for the period 1993-2012 (Senay et al., 2019). The stronger recruitment levels indicating new cohorts in the GSL may be a precursor to stock recovery and present hope for a potential future re-opening of the redfish fishery.
While population genetic structure in S. mentella and S. fasciatus has been intensively investigated using microsatellites, no study has yet used a genome-wide dataset to assess population genomic structure, spatial ecology, and demographic history to infer how contemporaneous and historical three-dimensional geographic drivers played a role in the formation of distinct genetic clusters, at the species, ecotype, and population levels. Here, our aim was to genotype S. mentella and S. fasciatus individuals sampled at different depths using genotype by sequencing (GBS) to (i) delineate the contemporary population genomic structure in Sebastes spp. located in the Northwest Atlantic Ocean, including the recent GSL cohorts observed, (ii) assess how space, fishing depth, and sampling years explain the contemporary distribution of species, ecotypes, and populations, and (iii) characterize the demographic history of speciation.

| Sampling and molecular techniques
From 2001 to 2015, we collected a total of 860 Sebastes spp. individuals, captured from 28 sampling sites in the Northwest Atlantic Ocean (Table 1 and Figure 1). Geographic coordinates of sampling sites are average latitude and longitude of sampling locations of all individuals considered in each site (see Table S1 for sampling location of each individual). Each site was sampled in only one year. For most sampling sites, adults were sampled except in the GSL and adjacent areas where juveniles of the new cohort were sampled ( Figure S1).
Individuals used for genotyping were selected based on their geometric morphometrics or their population classification using the 13 microsatellites developed by Valentin et al. (2014). Only individuals that were accurately identified at the species level and S. mentella ecotype level (i.e., shallow and deep only) were included in this study (Table 1 for detailed information). This sampling strategy consequently reduced the proportion of admixed individuals between the species and the ecotype in this study. This sampling bias is considered in the interpretation of the results. Genomic DNA was extracted using the DNA tissue kit protocol from Omega Bio-tek. DNA F I G U R E 1 Map of the 28 sampling sites of Sebastes mentella and S. fasciatus along the Northwest Atlantic. Delineation of the current management units specific to redfish species (i.e., Units 1 to 3) and NAFO divisions are presented. Each sampling site is represented by a black point, and adjacent colored circles represent genetic clusters detected at this site. A genetic cluster was indicated as detected when at least an individual shows at least 50% of ancestry to this genetic cluster. There are three ecotypes described in S. mentella: GSL (cyan), shallow (blue), and deep (dark blue). The unknown S. mentella genetic cluster is identified in light gray. There are five populations observed in S. fasciatus: I (purple), II (orange), III (red), IV (green), and IV (pink) quality and quantity were checked on 1% agarose gels. Libraries were prepared based on a GBS protocol modified from Mascher et al. (2013) using PstI andMspI (details in Perreault-Payette et al., 2017). Individuals were barcoded with unique sequences and pooled in multiplexes of 48 individuals per library. Libraries were then sequenced on Ion Torrent Proton P1v2 chips.

| Bioinformatics and genotyping
We used the process radtags module available in STACKS v.1.46 to demultiplex libraries (Catchen et al., 2013). Reads were truncated to 80 bp, and adapter sequences were removed with CUTADAPT (Martin, 2014). Then, reads were aligned to the reference genome of Sebastes aleutianus (id GCA_001910805.1_ASM191080v1) using BWA-MEM (Burrows-Wheeler Aligner-Maximal Exact Match) algorithm with default parameters (Li, 2013). Loci were required to have a minimum stack depth of three (m = 3) among reads with potentially variable sequences and identified using pstacks. The catalogue of consensus loci was created allowing up to three mismatches per sample. Genotypes were obtained after running populations module of STACKS v.1.46. SNPs were kept according to the following consecutive filtering steps and guidelines indicated in Benestan, Ferchaud, et al. (2016). First, SNPs genotyped in at least 60% of the individuals were retained. Potential paralogs were excluded by removing markers showing heterozygosity > 0.60 and −0.50 < F IS < 0.50 within sites, in more than 60% of the sites. Only SNPs showing a global minor allele frequency (MAF) > 0.01 or at least a MAF > 0.05 in one location were retained for the analysis. Only one SNP was kept per locus in order to avoid bias due to linkage disequilibrium Note: We used a code associated with the species identified (fas or men) and selected at the sampling location (01 to 28) to simplify graphic representations. We indicated the number of individuals genotyped per sampling site (Ngen) and the genetic group (Genetic group) inferred by Valentin et al. (2014) using 13 microsatellites.
TA B L E 1 Sampling sites analyzed in this study with year sampled, latitude, and longitude information as well as sampling depth (Table S2). This dataset was used for the analyses of population structure, the index of genetic differentiation F ST , the assignment tests, and the investigation of multivariate analyses. Another dataset of SNPs without the use of MAF filtering was also generated, with a total of 112,926 SNPs, to use with demographic models (more details in the next section). The resulting filtered VCF files were converted into the file formats necessary for the following analyses using PGDSpider v.2.0.5.0 (Lischer and Excoffier, 2012).

| Population structure and validation with assignment tests
Clustering among samples was inferred using the Bayesian clustering method implemented in the program ADMIXTURE v.1.23 (Alexander and Lange, 2011). Admixture analyses were run using 20,000 bootstraps, and the number of clusters was set from 1 to 29 (K) (n + 1 sampling sites). The support for different values of K was assessed according to the likelihood distribution (i.e., lowest cross-validation error) as well as a visual inspection of the coancestry values for each individual. Admixture analyses were also performed at the species level in order to investigate finer patterns of genomic structure. A PCA, implemented in the adegenet package (Jombart, 2008), was applied to the 834 individuals genotyped at both the 13 microsatellites (dataset from Valentin et al., 2014) and at the filtered 24,603 SNPs (see results), in order to compare the degree of resolution obtained using these two types of genetic markers. The first two principal components were visualized using ggplot2 package available in R (Wickham, 2011 polymorphic markers). The second run was achieved using the function assignment_ngs available in the package assigner (Gosselin et al., 2020) with the following parameters: sampling method = "ranked," thl = 0.3, iteration_subsample = 10, and assignment analysis = "gsi_sim." This function was designed to simultaneously test how many genetic markers were needed to accurately assign individuals at a species and population level while correcting for the problem of high-grading bias (Anderson, 2010).

| Drivers explaining genomic distances
We conducted distance-based redundancy analyses (db-RDA; Legendre and Andersson, 1999) to investigate the variables explaining redfish genomic variation among and within species. RDA is a direct extension of multiple linear regressions to model multivariate response data (Legendre and Gallagher, 2001 MEMs are derived from spectral graph theory and characterize a wide range of autocorrelation structures based on the survey design (i.e., distances between the n sampling locations; Dray et al., 2006).
Therefore, it is a spectral decomposition of the spatial relationships among the 59 GPS points of our entire dataset. This decomposition generates (n -1) eigenfunctions, which are orthogonal variables that can be used in statistical models as explanatory variables representing spatial structure at different scales among the studied locations.
Given the large spatial extent of the study area, distances among locations accounted for the earth curvature using the gcd.hf function in the codep package. To compute the MEMs, the truncation threshold was set to the length of the longest edge of the minimum spanning tree and only positive MEMs were retained using the dbmem function in the adespatial package. At three sampling sites (i.e., fas02, fas03, and fas12), depth was not recorded and therefore the median of depth distribution was used instead. Years were coded as dummy variables and treated as factors, therefore not assuming a linear relationship among years, but testing for annual differences in redfish genomic variation. To overcome multicollinearity among variables, we used the variance inflation factor with a threshold of 5, meaning that if a strong correlation was found between two variables, we retained only one of the two variables.
A global db-RDA model was run with all explanatory variables, and model probability and adjusted coefficient of determination (adjusted R 2 ) were calculated. The adjusted R 2 was quantified using the RsquareAdj function in vegan and accounts for the number of observations and number of degrees of freedom in the fitted model (Legendre et al., 2011;Peres-Neto et al., 2006). A selection of variables contributing to the explained variation was achieved by using both forward and backward selection, and a stopping criterion (or-diR2step function in the vegan package, Blanchet et al., 2008). This criterion limits overfitting by preventing selected variables to explain more variation than the full model developed with all explanatory variables. The same approach was performed within species, using a genomic distance matrix computed for S. mentella and S.
fasciatus separately, and subsets of explanatory variables associated with these individuals. Biplots were produced with a scaling of type 1 to illustrate distances among objects (i.e., individuals) and relationships with selected environmental variables. Species, ecotypes, and populations were identified by different colors in biplots to illustrate how genetic clusters were associated with selected variables. Maps of MEMs associated with redfish genomic variation were presented to visualize the most significant spatial scales.

| Demographic history of genomic divergence
We used the diffusion approximation framework implemented in dadi and modified by Tine et al. (2014) to compare the joint site frequency spectrum (JSFS) of our data to those simulated under different demographic models (Gutenkunst et al., 2009). These models incorporate the effects of selection at linked sites locally affecting population size (Ne) and the effect of differential introgression affecting the rate of migration (m) along the genome (Rougeux et al., 2017;Tine et al., 2014). During the divergence process, accumulation of barriers to gene flow tends to reduce the effective migration rate locally in the genome (Barton and Bengtsson, 1986), while selection at linked sites increases the effect of genetic drift through local reduction in effective population size (Charlesworth et al., 1993). Neglecting these effects can bias demographic inferences

| Two species with intraspecific geographic structure
Following the filtering steps shown in Table S2,  (ranging from 0.0004 to 0.1124)-than within S. fasciatus-on average 0.0257 (ranging from 0.0019 to 0.0715), which again was similar to ADMIXTURE finding a clearer population genomic structure in S. mentella than in S. fasciatus. Furthermore, the F ST heatmap showed that the sampling location men22 is genetically more similar to GSL than to shallow ecotype ( Figure S5), which is likely due to a mixed ecotype composition of GSL and shallow individuals at the men22 sampling site as indicated by the ADMIXTURE analysis (Figure 2c).

| Accuracy of assignment test at the species and ecotype level
We tested the validity of the clustering found by ADMIXTURE between and within species by performing individual assignment tests using all 860 individuals from both species. When considering the species and ecotype identity (shallow or deep) inferred from microsatellites and using the new 24,603 SNP markers genotyped, we obtained 100% assignment success at the species level. At the ecotype level for S. mentella, we assigned 100% of the individuals coming from GSL and deep ecotypes, while for shallow ecotype assignment, success varies from 68% to 97% (except for men16 where deep and shallow ecotypes were both sampled). The sampling locations men22 and men20 showed the highest percent of misassigned individuals-on average 32%-and these individuals were inferred to belong to S. mentella GSL, which corroborate the pattern seen in ADMIXTURE ( Figure 2c) and F ST analyses ( Figure S5).
Only a few percent of misassigned individuals sampled in shallow ecotype locations (i.e., men19, men20, men22, and men23) were assigned to deep ecotype (Figure 3a).  Figure 3b). The only population that matched a unique sampling site (fas12) was population IV (Figure 3b).
By selecting a subset of 258 "pure" individuals assigned with up to 90% of ancestry to their genetic cluster, using ADMIXTURE results, we then assess the number of SNP markers required to reach a high assignment success at the species and ecotype levels (i.e., > 90%). This analysis revealed that only two SNPs ( Figure S6) were needed to reach an assignment success of 100% at the species level while the use of a minimum of 2,000 highly divergent markers was required at the ecotype level for S. mentella (Figures 4 and S6).

| Inter-and intraspecific genomic distances explained by spatial distribution, fishing depth, and sampling years
As the geographic distributions of Sebastes species, ecotypes, and populations surveyed in this study seemed distinct but overlapping  (Table S3 and Figure 5B). The first db-RDA axis accounted for 59.3% of the   III, and V in orange, red, and pink, respectively, had a large distribution in the db-RDA space, and therefore, no clear relationships with explanatory variables could be identified (Figure 5c).
MEM1 and MEM4 were the spatial variables best explaining S. fasciatus population distribution. MEM1 outlined a large-scale east-west gradient. MEM4 also showed a large-scale gradient where northeast individuals were different from their southwest counterparts ( Figure S9).

| Evidence for secondary contact at both species and ecotype level
We tested four major different demographic models, and the model  Table 3 for uncertainties around parameter estimates). The secondary contact represented on average 10.5% of the total time of divergence for these SC2N models, indicating that the species remained isolated most of their divergence time. The effective population sizes (Table 3) and the ratio of effective population sizes under all SC models revealed that both species showed large and relatively similar Ne (i.e., ratio Ne Fas /Ne Men ~ 0.92). Estimates of migration rate were low (m 1 and m 2 < 8e-6; Table 3), suggesting that only a few migrants have been exchanged between species since secondary contact. We further found that approximately 50% of the genome displayed a reduced population size (Q = 0.50; Table 3) due to the action of linked selection. In these regions affected by linked selection, the Hrf indicated that the population size was 17% of the total population size.
For ecotype comparison, the median estimate of divergence time between S. mentella deep and shallow population was ~ 0.46 MYA under the SC2m models. The migration rate since secondary contact was higher than that observed between species and strongly asymmetric, with higher migration from the deep (m 12 = 0.0003) to the shallow population (m 21 = 9e-5 Table 3).
Overall, we found that ~28% (ranging from 26 to 38%) of the genome was undergoing reduced introgression since secondary contact (P parameter in Table 3). Finally, using haplotype data derived from the GBS loci, we found that the level of net divergence (Da) over all combined loci was an order of magnitude lower in TA B L E 3 Demographic parameter estimates obtained between the three major groups compared Note: AIC, Akaike's information criterion; log likelihood, maximum likelihood; theta: effective mutation rate of the reference population, which here corresponds to the ancestral population; Nref, size of the ancestral; Ne1 and Ne2, effective population size of the compared pair; m1 ← 2 and m2 ← 1, migration from population 2 to population 1 and migration from population 1 to population 2; me12 and me21, effective migration rate estimated in the most differentiated regions of the genome; Ts, time of split of the ancestral population in the two species; Tsc, duration of the secondary contact; P, proportion of the genome freely exchanged (1-P provides the proportion of the genome non-neutrally exchanged); Q, proportion of the genome with a neutral Ne (1-Q provides the proportion of the genome non-neutrally exchanged); hrf = Hill-Robertson factor representing the reduction of Ne in the region (1-Q) with reduced Ne.   (bottom panel) in the comparison between ecotypes (a) and between species (b and c). The suffix "2m" denotes heterogeneous migration and "2N" heterogeneous effective population size. AM = ancient migration. IM = isolation with migration, SC = secondary contacts, SI = strict isolation. All fittest models incorporate heterogeneous parameters, such as a heterogeneous rate of gene flow (2 m) and linked selection (2N). Each block corresponds to a population with effective population size of N1, N2, and θ. Horizontal dashed lines indicate the duration (T sc ) of the secondary contact period after the duration (T s ) of the allopatric period, while vertical dashed lines indicate change in effective size. Black arrows symbolize neutral gene flow, whereas red arrows symbolize heterogeneous introgression intraspecific genomic variations were best explained by large-scale spatial distribution and fishing depth. We also uncovered that postglacial secondary contact was the most likely scenario to explain the speciation history of species and S. mentella ecotypes.

| Highly structured Sebastes species in the northwest Atlantic
We observed a level of genetic differentiation spanning from weak to pronounced between the genetic clusters within both S. mentella and S. fasciatus, and all have distinct but overlapping distributions.
The largest estimate of genetic differentiation was between S. mentella and S. fasciatus species (F ST = 0.505; P-value < 0.001), supporting that these genetic clusters form independent demographic units.
Within S. mentella, we confirmed the presence of deep and shallow ecotypes and identified a new ecotype, GSL in the Northwest Atlantic with more genetic similarities to shallow than deep ecotypes. It is then possible that GSL ecotype may have diverged from the shallow ecotype. The shallow and deep ecotypes were previously identified and described as distinct genetic clusters inhabiting different depths along the continental slope in past studies across the North Atlantic (Cadrin et al., 2010;Saha et al., 2017;Valentin et al., 2014). Here, we suggest to use the term ecotype to describe these ecogeographically isolated groups. Ecotype was first defined and summarized as a genetic group specific to a particular habitat (Turesson, 1922). We considered that the three genetic groups of S. mentella described in this study correspond to that first definition and that results obtained from the db-RDAs have shown that the genetic distance between these ecotypes is explained by differences in habitats through space and depth.
Furthermore, the genetic composition of the S. mentella GSL ecotype shows consistent admixture with S. fasciatus (i.e., 18%) and it is the unique ecotype present in the GSL, a partially isolated area.
Within S. fasciatus, we identified five populations occurring from the Gulf of Maine to the Labrador Shelf with highly overlapping distribution. This contrasts with the more pronounced differentiation we observed between S. mentella ecotypes and consequently justifies the use of the term population (instead of ecotype) to describe these genetic clusters. Nevertheless, the five populations identified are mostly associated with distinct geographic areas along the Canadian coast (Figure 1), suggesting further investigation on the population structure of this species. The population I was the only genetic cluster that had up to 6% of S. mentella ancestry, which underlined that introgression may be the main cause of the distinctiveness observed in this particular population. Moreover, a higher genetic differentiation-in line with previous microsatellite observations (Valentin et al., 2014)-was detected for the semiclosed population located in Bonne Bay Fjord (fas12; population IV), for which we hypothesize that limited water exchange or limited physical connectivity may be responsible for this pattern, as previously proposed in the red algae Palmaria palmate (Li et al., 2015). Still, fas12 does not appear completely isolated and larvae may potentially migrate in the fjord since a mixture of two populations was observed at this location. Our results confirm that there is a unique population in Bonne Bay Fjord, which was in the past designated as "of special concern" because of its demographic isolation (COSEWIC, 2010); however, our results also suggest that the Bonne Bay Fjord could naturally be reseeded from other populations.
The  (Johns and Avise, 1998). Consequently, these genetically and ecologically distinct groups appear mostly closed (i.e., self-recruiting) and thus need to be individually protected for proper conservation (Hellberg et al., 2002). Although hybrids have been documented between these two species in the past, no F1 was documented in a nonbiased sampling design (Roques et al., 2001). Along with our result, this tends to indicate that introgression events may be ancient between the two species, and thus, introgression is not limiting the establishment of genome-wide barriers to gene flow between these species. Nevertheless, our current sampling design, which avoids potential hybrids by sampling only individuals clearly assignable to a priori groups, limits our interpretation of our results to properly answer to this question.
Although we identified pronounced intraspecific structure in S. mentella, some structure is still unexplained and more genotyping will be necessary to identify all the demographic units present in the Northwest Atlantic. Thus, an "unknown" genetic cluster was identified in S. mentella sampling zones (i.e., men21, men22, men23, and men24) that correspond to an area of contact with another Sebastes species, Sebastes norvegicus, which is present in the North Atlantic (Stransky, 2005). We also identified some admixture between the GSL and shallow ecotypes of the NAFO divisions from 3K to 2G.
These samples were from adults and suggest that past larval migra-

| Spatial distribution, depth, and other factors explaining inter-and intraspecific barriers to gene flow
Defining population structure in marine species remains challenging to investigate considering the absence of obvious barriers to gene flow (Palumbi, 1994). Relying on tools from spatial ecology along with population genomics, we found that spatial distribution explained most of the genetic differentiation in the redfish system while fishing depth and sampling years significantly explained from a moderate to a small proportion of the genomic distances between individuals of Sebastes spp. Spatial patterns were strong predictive variables for species, ecotypes, and populations, which is in concordance with the emerging results of recent seascape genomic studies (reviewed in Riginos et al., 2016) on the American lobster (Homarus americanus; Benestan, Quinn, et al., 2016), the sea cucumber (Parastichopus californicus; Xuereb et al., 2018), the red mullet (Mullus surmuletus; Dalongeville et al., 2018), and the Eastern oyster (Crassostrea virginica; Bernatchez et al., 2019). In particular, this north/south genomic structure may reflect an historical signature of a glacial refugium during the Pleistocene followed by population expansion throughout the northwest. Pleistocene glacial cycles associated with geological and climatic changes have indeed clearly played a role in shaping the contemporary distribution and abundance of marine organisms in the North Atlantic (Wares and Cunningham, 2001).

Depth may play a role in the isolation of these species and other
Sebastes species (Ingram, 2011;Shum et al., 2015). Despite the potential influence of bathymetry on marine organisms, no recent seascape genomics studies have already investigated and evidenced how this environmental feature can drive marine population genomics. Here, we uncovered that fishing depth was a significant explanatory variable for genomic distances between species, ecotypes, and populations of Sebastes spp. This outcome brings new evidence that habitat depth is a strong barrier to dispersal of Sebastes spp and contributes to the wide literature covering the topic specifically for S. mentella (Cadrin et al., 2010;Stefánsson et al., 2009) and other marine species such as the deep-sea coral (Desmophyllum dianthus; Miller et al., 2011), the common protobranch bivalve (Nucula attacellana; Jennings et al., 2013), and the small-spotted catshark (Scyliorhinus canicular; Kousteni et al., 2015). We also need to high-

| Physical separation followed by secondary contacts explains the speciation history for both species and S. mentella shallow and deep ecotypes
Our study is the first to extensively investigate the history of speciation for the Sebastes species inhabiting the North Atlantic using explicit demographic simulations. Here, we discovered that scenarios of secondary contact (SC) provided the best fit to the observed data for both species and ecotypes. Furthermore, modeling the process of speciation between S. mentella and S. fasciatus has never been done before, while for S. mentella deep vs. shallow comparison, Saha et al. (2017) previously identified the isolation with migration (IM) model as the best model to explain the divergence observed. However, the authors did not test for secondary contact and used a limited number of microsatellite markers, which are known to quickly converge to demographic equilibrium and therefore offer limited resolution to distinguish among models of divergence with gene flow (IM) and secondary contact (SC) . In contrast, SNP markers are more suited to perform demographic inferences given their lower mutation rate (Ellegren, 2000). Here, we used a much higher number of markers, thus reaching sufficient power to accurately distinguish these two models. Our parameter estimates suggested that the two species diverged first (~1 MYA) followed by a more recent split of the ecotype (~0.5 MYA). These divergence times correspond to a pronounced drop of sea level, suggesting that the redfish species may tend to become genetically isolated during periods of expansion/contraction associated with the Pleistocene glacial oscillations (2.6 million to 11,700 years ago), as previously suggested by Hyde and Vetter (2007 (Roux et al., 2016). Here, we found the Da estimate was ten times higher between species than ecotypes (Da = 0.0001 versus Da = 0.0015), which supports the idea that (i) ecotypes are more recently diverged than species, as well as the fact that (ii) species and ecotypes fall within the "grey-zone of speciation" (Roux et al., 2016). Overall, we highlighted that models integrating linked selection and gene flow provided the best fit to our data for both species and ecotypes.
These results supported a role for linked selection in shaping differentiation landscape, as is increasingly reported in the speciation literature (Burri et al., 2015;Cruickshank and Hahn, 2014;Rougemont and Bernatchez, 2018;Roux et al., 2016;Stankowski et al., 2019).

| Improving management of Sebastes spp. in the Northwest Atlantic Ocean
Despite the plethora of studies in ecology and genetics that have shed light on their limited movement and their large genetic differentiation (Roques et al., 2002;Roques et al., 2001;Valentin et al., 2014;Valentin et al., 2002), the two redfish species, S. mentella and S. fasciatus, are still often reported as "redfish"-due to their morphological similarities-without any differentiation in stock assessment, except for Unit 1 and 2 stocks where species-specific estimates are provided (DFO, 2020), and current Canadian management plans at some management units. Yet, reproductively isolated redfish aggregations can then be at risk of locally eliminating biological diversity and overfishing with the current large-scale management plan that does not reflect the biology of the species (Benestan, 2020). Since the collapse of the vast majority of the stocks in the Northwest Atlantic, the Canadian government has undertaken significant efforts to improve the recovery of redfish populations but more efforts may be required.
We identified unexpected asymmetric migration through the NAFO divisions 3K and 2G tend to then exchange more than 10% of individuals, which corresponds to a maximum threshold that can be used to infer demographic independence between two areas (Hastings and Harrison, 1994;Lowe and Allendorf, 2010;Waples and Gaggiotti, 2006). As the planning for possible re-opening of the redfish fishery in the GSL following its signs of recovery (Senay et al., 2019), the recovery of this stock may strongly influence adjacent areas to the GSL, such as the south of Labrador Shelf. We then urge for a conjoint discussion across Atlantic provinces and the south of Labrador Sea for the purpose of establishing a sustainable plan for this marine fish.
Genetic resources such as microsatellite markers were used for species identification of Sebastes species over the last two decades (Roques et al., 1999(Roques et al., , 2001(Roques et al., , 2002Valentin et al., 2014). Here, we found that only two SNPs were needed to accurately identify the two species. We also found that at least a subset of 2,000 SNPs could discriminate redfish S. mentella deep and shallow ecotype. This subset of 2,000 SNPs outperformed microsatellite markers for delineating ecotypes and improved assignment test success (i.e., ecotypes of S. mentella) as shown in a wide range of species with management or conservation concerns such as harbor porpoise (Phocoena phocoena; Lah et al., 2016), the brown trout (Salmo trutta; Lemopoulos et al., 2019), the Atlantic salmon (Salmo salar; Moore et al., 2014), and the great scallop (Pecten maximus; Vendrami et al., 2017). Therefore, a next step could be to develop a genotype assay (e.g., using Kaspar, Fluidigm, AmpliSeq, or Sequenom technologies)-based on the selection of these highly divergent SNP markers combined with an individual assignment test approach-in order to provide an efficient tool for traceability as previously validated for several marine resources in Europe (Nielsen et al., 2012). The DNA assay developed for Sebastes spp. would serve to distinguish the two species and three ecotypes within the fish catches in order to assess the true impacts of the redfish fishery in Canada.
Furthermore, a surprising high assignment success was achieved for