Original Article Genetic structuring in Atlantic haddock contrasts with current management regimes

The advent of novel genetic methods has made it possible to investigate population structure and connectivity in mobile marine ﬁsh species: knowledge of which is essential to ensure a sustainable ﬁshery. Haddock ( Melanogrammus aegleﬁnus ) is a highly exploited marine teleost distributed along the coast and continental shelf on both sides of the North Atlantic Ocean. However, little is known about its population structure. Here, we present the ﬁrst study using single-nucleotide polymorphism (SNP) markers to assess the genetic population structure of haddock at multiple geographic scales, from the trans-Atlantic to the local (fjord) level. Genotyping 138 SNP loci in 1329 individuals from 19 locations across the North Atlantic revealed three main genetic clusters, consisting of a Northwest Atlantic cluster, a Northeast Arctic cluster, and a Northeast Atlantic cluster. We also observed a genetically distinct fjord population and a pattern of isolation by distance in the Northeast Atlantic. Our results contrast with the current management regime for this species in the Northeast Atlantic, as we found structure within some management areas. The study adds to the growing recognition of population structuring in marine organisms in general, and ﬁshes in particular, and is of clear relevance for the management of haddock in the Northeast Atlantic.


Introduction
A biologically based fisheries management regime does not always coincide with current management regimes, which are often based on factors such as administrative boundaries, oceanographic features, spatial distribution of fisheries, and the most important target species in the respective areas (Reiss et al., 2009;Kerr et al., 2017). Inconsistency between population structure and spatial management units may result in the overexploitation of specific population components, causing loss of local genetic diversity and depletion, or even extinction, of local populations, which may go undetected if managed as a single unit (Allendorf et al., 2008). Loss of genetic diversity may hamper a population's ability to adapt to changing environmental conditions (e.g. Reed and Frankham, 2003). Depletion or extinction of local populations may lead to ecosystem changes. Hence, failing to identify existing population structuring can lead to over-harvesting of one or more separate stocks in a fishery (Hauser and Carvalho, 2008), potentially causing a collapse (Cadrin, 2020) as exemplified by the Atlantic cod (Gadus morhua) fisheries in Newfoundland (e.g. Hutchings and Myers, 1994).
Haddock (Melanogrammus aeglefinus) is a commercially exploited demersal marine teleost, distributed along the coast and continental shelf on both sides of the North Atlantic Ocean (Figure 1). In the Northeast Atlantic, the species is distributed from the Bay of Biscay through the British Isles, the North Sea and Kattegat (but not in the Baltic Sea), along the Norwegian coast into the Barents Sea. In the Barents Sea, the distribution extends to Spitsbergen in the Northwest and Novaya Zemlya in the Northeast. Haddock is also present in the waters around Rockall, Faroe, and Iceland and to a limited extent off the southern Greenland coast (Bigelow and Schroeder, 1953). Haddock in the Northeast Atlantic is assessed and managed by ICES (International Council for the Exploration of the Sea) as seven separate units (1/2, 4/6a/20, 5a, 5b, 6b, 7a, and 7b-k; Figure 1). In the Northwest Atlantic, haddock is distributed from Cape May, NJ, to the Strait of Belle Isle, north of Newfoundland (Bigelow and Schroeder, 1953) and fisheries statistical units are set by NAFO (North Atlantic Fisheries Organization), recognizing six haddock stocks (3LNO, 3Ps, 4TVW, 4X, 5Z, and 5Y).
Identification of spawning areas and characterization of population genetic structure is important for robust stock assessment. Figure 1. The global distribution area of haddock and locations of the samples used in this study, the management areas and the known spawning grounds. Grey-shaded areas indicate the distribution area, light green-shaded areas indicate known spawning grounds, and coloured dots indicate the position where the samples were collected. The borders for the ICES and NAFO management units for haddock are shown with thick black lines. The management divisions/subdivisions are shown within thin black lines where the numbering refers to the respective divisions/subdivisions. Regions described in the text are written in red colour.
In the Northeast Atlantic, haddock spawning areas have been identified southwest of Iceland, off the Faroe Islands, in the deep waters west of Scotland, the northern North Sea, and offshore along the Norwegian coast ( Figure 1). Knowledge of exact spawning locations for Northeast Arctic haddock is uncertain, but three spawning areas are currently recognized along the slope between continental shelf and the Norwegian Sea ( Figure 1): off the Møre coast (at $61-64 N), offshore of the Lofoten Islands/ Røstbanken/Vesterålen, and on the western side of Tromsøflaket (Saetersdal, 1952;Bergstad et al., 1987). Spawning has also been reported by fishermen farther north along the Norwegian coast and in several northern fjords (Giaever and Forthun, 1999) and it has been suggested that spawning takes place along the entire outer continental shelf from $62 to $70 N (Saetersdal, 1952;Bergstad et al., 1987). Following spawning along the Norwegian coast, eggs and larvae drift with the currents in a northern and north-eastern direction.
In the Northwest Atlantic, spawning has been reported in the Gulf of Maine and two spawning aggregations with phenotypic differences, and otolith morphometrics, have been described on Georges Bank (Begg, 1998;Begg et al., 2001), but there is still uncertainty about connectivity of haddock in the Gulf of Maine and on Georges Bank (NEFSC, 2014). Distinct spawning areas are also reported at Browns Bank and adjacent banks off the southwestern coast of Nova Scotia (Hurley and Campana, 1989) and Emerald-Western Banks on the Nova Scotian Shelf. Farther north, in Newfoundland waters, spawning also takes place at the Grand Banks and St. Pierre Bank (Begg, 1998).
Several features of haddock reproductive biology potentially influence genetic population structure. As for most other gadoids, haddock is a highly fecund species with large population sizes and a high dispersal potential. Pelagic eggs and larvae go through a metamorphosis to a pelagic juvenile stage after 2-3 months (Fahay, 1983), and the juveniles may not settle for a further 2-3 months (e.g. Bailey, 1975). Following settlement, haddock may move similar distances to that of drifting pelagic stages up to maturity (Wright et al., 2010). Due to the long pelagic stages during which eggs and larvae can drift over deep-water regions that adult haddock normally do not cross, gene flow between populations are likely to take place during these early life-history stages.
Despite extensive drift of pelagic early life stages, genetic differentiation may develop and persist in marine environments where geographic or oceanographic barriers reduce gene flow. Surface circulation patterns around banks and shelf areas can retain eggs and larvae in natal areas and potentially play a functional role in maintaining population structuring. Retention of haddock eggs and larvae due to eddies and gyres has been reported on both sides of the Atlantic Ocean, such as the Rockall Bank and near the Faroe Islands in the Northeast Atlantic (Raitt, 1936;Fraser, 1958;Lee, 1974) and Gulf of Maine, Georges Bank, and along the Scotian Shelf in the Northwest Atlantic (e.g. Begg, 1998;Lough and Manning, 2001).
Both morphological variability and genetic variability have been used to define stocks of haddock. In the Northeast Atlantic, previous analyses have come to divergent conclusions regarding the population structuring in the Barents Sea and the connectivity towards the north Norwegian coast. It is still uncertain whether the Barents Sea is occupied by locally spawning haddock in addition to migrating haddock of more southern origin (Raitt, 1936;Giaever and Forthun, 1999) and whether there is any population structuring in haddock in north Norwegian waters at all (Saetersdal, 1952). As the Barents Sea and North Norwegian waters are among the most important fishing grounds for haddock, a deeper understanding of the population structure in the area is needed. Also, within the North Sea, which is another important fishing area, contrasting conclusions have been drawn regarding genetic differentiation. Child (1988) reported no significant genetic differentiation in two isozyme loci while Jamieson and Birley (1989) found distinct transferrin allele frequencies east and west of the Greenwich median and for the Rockall Bank. In the Northwest Atlantic, morphometric and meristic analyses (e.g. Grosslein, 1962) and microsatellite analyses (Lage et al., 2001) revealed discrete haddock stocks on the different banks. This finding is contrasted by mtDNA analyses, revealing no significant genotype differentiation among Northwest Atlantic banks (Zwanenburg et al., 1992), in addition to indications of haddock from several geographic regions spawning on Georges Bank (Purcell et al., 1996). Hence, the population structure and the connectivity between the Northwest Atlantic banks are still mainly unresolved and warrant further population genetic investigations.
There are few investigations of population structure for haddock, and the few studies that have been performed shows contrasting results on both sides of the North Atlantic. To shed light on the fine scale genetic structure within or across management areas, new investigations using genetic markers with higher resolution are needed. The present study had two primary aims: first, to develop a set of single-nucleotide polymorphism (SNP) loci that could be used for population genetic analysis in haddock and, second, to use these markers to provide the first pan-Atlantic analysis of population genetic structure in this commercially important species. We discuss the observed population structure in light of previous findings and show that it does not conform to the current management regime for haddock.

Development of SNP markers
For the initial SNP detection, high molecular weight DNAs were isolated from haddock caught off Bear Island (Barents Sea; N ¼ 8) and Møre (Norwegian coast; N ¼ 8) using Qiagen blood and tissue kit (Qiagen), according to the manufacturer's recommendations (see Supplementary Text). A double digest RAD (ddRAD) library was constructed (Supplementary Text) and sequenced on the Illumina MiSeq platform. Sequence reads were demultiplexed and SNPs were identified/scored using STACKS 1.47 (Catchen et al., 2013), followed by stringent filtering (Supplementary Text) before SNP locus primer design, amplification, and genotype calling were performed using the MassARRAY iPLEX Platform (Agena Bioscience).

Samples, DNA extraction, and genotyping
A total of 1329 haddock grouped into 19 samples (N ¼ 32-116) were collected in 2014 and 2017 (Table 1). All samples were collected in spring, at or just prior to spawning time, except for the Varangerfjord, Altafjord, and Nordland samples, which were collected in October. Samples from the two Lofoten locations, the Viking bank, North Sea central, Moray Firth and, Iceland were collected at, or near, known spawning grounds ( Figure 1). Haddock from Norwegian waters were collected by research vessels and by IMR's coastal reference fleet. From ICES areas 1 and 2 (north of 62 N), 393 individuals were collected at six locations (including 97 juvenile individuals from the Barents Sea). From Population genetic structure in Atlantic haddock ICES areas 4, 6a, and 20 (south of 62 N), 356 individuals from five locations along the Norwegian coast (including Skagerrak) were collected in addition to 274 individuals from four locations in the North Sea, which were collected as part of the Marine Scotland Science contribution to ICES quarter 1 bottom trawl survey. From Icelandic waters (ICES area 5a), 99 individuals were collected by the IMR reference fleet. Finally, 207 individuals from three locations, representing three different NAFO management areas, were sampled from the southern part of the distribution range in the Northwest Atlantic by the Northeast Fisheries Science Center (NEFSC) bottom trawl survey. Our sampling did not cover locations west of the British Isles or north of the Scotian Shelf but covered three important management areas in the Northeast Atlantic and three in the Northwest Atlantic ( Figure 1). DNA was extracted from ethanol-preserved muscle tissue, using the E.Z.N.A Tissue DNA kit (Omega Bio-Tek). All samples were individually genotyped for 167 SNPs in 6 multiplexes (Supplementary Table S1), with the MassARRAY iPLEX system. Manual inspection and genotype clustering of all SNPs were performed in TYPER 4.1 (Agena Bioscience). Genotyped SNPs with 135-bp flanking sequence are listed in Supplementary Table S1. Flanking sequences of all 167 SNPs were mapped against the published haddock genome (Tørresen et al., 2018), with the online version of BLAST (blast.ncbi.nlm.nih.gov/Blast.cgi). We also mapped the flanking sequences to the gadMor1 Atlantic cod genome assembly (Star et al., 2011), using the BLAT search tool (www.ensembl. org/Gadus_morhua/) to identify potential genes associated with the SNPs. We checked for linkage disequilibrium between SNPs using an r 2 threshold of 0.4 with PLINK 1.07 (Purcell et al., 2007).

Population genetic analyses
We used ARLEQUIN 3.5.1.3 (Excoffier and Lischer, 2010) to calculate observed and expected heterozygosity (H o and H e ) within each sample and to test for locus by locus departure from Hardy-Weinberg equilibrium (HWE) in each sample separately (100 000 iterations and a Markov Chain of 1 000 000). Correction for multiple testing [false discovery rate (FDR)] was performed in R (R Core Team, 2012), using QVALUE (Storey, 2002) with a q-value of 0.05 as a threshold. Minor allele frequencies (MAFs) were calculated using PLINK. Global F ST and weighted average F ST values between all pairwise samples, between management areas, and between identified sample clusters were calculated in ARLEQUIN with 10 000 permutations. Corrections for multiple testing were performed according to the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995) with 0.05 as a threshold for significance. We tested all 138 SNPs for signs of selection, using FDIST2 implemented in ARLEQUIN with 100 000 simulations. A distance-based neighbour-joining tree was constructed for all 19 samples, using DARTR 1.2.0 (Gruber et al., 2018). An MDS (multidimensional scaling) plot based on the pairwise F ST values was constructed with the isoMDS function in R-based MASS (Venables and Ripley, 2002).
We used the correlated allele frequency and admixture model with the locprior option in STRUCTURE 2.3.4 (Pritchard et al., 2000) to identify major genetic clusters in the dataset, performing 10 independent runs for each value of K (burn-in of 100 000 MCMC iterations followed by 1 000 000 MCMC iterations). STRUCTURE was first run on the full dataset, followed by a hierarchical approach, where additional runs within each of the identified groups were run separately. Delta K, mean ln P(K), and the K estimates MedMeanK, MaxMeanK, MedMedK, and MaxMedK were identified with STRUCTURESELECTOR (Li and Liu, 2018). Pie-charts, based on the population estimates of the ancestry coefficients from STRUCTURE, were constructed with LEA 2.4.0 (Frichot and Franc¸ois, 2015) with the add.pie function in mapplots.
We used R-based ADEGENET 2.1.1 (Jombart and Ahmed, 2011) to perform discriminant analysis of principal components (DAPC) on the full dataset with all SNPs, followed by a hierarchical approach where additional DAPC analyses were performed within sample groups identified in the previous DAPC. The results from the DAPC were consistent across a wide range of retained PCs. The number of PCs included in the analyses was chosen at the point where the cumulative variance relative to the number of retained PCs starts to plateau. The loading plots from the DAPC were used to identify the main SNPs that drove the genetic divergence among the samples. Mantel's test with 999 permutations was used to assess the presence and magnitude of isolation by distance (IBD) in the full dataset and different subsets of the Northeast Atlantic dataset (see Supplementary Table S2 for details), testing the relationship between pairwise linearized genetic distance, F ST /(1 À F ST ), and Euclidean geographic distance (in km), using the functionality implemented in DARTR 1.2.0. The Euclidean geographic distance matrix was calculated from the average latitude and longitude coordinates (Table 1) using PointDistance in R-based GDISTANCE 1.2-2 (van Etten, 2017).

Results
A total of 6.02 M haddock paired-end (p.e.) reads, with phred scores >20, were demultiplexed from the MiSeq run (194-583 K p.e. reads per individual; mean ¼ 376 K p.e. reads). Initial STACKS analysis revealed 4746 ddRAD loci that were scored for at least 75% of the 16 haddock samples. Of these 3742 loci (79%) were scored as being polymorphic, with up to seven SNPs per ddRAD locus. To select a practical and robust set of SNPs for high-throughput genotyping, the dataset was stringently filtered with respect to number of SNPs per ddRAD locus, presence of repeats, and length of the flanking sequence (Supplementary Text). A stringent final set of 231 ddRAD loci were selected for SNP assay design, and 167 SNPs were successfully designed for the MassARRAY iPLEX system (Supplementary Table S1). Of the 167 genotyped SNPs, 14 were scored as monomorphic, and 14 others had call rates <0.95 and were removed from further analyses. In addition, two SNPs were relatively tightly linked (Mae029 and Mae206; r 2 > 0.4) and consequently, SNP Mae206, which had the lower call rate of the two, was removed from further analyses. Our final dataset, therefore, consisted of 138 polymorphic SNPs with an average call rate >0.95 and these were screened in 1329 individuals of haddock from both sides of the Atlantic Ocean ( Figure 1, Table 1). Of the 138 SNPs, 131 mapped to the published haddock genome (Tørresen et al., 2018), 130 to the Atlantic cod genome (Star et al., 2011), and 63 of the latter were located within annotated Atlantic cod genes (Supplementary  Table S1).
When comparing the SNP frequencies between the samples, only 56 of 2278 comparisons were significantly out of HWE after FDR correction (q < 0.05) (Supplementary Table S3). Of these, 48 were due to heterozygote deficit and eight to heterozygote excess. Notably, the SNP Mae024 displayed a heterozygote deficit in all samples, accounting for 19 of the 48 observations of a heterozygote deficit, whereas all 8 significant heterozygote excesses occurred at the Mae054 SNP (Supplementary Table S3). The number of polymorphic SNPs varied between 105 in Altafjord and 126 in Lofoten north, both in the Northeast Atlantic, whereas average observed heterozygosity ranged from 0.186 in the Barents Sea sample in the Northeast to 0.226 in the Gulf of Maine sample in the Northwest Atlantic (Table 1). As the SNPs were designed from a relatively small panel of individuals (N ¼ 16) from the Northeast Atlantic, ascertainment bias, mostly affecting missing rare alleles, could be an issue when applied to samples from other areas such as Iceland and the Northwest Atlantic. However, no evident trend was detected in the number of polymorphic loci or the observed and expected heterozygosity within the different areas and all SNPs with low MAF in the Northwest Atlantic or Iceland also had low MAF in the Northeast Atlantic (Supplementary Table S3).
The global F ST value over all SNPs and among all samples was 0.0108. For the management areas, the pairwise F ST values were all significantly different from zero, except for area 5Y relative to 5Z and 4X in the Northwest Atlantic, and the largest difference (F ST ¼ 0.0216) was observed between ICES area 1/2 in the Northeast Atlantic and NAFO area 4X in the Northwest Atlantic (Supplementary Table S4). Samples managed in area 5a (Iceland) had smaller F ST differences relative to the Northeast Atlantic than to the Northwest Atlantic. Pairwise F ST values between samples within the ICES 1/2 area were generally low, and not significantly different from zero among the Barents Sea, Varangerfjord, Altafjord, and Lofoten north samples. However, these four samples were all significantly different from the Lofoten south and Nordland samples ( The pairwise F ST values, visualized using MDS score distributions (Figure 2a), gave a clear indication of genetic clustering of samples into three major geographic regions. Briefly, a Northeast Arctic cluster extends south to northern Lofoten and is represented by the four northernmost samples (Barents Sea, Altafjord, Varangerfjord and Lofoten north), a western cluster included all samples from the Northwest Atlantic (Georges Bank, Gulf of Maine, and Western Scotian Shelf), and a Northeast Atlantic cluster that occupied most of the Northeast Atlantic, north to Lofoten south and west to Iceland, and included 12 samples. The latter appeared to form an approximate north-south gradient with the relatively isolated fjord sample from Osterfjord as an exception. These relationships corresponded with the Neighbourjoining tree (Figure 2b) in that a Northeast Arctic group (Barents Sea, Altafjord, Varangerfjord, and Lofoten north), a Northwest Atlantic group (Georges Bank, Gulf of Maine, and Western Scotian Shelf), and a distinct Osterfjord sample were evident. In addition, the NJ tree placed the Iceland sample between the Northeast Arctic group and a Nordland/Lofoten south cluster.
Bayesian cluster analyses, as implemented in STRUCTURE, based on all samples, indicated that the most probable K value, Population genetic structure in Atlantic haddock  based on mean ln P(K), was 5 and delta K was 2 (Evanno et al., 2005), while MedMedK, MedMeanK, MaxMedK and MaxMeanK (Puechmaille, 2016) Table S6). Visual interpretation of the plots from K3 to K5 (Figure 3a) showed clear trans-Atlantic differentiation (reflecting delta K ¼ 2) and a distinction of the four Northeast Arctic samples (Barents Sea, Varangerfjord, Altafjord, and Lofoten north), consistent with K ¼ 3, as identified with the Puechmaille method. By plotting the population estimates of the ancestry coefficients from the STRUCTURE analyses (at K ¼ 3) as pie-charts onto a map, the trans-Atlantic divergence was evident, where the north-western samples were identified by a major green component as well as  Table S6).
The DAPC on the 138 SNPs from all samples (Figure 2c) uncovered geographic structuring and identified a Northwest Atlantic cluster, a Northeast Arctic cluster, and a Northeast Atlantic cluster. This pattern became even clearer when the samples were grouped into these three regions prior to the DAPC (Figure 2d). Following the DAPC on the full dataset, DAPC were run on the Northeast and Northwest Atlantic groups separately (Supplementary Figure S2a and b). In the Northwest Atlantic, a separation of all three samples were observed, which contrasts with non-significant F ST values between the Gulf of Maine sample  Figure S2c), partial separation between the samples was observed, even though the F ST values were not significantly different from each other. Finally, DAPC was run iteratively on the Northeast Atlantic dataset, where identified clusters were removed consecutively to Figure 3. Structure plots of the assignment probabilities in the full dataset consisting of all 19 samples for the K values 1-5 (a) and pie-charts for all samples projected onto a map (b). In STRUCTURE, a vertical bar represents each individual. The plots are presented for all SNP markers based on the combined results from ten independent STRUCTURE runs. The pie-charts are based on the population estimates of the ancestry coefficients from ten independent STRUCTURE runs, for K ¼ 3. Note the trans-Atlantic divergence, where the north-western samples are dominated by a major green component, and divergence in the Northeast Atlantic, where the northernmost samples are dominated by an orange component and southern samples by a blue component.
identify finer structuring in the dataset. First, the Northeast Arctic cluster was removed, then the Osterfjord sample, and finally a cluster consisting of Lofoten south, Nordland and Iceland (Supplementary Figure S2d-f).
To summarize F ST , STRUCTURE, and DAPC analyses, we observed a genetic break between populations across the North Atlantic Ocean and also genetic divergence among northern and southern samples along the Norwegian coast around the Lofoten peninsula.
The loading plots, based on the DAPC analyses, showed the contribution of each SNP to the observed genetic differentiation in the total material and the three separate areas identified by the DAPC (Supplementary Figure S3). The five SNPs Mae200, Mae005, Mae017, Mae032, and Mae049, primarily drove the observed differentiation in the full dataset. Of these five SNPs, Mae200 and Mae032 mapped to LG3 in the cod genome and were not located within any genes. Mae005 lies within the arfgef3 gene and Mae017 lies within a pik3ap1 gene analogue, and both SNPs mapped to LG15 in the cod genome (Supplementary Table  S1). In addition, SNP Mae033 (also on LG15 in the cod genome) contributed significantly to the differentiation in the Northwest Atlantic comparisons while the Mae049 did not contribute significantly to the differentiation in the Northeast Arctic samples. All of the five SNPs driving the observed differentiation in the full dataset showed an allele-frequency gradient from north to south along the Northeast Atlantic coast (Supplementary Figure S4).
The presence and magnitude of IBD were assessed by Mantel's test of linearized multilocus F ST values, against Euclidean geographic distances. The results were highly significant, both for the study area considered as a whole and for the Northeast Atlantic separately, and the latter displayed the higher correlation (r 2 ¼ 0.565, p < 0.001) (Figure 4, Supplementary Table S2). Removal of the Osterfjord sample, which is physically close and significantly differentiated from all other samples (based on F ST values), increased the correlation (r 2 ¼ 0.696, p < 0.001), while removing the more geographically distant Northeast Arctic samples reduced the correlation (r 2 ¼ 0.311, p < 0.001).

Discussion
Here, we present an analysis of a large number of SNP markers to assess the genetic population structure of haddock throughout most of its range. We observed a generally low yet, in most cases, statistically significant evidence of population genetic structuring, at levels comparable to that found in other gadoids in the same areas (e.g. Saha et al., 2015;Dahle et al., 2018). The most substantial genetic differences were found across the North Atlantic Ocean, and among samples north and south of the Lofoten Islands, in addition to a clearly defined Osterfjord sample. These observations suggest a population structuring that partly contrasts with the current management regime for haddock in the Northeast Atlantic region. We also observe evidence for substantial IBD among the Northeast Atlantic haddock populations.

The Northeast Arctic
The Barents Sea is a major feeding ground for adult haddock during the summer months in addition to serving as an important area for juvenile haddock (Olsen et al., 2010). In winter, adult haddock undertake long spawning migrations from the Barents Sea towards spawning grounds along the northern Norwegian coast (Bergstad et al., 1987) and the "choice" of spawning grounds seems to vary significantly between years, partly  (Table 1). dependent on density-dependent processes (Langangen et al., 2018). Previous analyses have come to divergent conclusions regarding the population structuring in the Barents Sea and the connectivity with northern populations along the north Norwegian coast. It is still uncertain whether the Barents Sea is occupied by locally spawning haddock in addition to migrating haddock of more southern origin (Raitt, 1936;Giaever and Forthun, 1999) and whether there is any population structuring in haddock in north Norwegian waters (Saetersdal, 1952). In our study, we found no significant genetic differences between the Barents Sea and the two samples from the Finnmark coast (Varangerfjord and Altafjord), nor between the more geographically distant Lofoten north sample. However, significant differences were observed between the Barents Sea and all southern samples, including the Lofoten south sample. These findings suggest a genetic break between haddock north and south of the Lofoten Islands, at around 68 N, inconsistent with the current management regime in which haddock north of 62 N are managed as a single unit. All samples in the Northeast Atlantic shows a pattern of IBD. The Lofoten area represents an area with increased genetic differentiation, but samples both north and south of this area are still part of a continuum where genetic distance increase with geographical distance. These observations are in line with findings in other gadoid species, such as Atlantic cod (Dahle et al., 2018), saithe (Pollachius virens) (Saha et al., 2015) and Pacific cod (Gadus microcephalus) (Cunningham et al., 2009;Drinan et al., 2018), where genetic differentiation increases with physical distance.
We genotyped juveniles from the Barents Sea that presumably were in their nursery areas. They represented offspring of unknown spawning populations, possibly originating southward, somewhere along the Norwegian coast. It is not possible to discriminate between the source population(s) for these individuals, based on present data, except that they have genetic similarities with adults from both Lofoten north, Varangerfjord and Altafjord. In Atlantic cod, eggs and larvae from Northeast Arctic cod spawning outside the Lofoten Islands and along the north Norwegian coast drift northward with ocean currents into the Barents Sea (Bergstad et al., 1987;Brander, 2005). A similar scenario is likely for haddock (Bergstad et al., 1987) consistent with our findings. In a simulation study, (Castaño-Primo et al., 2014) found that eggs and larvae from haddock spawning in Vestfjorden (corresponding to the Lofoten south sample in our study) drifted towards the coast and not into the Barents Sea as the spawning products from outer Lofoten (corresponding to Lofoten north in our study) and Tromsøflaket do. A modelling study by Myksvoll et al. (2014) also suggested retention of a large portion of Atlantic cod eggs spawned within Vestfjorden (corresponding to Lofoten south in our study), whereas only a minor proportion of eggs were transported northwards through small straits. These modelling studies are consistent with our observation of a genetic break north and south of the Lofoten Islands and with genetic similarity between spawning individuals north of Lofoten and juveniles in the Barents Sea. A plausible scenario is that haddock spawning in the outer parts of the Lofoten Islands consists of migratory individuals from the Barents Sea, whereas haddock spawning in the inner parts of Lofoten consists of coastal "resident" individuals. Such a pattern would mirror the situation of migratory and non-migratory individuals of Atlantic cod in the Lofoten area (Rollefsen, 1933;Hylen, 1964;Nordeide and Båmstedt, 1998;Brander, 2005). We conclude that haddock populations north and south of the Lofoten Islands belong to different biological populations and should optimally be managed as separate stocks.

The Northeast Atlantic
Most F ST values south of the Lofoten south sample were small, in many cases not statistically significant, but displayed an IBD pattern in agreement with earlier findings (Giaever and Forthun, 1999). A notable exception to this pattern refers to the genetically distinct sample from Osterfjord, which most likely reflects a partially isolated fjord population. Such coastal populations are likely to be relatively small and could be subject to local over-fishing, as would other coastal haddock populations that may have gone undetected by our geographically coarse-scaled sampling in this area.
Our dataset contains Shetland and Moray Firth samples located west of the Greenwich meridian and the Viking Bank and North Sea central samples located east of the Greenwich meridian, in addition to the Skagerrak sample on the eastern outskirt of the North Sea. A distinction was previously made for haddock east and west of the Greenwich meridian (Jamieson and Birley, 1989). Based on the F ST values, the Shetland sample was significantly differentiated from all other samples in our dataset, except for the adjacent Viking Bank sample. Noteworthy, the Shetland and Viking Bank samples come from deeper (>100 m) and colder waters than did the Moray Firth and North Sea central samples (<100 m). Jamieson and Birley (1989) describe Shetland haddock as a genetically stable population, yet belonging to the "west of Greenwich" group (consisting of west of Scotland, east of Scotland and Shetland) but distinctly differentiated from the "east of Greenwich" group (consisting of the Viking Bank and the Fisher Bank in the central North Sea). The identification of Shetland as genetically distinct is also observed in our dataset, except for the low genetic differentiation relative to the Viking Bank, implying significant connectivity between these two samples. In contrast to Jamieson and Birley (1989), we observe a separation between the Shetland and Moray Firth samples which belong to the "west of Greenwich" group. It is difficult to explain the apparent reproductive isolation of Shetland from other biological evidence, as otolith chemistry indicates that juvenile haddock from the Moray Firth can recruit north to Shetland (Wright et al., 2010). However, as haddock larvae and eggs tend to be advected south from Shetland (Heath and Gallego, 1998), it is possible that the otolith chemistry results reflect a return migration after settlement, suggestive of philopatry. Based on the observation of a distinct differentiation of the Shetland sample, it is plausible that Shetland haddock forms a separate population, which should be considered in future management of North Sea haddock.

The Northwest Atlantic
The haddock stock on Georges Bank was previously thought to be distinct from the nearby Gulf of Maine and Canadian stocks, based on tagging, meristic, life-history and genetic analyses (Lage et al., 2001). Furthermore, haddock from Georges Bank were not significantly differentiated from northern haddock stocks in one study (Zwanenburg et al., 1992), while significantly differentiated from Nantucket Shoals (to the south) and the Grand Banks (to the north) in another (Lage et al., 2001). There is some evidence of haddock movement and mixing among US and Canadian management units, including movement between the Gulf of Maine, the Bay of Fundy and the Western Scotian Shelf (McCracken, 1960;Grosslein, 1962) and movement between the Gulf of Maine and Georges Bank (Needler, 1931;NEFSC, 2014). Recruitment synchrony has also been observed between Georges Bank and the Gulf of Maine (Clark et al., 1982;NEFSC, 2014).
In the present study, we found small but significant F ST differentiation between Georges Bank and Western Scotian Shelf samples, but not between the Gulf of Maine and the former two. The DAPC plot shows distinct clusters for all of these three samples, implying some degree of population structuring. The genetic distinction of the Georges Bank sample may in part be explained by the deep Northeast Channel, separating Georges Bank from Browns Bank and the Scotian Shelf hindering adult movement between the areas. In addition, such a pattern could be strengthened by the inshore currents leading from the Gulf of Maine outwards to Georges Bank combined with the gyres running clockwise around Georges Bank that acts as a transportation and retention mechanism for planktonic eggs and larvae (Lough and Manning, 2001).

Management implications
We observed genetic structuring in haddock at multiple geographic scales, from the largest (trans-Atlantic) to the local (fjord) level. We identified three main genetic clusters, consisting of a Northeast Arctic cluster, a Northeast Atlantic cluster, and a Northwest Atlantic cluster. In addition, we observed a genetically distinct fjord population and genetic structuring within the North Sea. The study adds to the growing recognition of population structuring in marine organisms in general and fishes in particular (Hauser and Carvalho, 2008;Salmenkova, 2011) and is of relevance for fisheries management (Reiss et al., 2009).
While the spatial pattern of genetic structure we observed for haddock matches with the current management regime in many instances, there are some notable exceptions. At an ocean scale, ICES and NAFO areas on either side of the North Atlantic capture the trans-Atlantic divergence and some of the regional patterns such as an Icelandic component and little distinct genetic structure along the Norwegian coast south of 62 N. In contrast, we found evidence for population separation and finer-scale structures within the more comprehensively sampled Northeast Atlantic. Of particular interest is the separation between samples collected north and south of the Lofoten Islands, the observation of distinct fjord populations, as well as structuring within the North Sea, which warrants a more detailed investigation. Hence, the current management of haddock partly contrasts with biological units identified herein and should be reconsidered. The existence of fine-scaled differentiation in coastal waters (fjords) is more difficult to implement directly in management, but technical measures such as spawning closures could be considered where deemed necessary. Such minor stock units may play an important role as reservoirs of genetic variants and aid in resilience of the species during environmental change. Population structure of haddock in the sampled areas of the Northwest Atlantic is more equivocal, as haddock from the Gulf of Maine were not significantly different than those on Georges Bank or the Western Scotian Shelf. In addition, the genetic differentiation along the Norwegian coast followed a pattern of IBD, combined with a genetic break in the Lofoten area. Patterns of IBD have also been observed in Atlantic cod (Dahle et al., 2018) and saithe (Saha et al., 2015) along the Norwegian coast. As a pattern of IBD implies a continuum where genetic distance increases with geographic distance, it is difficult to infer the distance that causes a level of reproductive isolation that can lead to demographic difference. Hence, it is challenging to incorporate such information into current fisheries management regimes. However, the identification of such patterns is an argument for not treating the entire coast into one management unit for this species.

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.