Combination of genetics and spatial modelling highlights the sensitivity of cod (Gadus morhua) population diversity in the North Sea to distributions of fishing

Michael R. Heath1*, Mark A. Culling2†, Walter W. Crozier3, Clive J. Fox4, William S. C. Gurney1, William F. Hutchinson5, Einar E. Nielsen6, Martha O’Sullivan7§, Katharine F. Preedy1‡, David A. Righton8, Douglas C. Speirs1, Martin I. Taylor2}, Peter J. Wright7, and Gary R. Carvalho2 Department of Mathematics and Statistics, University of Strathclyde, Glasgow G1 1XH, UK School of Biological Sciences, Bangor University, Bangor LL57 2UW, UK Agri-Food and Biosciences Institute, Newforge Lane, Belfast BT9 5PX, UK Scottish Association for Marine Science, Scottish Marine Institute, Oban PA37 1QA, UK Department of Biological Sciences, The University of Hull, HU6 7RX, UK Technical University of Denmark, National Institute of Aquatic Resources, DK-8600 Silkeborg, Denmark Marine Scotland Science, Marine Laboratory, 375 Victoria Road, Aberdeen AB119DB, UK Centre for Environment, Fisheries and Aquaculture Science, Pakefield Road, Lowestoft NR33 0HT, UK *Corresponding Author: tel: +44 141 548 3591; fax: +44 141 548 3345; e-mail: m.heath@strath.ac.uk


Introduction
Differences in genotype frequencies between population units account for a proportion of intraspecific genetic diversity, which is essential for the resilience of populations during periods of environmental change (Hilborn et al., 2003;Schindler et al., 2010). For exploited marine fish, the geographic boundaries which define barriers to gene flow are difficult to discern and do not necessarily coincide with the spatial configuration of management domains. Where management regions span multiple non-interbreeding subpopulations, this leaves less productive populations vulnerable to extirpation through overfishing (Hutchinson, 2008;Waples et al., 2008;Reiss et al., 2009). These concerns, among others, are behind proposed changes in fisheries management including the establishment of closed areas (Ralston and O'Farrell, 2008;Botsford et al., 2009;Armstrong et al., 2013). However, despite evidence of geographic discontinuities in gene flow within a range of fish species (Hauser and Carvalho, 2008), genetic studies have had little impact on fish stock assessments or the development of spatial models of fish population dynamics (Reiss et al., 2009). Partly, this maybe because early studies struggled to resolve statistically significant genetic structure with the markers then available (e.g. allozymes and microsatellite DNA). In addition, efficient numerical techniques for modelling the spatial demography of fish populations are a relatively new development and still not widely available.
Avariety of methods have previously been used to resolve spawning groups of cod in the northeastern Atlantic. Conventional and data storage tags deployed on cod caught and released at spawning locations have revealed that some fish roam over hundreds of kilometres, but return to spawn in the general area where they were tagged (Wright et al., 2006a;Righton et al., 2010), while others roam very little (Righton et al., 2010). Microchemical analyses of the juvenile component of adult otoliths (Elsdon et al., 2008) suggest that young of the year that settle in nursery areas have undergone a period of oceanographic drift as eggs and larvae and may have originated from distant origins (Wright et al., 2006a, b;Heath et al., 2008). Analyses of allele frequencies at microsatellite DNA loci in mature cod caught at spawning locations have revealed distinct groups of fish in the northeastern and the central North Sea (Hutchinson et al., 2003) and in the Irish and Celtic Seas (Hutchinson et al., 2001), though doubts have been raised concerning the significance of the results due to the inclusion of loci under selection (Nielsen et al., 2009). However, little of this information is recognized in management advice or measures to regulate cod fisheries in the region.
Existing management procedures for cod fisheries around the British Isles are based on four regions, the definitions of which long predate the availability of any information on genetics of population structure: North Sea, Eastern Channel, and Western Skagerrak (ICES Divisions IV, VIId, and IIIa); West of Scotland (Division VIa); Irish Sea (Division VIIa), and the Celtic Sea and Western Channel (Divisions VIIe -k; Figure 1). Total allowable catches (TACs) and national landing quotas for cod in each of these regions are set annually based on advice from scientific assessments and projections (e.g. ICES, 2010a) which assume that the stock in each region is self-contained with respect to recruits and is not subject to immigration or emigration. Evidence of genetic differentiation among fish from different assessment regions would uphold the spatial structure of management. Failure to detect differentiation would not necessarily undermine management, but might indicate the rates of exchange which could be demographically significant. On the other hand, genetic differentiation within an assessment region should be a cause for concern, since it would indicate the existence of small-scale populations which would be inadequately protected from overexploitation by the current management regime. So far, the evidence from microsatellite DNA analysis has neither dismissed nor upheld the existing management regions for cod, but has provided some indications of structure at smaller spatial scales.
In the study described here, we gathered data on singlenucleotide polymorphisms (SNPs) markers in Atlantic cod (Gadus morhua) sampled in the waters around the UK to identify the geographic barriers to interbreeding and how these compared with the management regions. SNPs have significantly higher allelic richness than the previously used microsatellite markers, resulting in greater discriminatory power at equivalent sample sizes. Using this information, we developed a spatial model of cod dynamics centred on the North Sea which resolved two reproductively isolated but nevertheless interacting subpopulation units. We then used the model to hindcast past changes in cod abundance and distribution and predict the consequences of varying spatial distributions of fishing.

Material and methods
The main geographic focus of the study was the North Sea. The first task was to analyse genetic data from a wider region of the northwest European shelf to identify reproductively isolated population units whose breeding distributions might include any part of the North Sea. Having established the geographic boundaries of such units, we then determined their biological properties from trawl survey data (growth rate and length-at-maturity, population length distributions, and biomass) and obtained data on fishery landings and discards. These data then formed the basis for parameterizing and tuning a spatial model of the interacting population dynamics of the selected units.

Genetic analysis
Fin clips from mature cod were preserved in 99% ethanol on 28 occasions during stock monitoring surveys and from commercial fisheries around the UK during spawning periods (March/April) in 1998-2007 (20 geographic localities with 11 temporal replicates; Table 1). In addition, to identify loci under selection during early life stages, cod eggs were harvested from plankton samples collected near three sampling sites for spawning fish (Table 1). Eggs visually identified as cod were manually picked out of the plankton samples at sea and preserved in 99% ethanol. Species identity was later confirmed using molecular probes (Fox et al., 2005).
DNA was extracted from preserved material by salt extraction (Aljanabi and Martinez, 1997), concentrations quantified using a Nanodrop spectrophotometer then "normalized" to 7 ng ml 21 . Approximately 40 adult individuals from each location, or 10 -21 eggs per sample, were analysed for 96 SNP loci using the KASPar system (KBiosciences, UK). Loci were chosen from previous studies (Moen et al., 2008;Hubert et al., 2010) based on their utility for discriminating geographic samples (Supplementary  Table S2). Allele frequencies and observed and unbiased expected heterozygosities were calculated for individuals with data on .85% of the loci using Geneclass2 (Piry et al., 2004), and the significance of departures from Hardy -Weinberg equilibrium tested in Genepop 4 (Raymond and Rousset, 1995;Rousset, 2008) using the Markov chain method (dememorization number ¼ 10 000, batches ¼ 100, iterations ¼ 10 000). Genotypes at all pairs of loci were also tested for linkage disequilibrium in Genepop 4 (exact test with significance levels determined by the Markov chain method with the same setting as above). To identify loci likely to be under positive selection, an independent dataset of 1200 SNP loci for five cod samples collected around the UK (Nielsen et al., 2012) was investigated using Bayesian (Foll and Gaggiotti, 2008). Six loci in the current dataset were identified as putatively under positive selection (Supplementary Table S2).
Two datasets were constructed (i) all loci and (ii) a subset of 90 "neutral" SNP loci, i.e. SNP loci not showing predicted deviations due to selection. Maximum likelihood geographic clustering of mature fish samples was estimated by the Bayesian method coded in BAPS 5 (Corander et al., 2003;ten where K is the number of clusters). Analysis of molecular variance (AMOVA) was then performed to investigate how genetic variation was partitioned within and among groups using Arlequin 3.5 (Excoffier et al., 1992). Fisher's exact tests were performed to identify differences in allele frequencies between eggs and adults collected from close geographic proximity using Genepop 4 with significance levels determined using the Markov chain method (dememorization number ¼ 10 000, batches ¼ 100, iterations ¼ 10 000). In all analyses, p-values were corrected for multiple tests (Benjamini and Hochberg, 1995).

Trawl survey data analysis
Data collected during annual quarter 1 (February -March) ICES International Bottom Trawl Surveys (IBTS) between 1980 and 2005 (ICES, 2010b) were used to construct time-series of mature and immature biomass and biological properties of cod. The IBTS survey data extend back in time to the 1965, but inconsistencies in sampling methodology and geographic extent make the pre-1980 data difficult to interpret. Post-1980, the surveys used of the Grande Overture Vertical trawl and are considered to provide a coherent record of the fish community (ICES, 2010b), despite further ongoing standardization especially during 1980-1990. Data post-2005 were not available when we carried out the analysis.
The raw data on cod from the surveys consisted of numbers of fish caught in length class intervals (usually but not always of 1 cm width), in tows of known duration and location. Distance towed and trawl geometry data (wing spread), or values imputed from statistical relationships with depth and speed (ICES, 2010b), allowed the catch rate perunittimetobenormalizedtoareaswept(km 2 ). The survey protocol provides for the replication of tows at the scale of ICES statistical rectangles (18 longitude × 0.58 latitude). Length-stratified data on the age (from otoliths) and maturity of individual fish are collected according to geographic strata comprising aggregations of 10-30 ICES rectangles, referred to as "Standard Roundfish Areas" (ICES, 2010b). However, individual length, age, and maturity records in the dataset are often not attributable to specific hauls or statistical rectangles but only to Roundfish Area.
There were six stages to processing these data to provide the annual summaries required for our analysis. For each survey year: 1. compilation of arithmetic mean catch density-at-length by statistical rectangle (normalized to swept-area); 2. construction of smoothed matrices of the probability of age given length, and maturity given length for each Roundfish Area, according to the methodology of Stari et al. (2010); 3. matrix multiplication of the outputs from (1) and (2) to produce vectors of proportion mature-at-length, and matrices of mature and immature catch density by length and age jointly, for each statistical rectangle (Stari et al., 2010)-the annual probability matrices for each Roundfish Area were applied equally to the catch densities-at-length in all their constituent statistical rectangles; 4. estimation of a length-dependent catchability correction function using the methodology of Fraser et al. (2007), and application to the outputs from (1) and (3); 5. conversion of numbers to biomass using a standard weightlength relationship (Coull et al., 1989); 6. integration of catchability-corrected data over geographic domains.
Regarding the catchability correction (stage 4 above), we followed Fraser et al., (2007) and estimate the spatially integrated annual abundances-at-age of cod in the whole survey region (numbers per age class) from the survey data and related these to the corresponding annual (1 January) abundances-at-age from the ICES stock assessments (ICES, 2010a). The ratio of survey abundance to the assessed abundance of each age class was an index of cachability-at-age. We then transformed these back into catchability-at-length [Q(L)] using a von Bertalanffy growth function for cod (see model description later) and fitted a piece-wise cubic model to the combined multiannual data: Regarding the spatial integration of survey data, mature and immature fish abundances and biomasses in statistical rectangles were estimated as the sum over all lengths, of the product of mean densities-at-length (km 22 ) and the sea surface area of cod habitat (given by water depths between 10 and 300 m). Because the surveys were carried out at around the spawning time of cod (February -March), we could reasonably assume that mature fish from reproductively distinct population units were spatially disaggregated according to their natal spawning regions. This would not necessarily be the case at other times of year. Statistical rectangles were therefore assigned to genetic population units according to the geographic groupings emerging from the SNP analysis, and the spawning-stock biomass of each genetic unit derived by summing mature biomass over all relevant rectangles. Immature biomass was not assigned to population units, since fish from different units were assumed to be intermixed throughout the year until attaining maturity. Annual recruitment to the survey domain as a whole was estimated as the population abundance of age 1 fish at the time of the survey. Domain-wide annual egg production was the sum over all length classes of the potential fecundity of mature female fish (all units combined, sex ratio assumed 50:50), with fecundity determined from a function of the gutted body weight per individual of each length class (Heath et al., 2008).

Fishery catch data
Regionally integrated estimates of the annual biomass of cod landed at ports, and discarded at sea, for the North Sea, and the west of Scotland, were available from ICES stock assessment working group reports (ICES, 2010a). Landings data spatially disaggregated by the ICES statistical rectangle in which catches were taken were obtained for the years 1997-2004 .

Spatial population dynamics model
The purpose of the population dynamics model was to simulate the responses of non-interbreeding population units of cod to changes in spatial patterns of fishing mortality. If there are no interactions between the population units then this is a trivial exercise and survival depends linearly on the mortality rate directed at each unit. However, the density-dependent processes which regulate cod population dynamics act during the early life and immature stages when fish from different units intermingle depending on interannually varying passive dispersal of eggs and larvae and active migration patterns. Since density-dependent processes presumably do not discriminate between fish from different units, the result is a competitive interaction between units and scope for non-linear responses to fishing. The key feature which the model was therefore required to represent was the time-dependent geographic distributions of demographic stages in each population unit. The model which we developed to achieve this was adapted from a single-population spatial scheme previously described by Andrews et al. (2006). Our adaptation was to represent multiple, reproductively isolated population units, with density-dependent mortality acting on the combined abundance of specific demographic stages from all units. A summary description of the model is given below, and the equations and their parameters are presented in Supplementary Material.

Model description
The model was discrete in time and space. The spatial granularity of the model corresponded to the ICES statistical rectangles, i.e. maritime spatial cells of a 18 longitude × 0.58 latitude grid, spanning 128W-108E, 498N-638N (Figure 1). At every spatial grid cell, each reproductively isolated population unit was represented by numbers in length classes for each of three life-history stages; pelagic eggs and larvae; settled immature; and mature fish. Growth in length was represented by progression through the series of classes at time-steps of 0.16 d. The length class intervals defining each population unit at each grid cell were constructed to represent equal intervals of development time (hence unequal intervals in length; Gurney et al., 2007). The length bounds of each class, given the number of classes (30 in the pelagic vector, 10 immature, 110 mature), were defined by the asymptotic length and growth velocity parameters of the von Bertalanffy growth equation estimated from analysis of length and age data collected during annual monitoring trawl surveys [Supplementary Equation (4) and Table S1].
Pelagic individuals at or exceeding a threshold size (minimum settlement size, 4 cm) were transferred to the corresponding immature demersal length class, if they were located in a spatial cell with appropriate settlement attributes and in which there was free capacity. Each settlement cell had a maximum carrying capacity for settled fish. Individuals from any population unit could settle in any suitable grid cell, so the carrying capacity referred to the combined settled fish biomass from all units. This meant that the population units were in competition with each other at settlement. Pelagic stages which failed to settle due to lack of available carrying capacity suffered an escalating mortality rate with time. This mechanism provided density-dependent control in the model. Sexual maturation was represented by the transfer of individuals from length classes of the immature stage to the equivalent length classes in the mature stage within a population unit. In the model, maturation occurred once a year, on the first development update after the beginning of October regardless of location. The proportions transferred were obtained from a Weibull distribution (Andrews et al., 2006) with year-specific parameters estimated by fitting to the proportions mature at length in the annual survey data (Supplementary Table S1).
Spawning was represented by the addition of numbers to the egg class (the first length class of the pelagic life stage) of each population unit. This only occurred at development time-steps within an assigned spawning season, in cells designated as spawning sites. The number of eggs produced at each time-step was determined by the annual potential fecundity of the mature stock in each cell, which was a function of individual body weights (Supplementary Material).
Movement was represented by relocating individuals in each length class of every grid cell, to the corresponding length classes in nearby destination cells at discrete transport time-steps. The relocation fraction for each source-destination pair was specified by transition matrixes for each development stage (pelagic, immature, mature) calculated offline by particle tracking simulations incorporating deterministic advection velocities and stochastic diffusion velocities (Gurney et al., 2001;Andrews et al., 2006). Pelagic stages were assumed to be planktonic, so for these length classes the Combination of genetics and spatial modelling advective terms in the particle tracking scheme were obtained from daily resolution hydrodynamic output of an ocean circulation model (POLOMS; Allen et al., 2001). Hydrodynamic data were available for the years 1980-2005, and transition matrices for movement updates were calculated at 7 d intervals. For mature fish, the advective terms in the particle tracking simulations represented a northerly post-spawning migration rate based in tagging data (Righton et al., 2010; Supplementary Table S1), and transition matrices were calculated at 28 d intervals. In addition, at the first movement time-step after start date of the spawning season, mature fish from each population unit, wherever they happened to be, were redistributed evenly over their respective spawning cells and remained there until the end of the spawning season. Particle tracks used to calculate transition matrices describing immature fish movements were based only to diffusion.
Domain-wide natural mortality rates were applied to all length classes of all life stages. Pelagic mortality rates were treated as yearspecific external forcing data and were derived from survey data as described below. Settled fish (immature and mature) were subject to a natural mortality rate which decreased exponentially with length up to a threshold (95 cm), above which mortality increased with length (for details, see Supplementary Material).
In addition to natural mortality, settled stages were subject to length-, time-, and space-dependent fishing mortality rates, which were external forcing factors. Deaths in the model due to fishing mortality represented catch, which was apportioned between landings and discards based on a minimum landing length (35 cm;EEC, 2006).

Model spawning and settlement locations
Each spatial cell in the model was assigned binary attributes to denote spawning and settlement locations. Assignment of cells as spawning locations (Figure 1) was based on a synthesis of historical egg surveys and trawl catches of ripe fish (Fox et al., 2008) and attributed to the population units based on results from the analysis of SNP data. Settlement cells (Figure 1) were identified based on ICES quarter 3 (July-September) trawl survey catches of age 0 cod (,12 cm; Heath et al., 2008).

Settlement carrying capacity
The carrying capacity of individual grid cells for settlement stages was a key parameter in the model. Carrying capacity was the product of an area-specific capacity parameter (ind. m 22 ) which was uniform over the model domain and constant with time, and the area (m 2 ) within each settlement cell which was both shallower than 200 m and classified as mud, sand, or coarse sediments (not rock or diamicton) from geological surveys data (British Geological Survey, 2013). This designation was based on published analyses of juvenile cod habitat preference (Borg et al., 1997;Juanes, 2007;Lough, 2010).
The area-specific capacity parameter was estimated from published data on the age 1 carrying capacity of ten subregions of the model domain (Heath et al., 2008). The origin of these data was as follows: (i) the annual proportion of regionally integrated abundance of age 1 cod resident in each of the ten subregions was estimated from the quarter 1 IBTS surveys; (ii) for each subregion, the proportions were averaged over the survey series 1983-2005; (iii) the average proportions were applied to the maximum of a smoother through the time-series of North Sea and West of Scotland numbers at age 1 from ICES stock assessment outputs; (iv) the maximum number at age 1 in each area was then inflated to numbers in August the previous year (i.e. at settlement time) using the natural mortality applied in the North Sea assessments, to provide an estimate of the carrying capacity for newly settled juveniles in each subregion. To estimate the domain-wide area-specific carrying capacity, we divided the capacity of each subregion by the sea-surface area of settlement cells within the subregion (Heath et al., 2008) and averaged over the 10 subregions (0.208, s.d. 0.093 ind. m 22 )

Mortality rates
Natural mortality of pelagic length classes was parameterized from analysis of trawl survey data, as the ratio of survey-wide egg production in a given year to survey-wide recruits (numbers age 1 in February/March) in the following year (for details, see Supplementary Material).
Forcing data on length-dependent fishing mortality rates were derived from the annual estimates of fishing mortality-at-age published by ICES assessment working groups (ICES, 2010a). For a given year, rates-at-age (F age ) were transformed to estimates of rates in model length classes (F q ), uniform over the whole model arena, using the annual age -length keys derived from trawl surveys. Mortality rates in length classes were then disaggregated into a length selectivity vector (S q ), and a term F which was the mean mortality rate over the length range 41 -79 cm (approximately equivalent to ages 2-5). Selectivity in a given length class q was given by S q = F q /F. Spatial disaggregations of model-wide fishing mortality rates F MW were carried out using trawl survey estimates of stock abundance in such a way as to preserve the global weight of fish caught in the model domain as a whole. If W N and W S are survey estimates of the harvestable biomasses of fish in two regions (N and S) with fishing mortality rates F N and F S respectively, then: If we hypothesize a fishing mortality rate F S relative to F MW , we can determine the implied rate F N required to conserve the model-wide catch: Alternatively, if F N and F S are known, then the resulting model-wide fishing mortality rate F MW is given by:

Usage of the model
We used the model in two ways, (i) to hindcast the spatial distributions of abundance and catch of the population units at 28 d intervals in response to circulation patterns and fishing mortality rates over the period 1980-2005 and (ii) to analyse the sensitivity of equilibrium states to spatial patterns of fishing mortality. For the hindcast simulation, the model was tuned to minimize the residual between modelled and survey estimates of immature and mature population unit abundances and regional landings and discards during 1980-2005 (ICES, 2010a, b). Tuning was by manually adjusting a parameter defining the relative spatial distribution of fishing mortality between two fishing subregions [F S and F N given annual forcing values of F MW in Equation (3); F N = F in ICES Area IVa (Figure 1), and F S = F in the rest of the model domain]. For all hindcast runs, the model was forced by yearspecific data on ocean circulation, maturity-at-length parameters, pelagic stage mortality, model-wide fishing mortality, and length selectivities. Sensitivity of the model to spatial patterns of fishing mortality was investigated by simulating the equilibrium state of the each population unit (landings, discards, biomass, etc.) for a matrix of 41 × 41 (i.e. 1681) combinations of fishing mortality rates in the range 0 -1.2 year 21 in each of the two fishing subregions (F N and F S ). Length selectivity by the fishery, pelagic dispersal, growth, maturation, and larval mortality rates were as pertaining in 2005. The model-wide fishing mortality rate (F MW ) for any given combination of F N and F S was then estimated by the biomass-weighted average [Equation (4)]. We then extracted F MW and the corresponding simulated data on landings, spawning biomass etc of each population unit from the matrix, along transects of constant spatial pattern of fishing [i.e. constant values of (F N /F S )].

Genetic results
No SNP loci showed consistent deviations from Hardy -Weinberg equilibrium across populations and there was no evidence of linkage disequilibrium between any pair of loci. Fisher's exact test demonstrated that early stage eggs from the northern and central North Sea plankton samples were genetically consistent with locally spawning adults and clustered with local populations. Eggs collected off northeast Scotland did not show significantly different allele frequencies, but were distinct from the local adults (but with low log-likelihood) in the cluster analysis (Table 1 and  Supplementary Table S2). Hence, there was no detectable selection occurring at any of the SNP loci between the pelagic egg and adult stages and that consequently, allele frequency differences among populations (at outlying loci) are likely to reflect actual levels of population structure rather than differential mortality between the egg and adult stage.
Maximum likelihood cluster analysis of the SNP data (see Electronic dataset) revealed three distinct spatio-genetic groups with high statistical significance (Figure 1), regardless of whether or not loci apparently under selection were included in the analysis (Table 1 and Supplementary Table S2). We refer to the groups as (i) the "Viking" unit, centred on Viking Bank in the northeastern North Sea, (ii) the "Dogger" unit, centred broadly on the Dogger Bank but found across a broad swathe of the southern and western North Sea and west of Scotland, and (iii) the "Celtic" unit, found in the western Channel, Celtic Sea, Irish Sea, and Firth of Clyde. There was no evidence that the breeding domain of the Celtic unit extended into any part of the North Sea.
The AMOVAs of both the full 96-locus and the 90-locus neutral marker datasets revealed significant differences among the three groups identified by the cluster analysis (Table 2). Genetic differentiation between the clusters was so strong as to imply effective barriers to inter-breeding. There was some evidence of finer scale differentiation of subpopulations within each cluster, but only at much lower levels of significance (Table 2). From these results, we can assert natal fidelity of cod at the geographic scales of the Viking, Dogger, and Celtic units, but not at finer spatial scales.
Bathymetry offered an effective delineation of the barrier between the Viking and Dogger units, with Viking spawning sites being confined to waters deeper than 100 m in the North Sea. Boundaries between the Celtic and Dogger units off the north of Ireland and in the English Channel were not obviously related to bathymetry.

Selection of population units to include in the model
Based on the results from the genetic sampling, we confined the spatial population model, which was focused on the North Sea, to representing the Viking and Dogger units but not the Celtic unit. This was because the spawning ranges of both the Viking and Dogger units included parts of the North Sea, whereas the range of the Celtic unit did not. In addition, evidence from tagging and otolith microchemistry indicated a low likelihood of intermingling between Celtic unit juveniles and adults and those from the other units (Wright et al., 2006a, b;Righton et al., 2010). Nevertheless, our genetic data indicated that the geographic extent of the model needed to extend beyond the confines of the North Sea to include waters west of the UK to include the full breeding and distributional range of both the Viking and Dogger units.

Biological properties of the Viking and Dogger population units
Data from the IBTS surveys of the North Sea and west of Scotland regions (ICES, 2010b) encompassed almost the entire geographic ranges of the Viking and Dogger units as identified from the SNP data, except for the eastern English Channel. Catchability-corrected abundances of mature fish at age and length were integrated over statistical rectangles deeper than 100 m to represent the Viking unit and over rectangles shallower than 100 m to represent the Dogger unit (Figure 1). Rectangles representing the Viking unit were mainly contained within a single Standard Roundfish Area (Area 1; ICES, 2010b) which also encompassed a few rectangles shallower than 100 m especially around the Shetland Isles. Hence, the underlying age and maturity data upon which the Viking fish were characterized may have been contaminated by Dogger fish, but it was not possible to determine the extent of this, since the exact locations of individual age and maturity samples within Roundfish Areas were often not recorded.
There was no statistically significant difference in length-at-age between mature fish assigned to the Dogger and Viking units (Figure 2). At age 6, the median length of mature Dogger fish slightly decreased over time (linear trend 20.064 cm year 21 ), whereas Viking fish increased (+0.546 cm year 21 ), but neither of these trends were significant (p .. 0.05). In contrast, there was a highly significant (p , 0.01) trend in maturation characteristics over the 1980-2005 survey period. Length at 50% mature declined from 65-70 cm in 1980 to 45 -50 cm in 2005 ( Figure 2). The data suggested that Viking fish matured on average 2-3 cm larger than Dogger fish in the early years of the series, but the difference has diminished over time and was small relative to the overall timetrend in maturation rates. Based on these results, we assumed the same fixed von Bertalanffy growth parameters for both population units but resolved a linear trend in maturity-at-length in each unit (Supplementary Table S1).

Trends in Viking and Dogger population units and model hindcast
The trawl survey data showed that biomass of mature fish in waters shallower than 100 m during February/March (Dogger unit) halved over the period 1980-2005, but there was no significant trend (p . 0.05) in mature biomass in waters deeper than 100 m deep (Viking unit; Figure 3). Also, the observed length distributions of mature fish in the Dogger unit contained a lower proportion of large fish (.50 cm) compared with those for the Viking unit ( Figure 4) implying greater mortality rates. The tuned model provided a credible hindcast of all the observed time-dependent features of the regional stock and fishery, including the dynamics of the Dogger and Viking population units (Figure 3), and the population unit length distributions (Figure 4, Supplementary Figures S1 -S3). The optimal model fit was obtained when fishing mortality rates in the southern North Sea and west of Scotland (ICES Areas IVb, IVc, Via, and VII) were on average 1.2 times higher than in the northern North Sea (ICES Area Iva; Figure 3).

Model scenario analyses
The global maximum sustainable yield from the system (given 2005 environmental and maturity conditions) was achieved by the extreme scenario of harvesting cod only in the northern North Sea (north of 57830 ′ N; F N ≈ 0.8, F S = 0; Figure 5). In this regime, the Dogger unit was harvested mainly during its migration phase while the majority of its spawning sites were protected from Figure 2. Quarter 1 (January -March) survey estimates of lengths at age 2 and 6 of mature fish (upper) and length at maturity (lower) for fish in the Dogger and Viking population units by year of sampling. Ages 2 and 6 are chosen to illustrate the growth pattern of fish. Left panels, fish assigned to the Dogger population unit, right to the Viking unit. Thick lines show the median (50th centile) of length-at-age (upper) or length at 50% mature (lower); dashed lines show 25th and 75th centiles. t-test statistics showed no significant trends in lengths at age 2 and 6 or differences between the units. Trends in length at 50% maturity were highly significant (p , 0.01), with a marginally significant difference between units (p ≈ 0.05).
fishing-a classic spill-over harvesting scenario (Gerber et al.,2003). However, under these conditions, the population diversity was lost because the Viking unit was driven to extinction. With more realistic spatial distributions of fishing (0.667 ≤ F N /F S ≤ 1.5), maximum equilibrium landings resulted from model-wide fishing mortality rates between 0.3 and 0.4 year 21 (Figures 5 and 6). However, at these model-wide rates of fishing, the Viking unit was predicted to become extinct if fishing was spatially biased towards the northern North Sea by a factor of ≥1.5. On the other hand, the Viking unit persisted at 30 -40% of unfished biomass if Model-wide fishing mortality rate (year 21 ) was an external driving series; the ratio of rates in the northern North Sea (ICES Area IVa) and everywhere else in the model domain (ICES Areas IVb, IVc, Via, and VII) was a tuning parameter. The pelagic stage mortality rate index was also an external driving series derived from survey data. Landings (kt year 21 ), discards (kt year 21 ), immature biomass (kt), and recruits (×10 6 individuals at age 1 in quarter 1) refer to the model domain as a whole. Spawning-stock biomass (kt, mature males and females combined) refers to the Dogger and Viking population units. The cumulative proportion of annual landings shows the relative spatial distribution between regions-the area between each line represents the simulated proportion of total landings taken from ICES Fishing Areas (Figure 1); point values represent the proportions from each area compiled from national statistics (only available from 1997 to 2004).

Combination of genetics and spatial modelling
fishing mortality was concentrated in southern North Sea by a factor of 1.5. These results were only marginally affected by the substitution of the contrasting 1980 maturity-at-length schedules in the model instead of those applicable in 2005 ( Figure 6).
To diagnose the mechanism governing the response of the Viking unit to spatial distributions of fishing, the matrices of 41 × 41 combinations of fishing mortality were recreated, first with only the Viking unit present in the model, and second with only the Dogger unit present. The results ( Figure 6) showed that the abundance of the Dogger unit was only weakly sensitive to the presence of the Viking unit. However, the Viking unit was highly sensitive to the presence of the Dogger unit, with the unfished spawningstock biomass of the Viking unit being approximately four times higher with the Dogger unit absent. In addition, spatial variations in the fishing mortality rate had almost no effect on the Viking unit when the Dogger unit was absent from the model (Figure 6). We conclude that the sensitivity of the Viking unit to spatial distributions of fishing was not simply a direct effect of fishing mortality, but included indirect effects from changes in the abundance of the Dogger unit acting through the density-dependent mortality enhancement of settling juveniles. The Dogger unit, being the more abundant at equilibrium, was relatively unaffected by changes in the Viking unit.

Distribution of population units relative to stock management regions
The population units of cod identified by our analysis of SNP data do not map comfortably onto the existing cod management regions. The genetic evidence, supported by tagging data (Righton et al., 2010), clearly show that cod spawning in the Firth of Clyde, Irish Sea, Celtic Sea, and western Channel form a single breeding unit (Celtic unit). Hence, the scientific basis for separate assessment and management of the Irish Sea (ICES Division VIIa) and the Celtic Sea/Western Channel (Divisions VIIe -k) seems weak. Similarly, the rationale for treating the west of Scotland (Division VIa) as a distinct stock for management purposes is not sustained by the genetic data. Mature cod caught off the west of Scotland showed a highly significant genetic affinity with cod found across a broad swathe of the North Sea (Dogger unit), whereas fish in the Firth of Clyde, which is also technically included in the west of Scotland management region, are clearly part of the Celtic unit. Tagging data also show extremely limited exchange of cod between the Clyde and the open shelf west of Scotland (Wright et al., 2006b).  Based on these results, there would seem to be a strong biological reason for combining assessments of cod in the North Sea and west of Scotland (excluding the Firth of Clyde). However, from a management perspective, there are clearly strong reasons for setting TACs at smaller spatial scales which resemble the underlying geopolitical and economic structures in the fishing industry. This is one of the issues with which spatially resolved population dynamics models, such as we describe here, can help-by simulating the effects of subregional spatial patterns of harvesting rate. The other key issue with which spatial models can help is the situation which occurs in the North Sea, where multiple genetic units are contained within the same management region (Viking and Dogger units).
Here, the challenge is to conserve the population diversity by managing the overall and/or spatial distribution of harvesting to protect the least productive genetic units.

Differentiation of population units
The key to specifying our model structure was unequivocal evidence of geographic barriers to gene flow between population units. Our SNP data clearly differentiated the Celtic, Dogger, and Viking population units, indicating a high degree of natal fidelity at these geographic scales. However, the emergent biological differences between fish from the Viking and Dogger units were slight according to the data available to us, with the main difference being in length at maturity-Viking fish matured at a larger size (≏5 cm) than Dogger fish.
Gene flow sufficient to maintain genetic homogeneity within each of these units does not necessarily preclude biological differentiation at smaller spatial scales as a result of phenotypic plasticity. If we regard the environmentally determined component of the phenotype (as opposed to the genetically determined) as being the response of traits to environmental conditions (e.g. growth rate in relation to temperature), then environmental gradients can lead to spatial patterns in, e.g. length at age, or maturation depending on migration patterns within a given genetically homogeneous population unit. Hence, other investigators have found subregional variations in length-at-maturity of cod within the domain of the Dogger unit (Wright et al., 2011;Neuheimer and Grønkjaer, 2012). In our model, we attenuated the free exchange of fish between spawning locations in southeastern and northwestern parts of the Dogger unit domain to reflect migration patterns previously demonstrated from tagging data (Righton et al., 2010). However, in the absence of any information at the level of the population units on responses of growth, maturation and fecundity to environmental conditions, we assumed that each unit was homogeneous with respect to biological properties, although this is apparently not strictly true in reality.
We found no temporal trends in length-at-age of mature cod from either the Viking or Dogger units in our analysis, in apparent contradiction to Neuheimer and Grønkjaer (2012) who found a North Sea-wide increasing trend in length-at-age, since the mid-1990s related to the temperature history of successive cohorts. However, Neuheimer and Grønkjaer compiled the data for mature and immature fish combined, whereas we consider only mature fish disaggregated by population units. Hence, it is rather difficult to make a direct comparison between the two results. Nevertheless, the major temporal trend in the data is the declining length at maturity. We can speculate that this may be an evolutionary response to intense harvesting.
Apart from length at maturity, the other major differentiation between the Viking and Dogger units is in the dispersal patterns of eggs and larvae from the spawning areas. According to the particle tracking simulations, a consistently lower proportion of eggs shed in the Viking area reach suitable settlement habitat than from the Dogger area (Heath et al., 2008). This is the main factor which leads to lower productivity of the Viking unit.

Interaction between population units and the regulation of population dynamics
The fact that less productive population units in a mixed assemblage are more prone to depletion by harvesting is not especially surprising. However, our model results show unexpectedly complex responses of population diversity to variations in spatial patterns of fishing due to interactions between units arising from densitydependence. Competition between population units was by far the dominant factor compared with simply the direct effect of a redistribution of exogenous fishing mortality. Density-dependent controls on the population dynamics of cod must act during the late-pelagic and early demersal stage, because regional recruitment is set by age 1 (ICES, 2010a), and field observations show strong habitat dependence of mortality rates during this phase (Juanes, 2007;Lough, 2010). We can assume this mortality to be indiscriminate with respect to fish from different population units, creating the competition which we represented in our model. Even if the details of the process in the model are not strictly true to nature, and we do not have precise information on this, the effect is to capture the essence of the multipopulation control mechanism with an emergent density-dependent mortality which is patchy and transient in space and time depending on the dispersal of pelagic stages. This contrasts with previous modelling studies where density-dependent regulation has been implemented at the level of individual populations through parameterized relationships between unit recruitment and unit spawning-stock biomass (Cope and Punt, 2011;Ying et al., 2011), which take no account of the combined resource requirements of spatially overlapping units.
Competition between population units in our model leads to marked divergence between the traditional fisheries management controls required to ensure sustainable yield, and those required to also ensure long-term persistence of population diversity. Both our data analysis and the hindcast model fit indicate that fishing mortality rates on cod have historically been lower in the northern North Sea compared with the shallow southern waters, which has fortuitously protected the Viking unit from depletion. However, this situation is threatened by a shift of mixedspecies demersal fishing effort into deeper water over recent years (Greenstreet et al., 1999(Greenstreet et al., , 2009, which is unregulated since it does not imply any change in overall effort within the North Sea.

Implications of our model for stock management
The problems caused by mismatch between population units and fisheries management areas are well recognized, but until now, integrating population genetics into fisheries management has been difficult (Reiss et al., 2009;La Valley and Feeney, 2013)a s current harvest management regions are long-established and so enshrined in national and international legislation that any revision would be extremely complicated to negotiate. However, our combined genetic and modelling approach could provide the basis for a generic solution, by providing a means of retaining the existing regional management structures, but estimating a precautionary regional harvesting rate which would ensure that all population units remain at safe levels of abundance, given prognoses of the likely spatial distribution of fishing. For example, our model indicates that achieving regional scale maximum equilibrium landings of cod-referred to as maximum sustainable yield which is a key target for fisheries management (Caddy and Mahon, 1995;Hilborn, 2007)-is commensurate with a regional spawning-stock biomass at 20% of the unfished state. With hydrodynamic, pelagic survival and maturity schedules applicable in 2005 this corresponds to a regional fishing mortality rate of 0.35 year 21 (Figure 6) which is within the range of fishing mortalities considered by ICES to be consistent with maximum sustainable yield (0.16-0.42 year 21 ). However, if the prognosis was that the spatial distribution of fishing was likely to be 1.5 times higher in the northern North Sea than the south, then the regional harvesting rate required to ensure that the Viking unit does not fall below 20% of unfished biomass would need to be reduced to ≏0.2 year 21 .Our continuing research will investigate the sensitivity of these prognoses to expected environmental conditions (pelagic stage mortality and ocean circulation) and updating the hindcast runs beyond 2005.

Conclusions
Using new SNP data, we show that cod stocks in the North Sea and waters west of the UK include three genetically distinct population units, thereby resolving prior uncertainty based on microsatellite data. In the North Sea, fish from two of these units (Dogger and Viking) spawn at different locations, but intermingle at other times. Our data analysis shows differing trends of abundance over time for these units, and our model indicates that this can arise from a combination of spatial patterns of harvesting and competition between the population units. Management restraints on cod fishing aim to deliver maximum sustainable yields at the scale of the whole North Sea, but our model shows that this poses a quantifiable threat to maintenance of the underlying population diversity if fishing on one unit is allowed to become dominant. Redefining the boundaries of management regions is probably politically unattainable, but our model can advise on precautionary rates of North Sea-wide fishing required to protect extant population diversity.

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