Analysis of coastal cod (Gadus morhua L.) sampled on spawning sites reveals a genetic gradient throughout Norway’s coastline

Background Atlantic cod (Gadus morhua L.) has formed the basis of many economically significant fisheries in the North Atlantic, and is one of the best studied marine fishes, but a legacy of overexploitation has depleted populations and collapsed fisheries in several regions. Previous studies have identified considerable population genetic structure for Atlantic cod. However, within Norway, which is the country with the largest remaining catch in the Atlantic, the population genetic structure of coastal cod (NCC) along the entire coastline has not yet been investigated. We sampled > 4000 cod from 55 spawning sites. All fish were genotyped with 6 microsatellite markers and Pan I (Dataset 1). A sub-set of the samples (1295 fish from 17 locations) were also genotyped with an additional 9 microsatellites (Dataset 2). Otoliths were read in order to exclude North East Arctic Cod (NEAC) from the analyses, as and where appropriate. Results We found no difference in genetic diversity, measured as number of alleles, allelic richness, heterozygosity nor effective population sizes, in the north-south gradient. In both data sets, weak but significant population genetic structure was revealed (Dataset 1: global FST = 0.008, P < 0.0001. Dataset 2: global FST = 0.004, P < 0.0001). While no clear genetic groups were identified, genetic differentiation increased among geographically-distinct samples. Although the locus Gmo132 was identified as a candidate for positive selection, possibly through linkage with a genomic region under selection, overall trends remained when this locus was excluded from the analyses. The most common allele in loci Gmo132 and Gmo34 showed a marked frequency change in the north-south gradient, increasing towards the frequency observed in NEAC in the north. Conclusion We conclude that Norwegian coastal cod displays significant population genetic structure throughout its entire range, that follows a trend of isolation by distance. Furthermore, we suggest that a gradient of genetic introgression between NEAC and NCC contributes to the observed population genetic structure. The current management regime for coastal cod in Norway, dividing it into two stocks at 62°N, represents a simplification of the level of genetic connectivity among coastal cod in Norway, and needs revision. Electronic supplementary material The online version of this article (10.1186/s12863-018-0625-8) contains supplementary material, which is available to authorized users.


Background
The Atlantic cod (Gadus morhua L.) is a demersal fish found in the northern waters across the Atlantic. Due to its large size and historical high abundance, cod has formed the basis of some of the most commercially important fisheries in the Atlantic for centuries. However, a legacy of over-exploitation, potentially exacerbated by climate-driven changes in distribution, has left several cod populations depleted [1]. In turn, this has also resulted in the decline of many of the commercially significant fisheries [2,3]. The best-known example of this is the total collapse of the northern cod fishery off Newfoundland which fell from8 00,000 tons around 1970 to < 1000 tons by 1992 [4]. Cod is one of the best studied marine fishes, and as for many marine species where it has been investigated, statistically significant spatial population genetic structure has been observed (reviewed by [5][6][7]). The first studies of population genetic structure of cod were conducted in the 1960's using haemoglobin (e.g. [8][9][10][11]) and transferrin [12]. Shortly after, population genetic studies were performed using allozymes [13][14][15], mtDNA [16,17], Pantophysin (Pan I) [18][19][20][21] and microsatellites [7,[22][23][24][25][26][27][28]. More recently, single nucleotide polymorphisms (SNPs) [29][30][31][32][33][34] and full-genome sequencing [35,36] have been used to investigate the evolutionary relationships among cod populations. While not all studies on cod have revealed statistically significant population genetic structure [15,16,37], the majority have, and collectively, they reveal a species displaying population genetic differentiation throughout its native range.
Studies have revealed genetic differences between cod sampled on the eastern and western sides of the Atlantic [7,30,38], and between cod sampled in different ecosystems: in the Baltic vs. North Sea [25], North Sea vs. Norwegian coast [24], Canadian coastal vs. oceanic [38,39], coastal and oceanic in Gulf of Maine [40], and Greenland coastal vs. offshore with link to the Icelandic offshore cod [41,42]. Also, genetic differences have been observed between cod sampled in different areas within countries, including the North Sea [23], within UK [7], Iceland [43], North America [38,40], and Norway (North East Arctic Cod (NEAC) vs. Norwegian Coastal Cod (NCC)) [27,31,32,44]. Finally, even on small spatial scales, such as among neighbouring fjords in Norway, genetic differences have been observed [24,26]. These genetic studies therefore demonstrate that the species consists of multiple populations displaying varying levels of connectivity.
In addition to spatial genetic structure, genetic differences have been observed among cod "ecotypes" displaying stationary and long-distance migratory behaviours. For example, large genetic differences have been reported between NCC and NEAC, which display stationary and migratory behaviours respectively [29,31,32]. The same has also been observed between migratory North Sea cod (NSC) and coastal cod in southern Norway [45], between migratory and stationary in Iceland [31,46], and Canada [47]. Furthermore, the population genomic approaches [48] utilised in many of the most recent studies have also revealed large genomic inversions, notably in linkage groups 1, 2 and 7 that are assumed to be responsible for much of this strong divergence between ecotypes [29,32]. However, genomic islands of divergence, probably also caused by genomic inversions, have also been reported between low-salinity adapted cod in the Baltic, and the North Seas [49].
Despite the considerable number of population genetic studies on cod, there are still a number of remaining questions. This is the case for Norway, which is the country with the largest remaining commercial catch of cod in the Atlantic. In Norway, coastal cod are currently divided into two management stocks, defined as NCC north of 62°, and coastal cod south of 62°(from hereon we refer to both of these components as NCC for simplicity). While the NEAC stock is currently at a high level [50,51], NCC is depleted and a rebuilding plan has been recommended by ICES [51,52]. Importantly, while studies of NCC sampled in neighbouring fjords have been conducted in limited geographic areas [24,26], population genetic structure has not been investigated in detail along the entire Norwegian coastline spanning.
In 2002, a project was initiated to map the population genetic structure of NCC along the entire Norwegian coast. This included sampling cod from 55 spawning locations spanning the Russian to Swedish borders with Norway ( Fig. 1). At the same time, biological data and brood-stock fish were sampled [53], some of which have formed the basis of common-garden experiments [54,55]. Here, we present the results of the genetic analysis based upon data from microsatellite loci and the Pan I locus.

Sampling
In the period from 2002 to 2007, 4422 cod were collected from 55 locations along the coast of Norway (Fig.  1). Eleven of those sites were sampled more than once during the six-year period. Sampling was conducted during the spawning season (late winter and early spring) by taking samples from the catch of local fishermen. All of these fish were collected at known spawning sites. Most of the cod were sampled using gill nets, although for some individuals around the Lofoten islands, demersal trawls were also used. Biometry data such as length, weight and sex of all individuals were recorded, and the gonads were visually inspected to determine the stage of maturation. Fin clips from all the individuals were taken and stored in 96% ethanol prior to DNA extraction.
Although genetic analyses have also been used for identification of NEAC and NCC in Norwegian fisheries [56,57], in routine Norwegian stock assessment north of 62°, cod are assigned to NEAC and NCC by otolith category. Otoliths from all the fish sampled in this study were analysed to determine age, age at first time of spawning and the number of spawning periods, as well as to identify individuals as NEAC or NCC. Otolith type 1 and 2 is assigned NCC whereas 4 and 5 are assigned to NEAC, as described for the first time by Rollefsen [58] and modified by Berg & Albert [59]. Assignment otolith category is in strong agreement with results from genetic markers [27,59], and particularly Pan I. At this marker, NCC show a high frequency of the Pan I AA genotype and NEAC show a high frequency of the Pan I BB genotype [26]. In the present study, all individuals with otoliths belonging to types 4 and 5, were removed from the majority of the statistical analyses in order to exclude any NEAC from the samples, which could influence estimates of NCC population genetic structure (76 such individuals were removed in total).

Genotyping
DNA extraction was performed using the Qiagen DNea-syH96 Blood & Tissue Kit; each of which contained two or more negative controls. All 4346 individuals were genotyped at six microsatellite loci: Gmo2, Gmo3, Gmo34, Gmo35, Gmo132 and Tch11 [60][61][62], together with the Pantophysin locus Pan I [18]. In addition, a subset of 1295 individuals from 17 of the sites were genotyped for an extra set of nine microsatellites (bringing the total microsatellite loci up to 15): GmoC18, GmoC20 [63]; GmoG13, GmoG18 [64]; GmoG25, GmoG40, GmoG43, GmoG45 [64], and Tch22 [62]. Thus, the present study includes two overlapping data sets: one with all 4346 individuals genotyped for six microsatellites and Pan I (55 locations, hereon referred to as Dataset 1), and the other with a sub-set of individuals 1295 individuals genotyped for 15 microsatellites and Pan I (17 locations, hereon referred to as Dataset 2). The PCR conditions are available from the authors upon request. PCR products were analysed on an ABI3130XL sequencer (Applied Biosystems), whereas Pan I was genotyped on 2.5% MetaPhore gels. Microsatellite alleles were scored using GeneMapper v4.0 (Applied Biosystems).

Statistical analysis
Statistical analyses were conducted separately for microsatellites and Pan I. The total number of alleles, and allelic richness, were both calculated with MSA [65], whereas observed (H O ) and unbiased expected (uHe) heterozygosity, inbreeding coefficient (F IS ) and deviations from the expected Hardy-Weinberg distribution, were computed with GenAlEx [66]. Possible linkage (LD) between all locus pairs per population was tested using the program GENEPOP on the web [67] with significance based on the Markov chain method with 10,000 dememorizations, 20 batches and 5000 iterations per batch. Effective population size (Ne) per sample was computed using LDNE [68], implementing the threshold values of lowest allele frequency of 0.05 and 0.01. Where applicable, signification was corrected by multiple comparisons by sequential Bonferroni correction [69] implemented in the calculator developed by Justin Gaetano (2013).
To test if loci deviated from neutrality, outlier analyses were conducted with LOSITAN [70] under a stepwise model and the following settings: 1000000 simulations, 99.5% confidence interval, forced mean F ST , and with a 0.01 false discovery rate. Genetic differentiation among sampling sites was tested using the Analysis of Molecular Variance (AMOVA) as well by pairwise F ST. Both analyses were implemented in ARLEQUIN v.3.5.1.2 [71] and significance was calculated after 10,000 permutations. Likewise, hierarchical AMOVA was conducted by pooling populations according to geographic areas depicted in Table 1.
Several parameters such as the number of alleles, allelic richness, allele frequency, Ho, uHe and pairwise F ST were tested for trends in the geographic north-south gradient using the non-parametric Kendall measure of rank correlation [72], which measures the similarity of the orderings of the data when ranked by north-south gradient or by the value of the variable tested [73], and implemented in the R Package 'Kendall' [74].
Population genetic structure and connectivity were investigated using two approaches. First, through BAR-RIER 2.2 software [75], which aims to reveal genetic barriers between populations using Monmonier's [76] algorithm. The significance for this analyses was tested by bootstrapping 1000 matrices computed with Nei's DA genetic distance [77]. In addition, STRUCTURE v. 2.3.4 [78] was used to identify genetic groups under a model assuming admixture and correlated allele frequencies without population information. STRUCTURE analyses were parallelised with the program ParallelStructure [79] to speed up computation time. After multiple runs of both programs, and probably due to the nature of the genetic structure observed (see results), the software failed to find clear barriers (data not presented for BAR-RIER). Therefore, in order to complete some of the population genetic analyses (i.e., AMOVA as described above), we subjectively chose six geographic regions for some of the hierarchical analyses, which are not advocated as Management Units.
In the six microsatellites common to both datasets, the total number of alleles observed per locus showed extremely similar values, despite that the number of individuals genotyped displayed a 3-fold difference between the two data sets (Table 3). In five of the six microsatellite markers used in Dataset 1, global F ST per locus significantly differed from zero (Table 3). In eight of the fifteen microsatellite markers used in Dataset 2, global F ST per locus was significantly different from zero. The locus Gmo132 clearly displayed a much higher global F ST than any of other loci in both data sets, and was reported to be under directional selection by LOSITAN (P = 1.0).
Global F ST over all microsatellites was low but statistically significant in both datasets (Dataset 1, F ST = 0.0075, P < 0.0001, and Dataset 2, F ST = 0.0042, P < 0.0001). Global F ST did not change when removing the individuals showing the genotype Pan I BB, which is the dominating genotype in NEAC and observed in very low frequency in NCC. This meant excluding 104 and 39 individuals Table 1 Summary statistics for Dataset 1 (66 samples genotyped at six microsatellites): geographic region the sampling site belongs to, number of sample, sample name, total sampling size (i.e. NEAC and NCC), sampling size for NCC (i.e. individuals with otoliths 1 and 2); % of NEAC in the sample (i.e. individuals with otoliths 4 and 5); number of alleles; AR, allelic richness (based on a 25 diploid individuals); observed (Ho) and unbiased expected heterozygosity (He), inbreeding coefficient (F IS ), number of deviations from Hardy-Weinberg equilibrium (HWE) and from Linkage Disequilibrium (LD) at α = 0.05 The genetic matrix for Dataset 1 revealed that 2.7% of the pairwise F ST comparisons within the seven geographically-determined regions were significant (which equates to 19% of the total combinations), in contrast with 50.5% of the pairwise comparisons among regions (i.e. 59% of the total, Additional file 2: Table S1). However, when using Bonferroni correction for multiple comparisons, the corrected critical value dropped to 0.00005, which involved a decrease from 2.7 to 0.1% of significant pairwise F ST within regions, and from 50.5 to 20.5% among regions. Similarly, the genetic matrix for Dataset 2 revealed that 4% of the pairwise F ST comparisons within the seven geographically-determined regions were significant (i.e. 24% of the total combinations), in contrast with the 55% found among regions (65.6% of the total, Table 4). The corrected critical value dropped to 0.0003 after Bonferroni correction resulting in no significant comparisons within regions, and a drop of 55 to 27% of significant pairwise F ST among regions. These data suggest isolation by distance, in agreement with the fact that pairwise F ST values were found to be strongly correlated with the ordering of the samples in the Table 1 Summary statistics for Dataset 1 (66 samples genotyped at six microsatellites): geographic region the sampling site belongs to, number of sample, sample name, total sampling size (i.e. NEAC and NCC), sampling size for NCC (i.e. individuals with otoliths 1 and 2); % of NEAC in the sample (i.e. individuals with otoliths 4 and 5); number of alleles; AR, allelic richness (based on a 25 diploid individuals); observed (Ho) and unbiased expected heterozygosity (He), inbreeding coefficient (F IS ), number of deviations from Hardy-Weinberg equilibrium (HWE) and from Linkage Disequilibrium (LD) at α = 0.05 (Continued)   (Fig. 3a-b). Hierarchical AMOVAs were conducted for datasets 1 and 2 using the aforementioned geographic approach to define regions. In both cases, the three levels of division (among regions, among populations within regions and within populations) showed significant values ( Table 5). The differentiation observed among regions was higher than the differentiation observed among populations within regions, however, most of the observed genetic variance was hosted within populations (> 99% according to the microsatellites in Datasets 1 and 2).
Hierarchical AMOVAs per locus showed that the locus Gmo132 was the only one significant at all levels of division in both datasets (Table 6). This marker, which was depicted to be under directional selection by LOSITAN analysis, showed the highest differentiation within populations (average F ST of 0.034). Locus Gmo34 showed significant F CT and F ST (average F ST = 0.005) in all analyses. The extra set of nine microsatellites used in Dataset 2 showed one locus also significant at all levels (Gmo45), but also revealing a weak degree of structuring (average F ST = 0.006). Interestingly, the most frequent alleles for both Gmo34 and Gmo132 displayed a strong and highly significant trend in the north-south gradient (Fig. 4). No such trend was observed for the most common alleles in Gmo45 (data not presented). Allele 96 from locus Gmo34 ranged from 0.739-0.375 (τ = − 0.448, P = 1.174 e-07) and allele 115 from locus Gmo132 ranged from 0.681-0.105 (τ = − 0.626, P = 1.111 e-13) in the direction north to south.
Looking at the results for the Pan I locus, that was analysed for all samples, the three different genotypes  Effective population size (Ne) was calculated using two values of lowest allele frequency (0.05 and 0.01). Sites have been ordered from north to south and information about the geographic region where they are placed Table 3 Summary information for the markers used in this study: linkage group (LG) to which they belong to, whether the linkage group displays a known inversion, global F ST (and associated P-value after 10,000 permutations) and number of alleles. n/a stands for "not available" Tch11 n/a n/a 0.0020 (P = 0.0003) 0.0015 (P = 0.0788) 29 27 GmoC18 n/a n/a 0.0002 (P = 0.4351) 19 GmoC20 n/a n/a 0.0000 (P = 0.5175) 21 GmoG13 n/a n/a 0.0021 (P = 0.1768) 18 GmoG18 n/a n/a ***** 5 GmoG25 n/a n/a ***** 63 GmoG40 n/a n/a 0.0035 (P = 0.1111) 16 GmoG43 n/a n/a 0.0011 (P = 0.1241) 30 GmoG45 n/a n/a 0.0056 (P = 0.0000) 47 Tch22 n/a n/a ***** 6 Pan I 1 Yes 0.107 (P < 0.0001) 2 LG is in accordance with nomenclature from Hubert et al. [48]. Microsatellites were blasted against the gadMor2 assembly [101] by Per Erik Jorde (pers. comm) P-values in boldface type were significantly different from zero displayed clearly distinct frequencies throughout the NCC samples (Fig. 5). Genotype AA was the most frequent and experienced an increasing trend towards the south, whereas genotype BB, present in low or very low frequencies in all populations, showed a decreasing southwards trend. The frequency of allele Pan I A revealed a negative trend towards the south (τ = 0.23, P = 0.0068, Fig. 4). Overall genetic structure for Pan I locus was highly significant (F ST = 0.107, P < 0.0001) and hierarchical AMOVA was also significant at the three levels of grouping (Table 5). Likewise, pairwise F ST for Pan I showed a highly significant north-south trend (τ = 0.871, P < 2.22 e-16). In order to investigate the potential admixture between both types of cod, NEAC was used as an outlier group for NCC in STRUCTURE. These analyses showed no clear clustering of populations in the full dataset genotyped at 6 microsatellites (Fig. 6a) in agreement with BARRIER (results not shown). However, the two outlier loci (Gmo32 and Gmo134) displayed a north-south gradient of admixture (Fig. 6b) whereby NCC showed greater genetic similarity to NEAC in the north than in the south. A similar trend could also be observed in the restricted dataset genotyped at 15 microsatellites (Fig. 6c).

Discussion
Norway is the country with the largest remaining commercial catch of Atlantic cod, and this is the first study to investigate population genetic structure of coastal cod (NCC) along the entire Norwegian coastline. Following the analysis of > 4000 cod sampled from 55 locations, we obtained the following main results: 1. Statistically significant population genetic differentiation was revealed among most of the samples, 2. The observed genetic differentiation followed a pattern of isolation by distance along the north-south gradient, without any clear "breaks", 3. No distinct change in genetic variation (i.e. allelic diversity or heterozygosity) was observed among the samples in the north-south gradient. Based upon these results, we conclude that NCC displays statistically significant population genetic structure along the Norwegian coastline. Consequently, these results demonstrate that the current management regime, dividing coastal cod in Norway into two management groups, one north and one south of 62°, represents an over-simplification of the true level of population genetic structure, and as such, requires re-evaluation.
As detailed in the introduction, Atlantic cod is one of the best-studied marine fishes, and a large number of population genetic studies have revealed a species displaying extensive population genetic structure (reviewed by [5][6][7]). While several of the newest studies in this species have utilised genomics methods to identify and characterise population genetic structure (e.g. [29,32,38,41,45]), our study based upon microsatellites and Pantophysin, has nevertheless made a significant contribution to our understanding of population genetic structure in this species, especially in Norway. First the biological, and then the management implications of our results are discussed below.

Population genetic structure of NCC
The primary aim of the present study was to investigate population genetic differentiation among coastal cod sampled from the entire Norwegian coastline, which stretches over 2500 km from the south in Oslofjord, to the north in Porsangerfjord (Fig. 1). In some areas of Norway, in particular the north, NEAC are found using the same spawning areas as NCC [80,81]. Consequently, some NEAC were inadvertently collected during the sampling process, especially in the north (Table 1). In   GmoG18 ******* ******* ******* ******* ******* ******* GmoG25 ******* ******* ******* ******* ******* ******* Tch22 ******* ******* ******* ******* ******* ******* P-values in boldface type were significantly different from zero Fig. 4 Frequency of the most common alleles within loci Gmo132, Gmo34 and Pan I in Dataset 1. Populations are ordered from north to south and the first one, called NEAC, consists of the 76 Pan I BB individuals with otolith types 4 and 5 that were purged from the dataset. Alleles Gmo132_115 and Gmo34_96 experienced a highly significant negative trend southwards (τ = − 0.626, P = 1.11 e-13, and τ = − 0.448, P = 1.17e-07, respectively), whereas for Pan I _A, this tendency was reversed (τ = 0.23, P = 0.007). The allele frequency for NEAC assessed in the Barents Sea was 0.745 for Gmo132_115, and 0.959 for Gmo34_96 [28] order to eliminate the potential bias caused by NEAC in the present analysis of NCC population genetic structure, we used otolith categories 4 and 5 to exclude any potential NEAC from our biological samples prior to genetic analysis [58]. The Pan I BB genotype is nearly fixed for NEAC [26], and is almost completely diagnostic between NCC and NEAC [18,26,27]. In addition, Pan I genotype and otolith type both clearly differentiate between NEAC and NCC [27,82]. Despite purging otolith category 4 and 5 individuals from the current data set, 104 fish displaying an otolith structure typical for coastal cod (i.e., category 1-2) with the Pan I BB homozygote genotype (i.e. characteristic for NEAC) remained. Importantly however, and in the context of the present study, both including and excluding these 104 BB homozygote individuals did not influence the overall picture of population genetic structure. Therefore, it is concluded that we have effectively excluded NEAC from these analyses, and that the differences reported here primarily reflect genetic differences between NCC populations along the Norwegian coastline. Given that genetic differentiation did not show any clear genetic groupings, nor clear breaks in population connectivity, we suggest that NCC populations belong in a genetic gradient characterised by isolation by distance. This pattern is potentially driven by the generally limited migratory behaviour of NCC as has been seen in tagging experiments ( [24,81,83], but see [84]), spawning site fidelity [41,85], and potentially retention of eggs in certain fjord areas [83,86].
A distinct north-south cline in the frequency of the allele 115 for the microsatellite Gmo132 was observed here (Fig. 4). This pattern may be explained using two different scenarios: 1) Since this genetic marker has been tagged as an outlier in previous genetic analyses [26,87], the observed north-south cline at this locus might have been shaped under the influence of environmental factors. In our data however, significant allele frequency differences for this locus also existed among samples within the same region, thus it is unlikely that selection alone has created the observed north-south allelic cline (see [28]). 2) Although NEAC were removed from most of the analyses of the present study, this does not preclude the possibility that genetic introgression and admixture between NEAC and NCC contributes to the observed pattern of population genetic structure in NCC. Supporting this suggestion is the fact that the frequency of the 115 allele in Gmo132, in all of the northern samples investigated here, displayed a frequency around 0.5-0.6, which is very similar to the frequency of this allele in NEAC [~0. 7 28]. Therefore, admixture between NEAC and NCC, primarily in the north and following a decreasing gradient of presence towards the south, may contribute to the pattern of population genetic structure observed in NCC. This suggestion is consistent with simulations [88] illustrating that introgression from a genetically distinct source (e.g., NEAC) may generate gradients in allele frequencies along a geographic axis originating at the edge of the contact zone (i.e., north in the current study). These "tails of introgression" may extend to large distances beyond the contact zone itself (i.e., towards the south in the cod case). Importantly, the gradient only appears where dispersal and gene-flow are spatially limited. In our case, the limited gene-flow among NCC populations in Norway results in a gradient in the allele frequency, which starts in the north due to introgression of NEAC.
Population genetic differentiation was observed among samples of NCC throughout all regions of Norway, including among samples from southern Norway where NEAC was hardly present (as determined either by Pan I genotype or otolith category). In the far southern region, it is unlikely that NEAC has played a direct role in shaping the evolutionary relationships among NCC populations (although it is possible that the tail of potential introgression from NEAC in the north could still extend to the far south and thus play a minor role). However, the interaction between migratory North Sea cod and NCC populations in the south of Norway also remains a potential source of influence on NCC structure in this region [89]. Clearly, the evolutionary relationship between NEAC and NCC, and North Sea cod and NCC, needs to be resolved in order to evaluate the degree of influence migratory ecotypes have on the population genetic structure of NCC throughout Norway.
Current evidence suggests that the long-distance migratory "ecotype" of Atlantic cod, including NEAC, was derived from the stationary "ecotype" [29], and that this divergence occurred prior to the split between the Northeast and Northwest cod populations [90], which  Table 1 has been estimated to have occurred approximately 100,000 years ago [91,92]. As stated above, NEAC and NCC exhibit large genetic and genomic differences between them, and while several potential mechanisms have been proposed [80,93,94], the ecological processes leading to and maintaining segregation of these two ecotypes, and the degree of genetic exchange between them in both time and space, remain unresolved. The (~2%) BB homozygotes observed in the present data set, that displayed an otolith structure typical for coastal cod, could fit in one of these scenarios: 1 -"true" NEAC that for some reason did not migrate to and from the Barents sea as is characteristic for NEAC, 2hybrid and or admixed individuals between NEAC and NCC that did not migrate to and from the Barents sea, possibly due to only carrying one set of chromosomes with the "inverted supergene" and therefore displaying reduced propensity for long-distance migration, 3 -"true" NCC cod that display the BB genotype as the B allele is observed in NCC although at low frequencies. Quantifying admixture between NEAC and NCC remains a challenge that even recent papers focussing on the genome-wide differences between NEAC and NCC have not completely resolved [29,32]. Pan I lies within one of the inverted regions of chromosome 1, which together with inversions in linkage groups 2 and 7 are responsible for nearly all of the genomic divergence between NEAC and NCC [29,32]. As these inversions block recombination in their respective locations on the genome for NEAC, this means that if NEAC and NCC hybridise, their offspring will contain one copy of the inverted parts of linkage groups 1, 2 and 7, and one copy of the ancestral collinear form. We suggest that a genomic analysis of BB homozygotes displaying otolith categories 1-2 may provide an important resource in resolving the issue of hybridisation between NEAC and NCC, and the degree to which this may or may not influence population genetic structure of NCC [95]. Finally, the discrepancy in the allele frequencies in Pan I (linkage group 1) and Gmo132 (linkage group 7) reported between NEAC and NCC along the north-south gradient require further investigation.

Management implications
Applying the same management strategy to multiple populations or stock components that vary in their abundance and/or resilience to exploitation inevitably results in overfishing and likely collapse of the weaker component [96,97]. The extent of population genetic structure revealed in the present study, irrespective if it is influenced by a gradient of NEAC admixture or not, and divergent selective forces or not, strongly suggests that the current division of coastal cod in Norway, above and below 62°north, will not be sufficient to ensure sustainable management of NCC throughout Norway. Our data clearly illustrate population genetic structure within all areas of Norway, and as such, needs to be taken into consideration when renewing management plans. This conclusion is supported by earlier studies of population genetic structure in the south of Norway, where temporally-stable population genetic differentiation has been observed between neighbouring fjords over relatively small distances, effectively demonstrating limited connectivity at least between some fjord systems [24,98,99]. Additionally, earlier data from the north of Norway, revealing genetic differences between coastal cod sampled in fjords, supports the main results from the present study [14,19,26,28]. It is nevertheless acknowledged that the lack of any clear "breaks" in population connectivity along the Norwegian coastline means that identifying appropriate management units is challenging.
An additional management contribution from the present analyses is the apparent lack of any clear differences in genetic variation, as reported by numbers of alleles, allelic richness, heterozygosity or the effective population size (Fig. 2-4; Additional file 3), in the north-south gradient. In all areas of Norway, NCC has been documented in decline since 1990s, to which overfishing has probably played a significant role. However, in the southern regions, populations also appear to be influenced by climate-driven recruitment challenges as has been illustrated by the beach-seine survey conducted in this region since 1919 [100]. Despite this, based upon the samples analysed here, there are no indications of severe genetic bottlenecks in any of the samples, and no clear differences in genetic diversity estimates between the north and south.

Conclusions
Norwegian coastal cod displays statistically significant population genetic structure along the Norwegian coastline, which follows a genetic gradient characterised by isolation by distance. Although not fully resolved, we suggest that genetic introgression and admixture between NEAC and NCC, most in the north and least in the south, together with limited gene-flow among NCC populations, may contribute to the observed population genetic structure of NCC. In turn, we conclude that the current management regime in place, dividing coastal cod in Norway into two management groups, north and south of 62°, represents an over-simplification of the true level of population genetic structure, and as such, requires re-evaluation.

Funding
The study was financed by the Norwegian Ministry for Industry and Fisheries.

Availability of data and materials
Supporting information is available in the additional files and further information is available from the authors on request.
Authors' contributions KEJ conceived the study, GD conducted the microsatellite and panthophysin analyses, and drafted the paper. FB, MQ and JIW conducted the statistical analysis, KEJ, TJ and AA provided background material to the study. KAG led the process of data interpretation and the final version of the manuscript, with scientific contributions from all other authors. All authors read and approved the final manuscript.
Ethics approval and consent to participate Sampling (fin clippings) where conducted from dead fish caught in the legal commercial cod fisheries along the Norwegian coast.

Consent for publication
There is no material in the manuscript (personal details, images or videos) that would require any consent to publish.

Competing interests
The authors declare that they have no competing interests.