Population assignment and local adaptation along an isolation‐by‐distance gradient in Pacific cod (Gadus macrocephalus)

Abstract The discernment of populations as management units is a fundamental prerequisite for sustainable exploitation of species. A lack of clear stock boundaries complicates not only the identification of spatial management units, but also the assessment of mixed fisheries by population assignment and mixed stock analysis. Many marine species, such as Pacific cod, are characterized by isolation by distance, showing significant differentiation but no clear stock boundaries. Here, we used restriction‐site‐associated DNA (RAD) sequencing to investigate population structure and assess power to genetically assign Pacific cod to putative populations of origin. Samples were collected across the species range in the eastern Pacific Ocean, from the Salish Sea to the Aleutian Islands. A total of 6,425 putative biallelic single nucleotide polymorphisms were identified from 276 individuals. We found a strong isolation‐by‐distance signal along coastlines that mirrored previous microsatellite results and pronounced genetic differentiation between coastal samples and those from the inland waters of the Salish Sea, with no evidence for hybridization between these two populations. Individual assignment success based on two methods was high overall (≥84%) but decreased from south to north. Assignment to geographic location of origin also was successful, with average distance between capture and assignment location of 220 km. Outlier analyses identified more loci potentially under selection along the coast than between Salish Sea and coastal samples, suggesting more diverse adaptation to latitudinal environmental factors than inshore vs. offshore environments. Our results confirm previous observations of sharp genetic differentiation of the Salish Sea population and isolation by distance along the coast, but also highlight the feasibility of using modern genomic techniques to inform stock boundaries and fisheries management in a low F ST marine species.


| INTRODUC TI ON
The identification of essentially self-recruiting populations and the estimation of dispersal rates have been the primary aim of population genetics and management for many years (Palsboll, Berube, & Allendorf, 2007). In terrestrial and freshwater systems, the expectation of sharp boundaries between genetically differentiated populations is often confirmed by genetic data, allowing unequivocal delimitation of management units (Palsboll et al., 2007). In marine species, on the other hand, with their large population sizes and comparatively high rates of gene flow, population boundaries are often subtle (Waples, 1998), and populations may be structured by limited dispersal distances rather than clear barriers to gene flow (Wright, Bishop, Matthee, & von der Heyden, 2015). The resulting isolation-by-distance (IBD) patterns (Rousset, 1997;Slatkin, 1993) are of limited use to fisheries managers, because statistical areas cannot be clearly defined (Spies, Spencer, & Punt, 2015), and arbitrary delineation of stock boundaries is unacceptable to stakeholders. Genetic markers under selection that are often detected by next-generation sequencing approaches may provide better resolution than more traditional markers (Hauser & Carvalho, 2008) and may also allow for the assessment of dispersal and seasonal migration based on individual assignment to population of origin, a technique that is rarely used in low F ST marine species.
Pacific cod (Gadus macrocephalus) in the northeast Pacific Ocean is such a species, where restricted dispersal distances result in a clear IBD signal (Cunningham, Canino, Spies, & Hauser, 2009), thus complicating the identification of clear management units. Only cod in the Salish Sea have been identified as a clearly differentiated and spatially segregated population (Canino, Spies, Cunningham, Hauser, & Grant, 2010;Cunningham et al., 2009), which together with the declining abundance of cod in that region led to the listing of Pacific cod in 2010 as a NOAA Species of Concern in the Salish Sea. Along the Aleutian Islands, Gulf of Alaska, British Columbia, and Washington State, there was a tight IBD pattern without clear stock boundaries (Cunningham et al., 2009). A follow-up study found an almost identical IBD pattern within Alaska but also identified subtle population boundaries (Spies, 2012), in part leading to the separation of the eastern Bering Sea and the Aleutian Islands management areas, which had been managed as a single unit until 2015. However, cod is still managed in three very large geographic management areas: eastern Bering Sea (EBS), Aleutian Islands (AI), and Gulf of Alaska (GOA), despite strong evidence for short generational dispersal distances <100 km (Cunningham et al., 2009) and consequent potential for population structure and local adaptation on much smaller geographic scales.
Pacific cod is an excellent system to test the power of nextgeneration sequencing approaches to detect management units along an IBD gradient, not only because of the species' limited distribution along the continental shelf, but also because of its significant economic, ecological, and conservation interest. Pacific cod is a demersal species distributed along the continental margin of the North Pacific to approximately 400-m water depth (Cohen, Inada, Iwamoto, & Scialabba, 1990). There are two major phylogeographic groups in the northeastern and northwestern Pacific, which appear to represent recolonizations from different refugia with a boundary in the western Aleutians (Canino et al., 2010;Cunningham et al., 2009). In Alaska, Pacific cod move into shallower waters between 100-m and 200-m depth to aggregate for spawning in early spring (Neidetcher, Hurst, Ciannelli, & Logerwell, 2014). Spawning is widely distributed, with some concentrations associated with island topography . Despite seasonal migrations, adult Pacific cod are noted for their spawning site fidelity and homing (Rand, Munro, Neidetcher, & Nichol, 2014;Shi, Gunderson, Munro, & Urban, 2007;Shimada & Kimura, 1994), which supports short lifetime dispersal estimates from genetic data (Cunningham et al., 2009). Fertilized Pacific cod eggs are negatively buoyant and adhesive until hatching (Fredlin, 1985), providing an additional mechanism to restrict dispersal distances. This distribution of spawning sites along a relatively narrow band along the shelf of the North American continent is expected to cause a linear IBD pattern (Slatkin, 1993), an expectation confirmed in two independent microsatellite studies (Cunningham et al., 2009;Spies, 2012).

Recent advances in next-generation sequencing technologies
have allowed for a more thorough survey of genetic variability and differentiation across the genome (Allendorf, Hohenlohe, & Luikart, 2010). One such technique, restriction-site-associated DNA (RAD) sequencing, has revolutionized population genetics in nonmodel organisms by allowing discovery and genotyping of thousands of single nucleotide polymorphisms (SNPs) in multiple individuals at relatively low cost (Baird et al., 2008;Miller, Dunham, Amores, Cresko, & Johnson, 2007). This large number of genetic markers not only affords much higher power in statistical tests of population differentiation and population assignment, but may also provide evidence for local adaptation from outlier loci that show particularly high differentiation among populations (Hauser & Seeb, 2008). Such outlier loci may be particularly valuable in discriminating populations where traditional markers show a simple IBD pattern without clear population boundaries (e.g., Atlantic cod in Norway, Hauser & Carvalho, 2008). Furthermore, outlier loci may provide valuable information on local adaptation of Pacific cod in the Salish Sea near the southern edge of the distribution, possibly explaining its continued decline despite protection for the past 30 years.
In this study, we used RAD sequencing of Pacific cod to (i) increase resolution of stock boundaries along the coastal IBD gradient in the northeast Pacific Ocean from Washington State to the Aleutian Islands, (ii) evaluate the performance of methods to accurately assign individuals geographically to putative source populations, and (iii) address the potential for selective differentiation, both along the latitudinal gradient from Alaska to Washington and between Salish Sea cod and nearby coastal populations.

| Sample collection and RAD sequencing
Samples were collected from spawning and prespawning aggregations of Pacific cod at seven locations across the northeastern Pacific Ocean (Table 1, Figure 1). An eighth sample of fish from the Salish Sea region was collected in 2012 in the United States and 2013 in Canadian waters (Figure 1 inset)

| Data analysis
Raw data were converted to genotypes using Stacks v1.21 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013) according to the methods of Gruenthal et al. (2014), with minor modification. Briefly, catalogs created in the cstacks subprogram were generated from the five most data-rich individuals from each sample. Flags (m = 3, M = 2, N = 4, n = 3, max_locus_stacks = 3) associated with increasing the number of loci, while reducing SNP and allele calling error rates, were set according to Mastretta-Yanes et al. (2015). A genotype file containing putative polymorphic SNPs present in ≥80% of fish per sample was filtered to include one SNP per RAD tag (flag: write_ran-dom_SNP) to minimize physical linkage, and SNPs with more than two alleles were discarded. Final filtering removed the last nucleotide position on the tag (base pair 94) and loci with minor allele fre-  (Rousset, 2008) using the program default parameters. Population pairwise F ST and associated p-values were also estimated using GENEPOP v4.2. The genetically effective size (N e ) of each population was estimated using NeEstimator v2.01 (Do et al., 2014) under the random mating model using the linkage disequilibrium method (Waples & Do, 2008), with an MAF cutoff of 0.05 (R. Waples, NOAA, personal communication (Meirmans & Hedrick, 2011) was estimated using diveRsity (Keenan, McGinnity, Cross, Crozier, & Prodöhl, 2013) within R for both microsatellite and RAD data, and pairwise linearized genetic distance (G ′′ ST /(1-G ′′ ST )) was plotted against shortest geographic distance along the continental shelf estimated from Google Earth.

| Outlier testing and alignment to the Atlantic cod genome
Outlier tests to identify putative candidate loci under selection were performed using two different methods. First, a Bayesian framework with a differentiation-based method was employed in BayeScan v2.1 (Foll & Gaggiotti, 2008). Testing was conducted with 20,000 iterations, prior odds of the neutral model of 100, and a type II false discovery rate of α = 0.05. Second, OutFLANK (Whitlock & Lotterhos, 2015) was used to identify candidate loci under selection.
OutFLANK estimates the distribution of F ST values at neutral, or nearly neutral, loci by fitting the empirical data to a chi-square distribution after trimming excessively high and low F ST values, as these loci may be under diversifying or balancing selection. The distribution is then compared to the empirical data, and outliers are identified as those outside the expected distribution. In this analysis, five percent of F ST values were trimmed from each tail of the distribution before estimating the null chi-square distribution. Samples from the coastal distribution in the northeastern Pacific (i.e., excluding the Salish Sea) were first analyzed to examine coastal patterns of selection, and then, a subset of samples (SS12/13, WC05, and HS04) were analyzed to examine differentiation associated with Salish Sea samples versus coastal samples at similar latitude.
To assess whether the outlier loci colocalized with regions identified as under selection in the congeneric Atlantic cod (G. morhua) by Hemmer-Hansen, Therkildsen, Meldrup, and Nielsen (2014) and identify putative candidate genes involved in selection, Pacific cod sequences identified as outliers were aligned to the Atlantic cod genome (gadMor1 genome assembly, http://www.ensembl. org/Gadus_morhua) using BLASTN (Altschul, Gish, Miller, Myers, & Lipman, 1990). Scaffolds or contigs that aligned uniquely to each locus with an e-value less than 10 −30 were retained. Genome scaffolds were then assigned to linkage groups by alignment to an Atlantic cod linkage map (Borza, Higgins, Simpson, & Bowman, 2010) using Bowtie (Langmead, Trapnell, Pop, & Salzberg, 2009).

| Population assignment tests
To assess the power of the data set to correctly assign individuals to their putative population of origin, we carried out three different  (Rannala & Mountain, 1997) and a leaveone-out procedure to avoid high-grading bias (Anderson, 2010). All 6,425 loci were included in analyses, and results were compared to population assignment from microsatellite data (Cunningham et al., 2009).
Second, we tested our ability to assign individuals to their sample of origin with subsets of loci using Assigner (Gosselin, Anderson, & Bradbury, 2016). Assigner uses a training data set to identify highly discriminatory loci (based on F ST ), followed by a leave-one-out method on an independent test data set to test the assignment (Anderson, 2010). Half the individuals in each sample were used as a training data set, and the remaining individuals were used to assign individuals to populations. Data sets with the 10, 50, 100, 200, 500, 1,000, and 6,425 (all) loci with the highest F ST were used, with ten iterations per data set.
Third, we used an assignment test that assigns individuals to a geographic location within an IBD pattern rather than to a pop- Each sample was assigned to a geographic location three times using different seed values to assess whether model parameters facilitated the convergence of results. The median value of the three assignments was used as the final assignment position.
Correlations between assigned geographic locations of origin and capture locations were determined using a randomization test.
Individuals were randomly assigned to capture locations, and the distance between the actual capture location and assigned geographic location of origin was compared between the randomized data set and the empirical data, with associated p-values calculated in R (R Core Team, 2017).

| Individual ancestry
Individual ancestry was estimated for all southeast samples (SS12/13, JDF12, WC05, and HS04) to further resolve population structure and gene flow between the Salish Sea and coastal samples. Individual ancestry was estimated using Structure (Pritchard, Stephens, & Donnelly, 2000). We used 50 000 and 100 000 burn-in and MCMC iterations, respectively, and K-values ranging from one to six run three times each. The optimal K-value was determined using the method of Evanno, Regnaut, and Goudet (2005) implemented by Structure Harvester (Earl & Vonholdt, 2012). Individuals with a Q < 0.9 of any one population were considered of mixed ancestry (Vähä & Primmer, 2006). Similar analyses in coastal samples only reproduced the isolation-by-distance pattern shown by other analyses and are not shown.

| RE SULTS
A total of 6,425 putative biallelic SNP loci and 297 individuals were retained after filtering. Pairwise coefficients of relatedness (r) ranged from 0.00 to 0.99, with only 19 comparisons showing evidence for a genetic relationship other than unrelated (r = 0.32 ± 0.251). Sixteen pairwise comparisons were at the half-sibling level (r = 0.22 ± .065) and three at full-sibling status (r = 0.84 ± 0.229; Table S1). Two putative full-sibling pairs in SS12/13 generated r values in excess of 0.95, which was consistent with the same fish being sampled twice.
The individual in each pair with fewer data was removed from fur-  Table S2. Global H E averaged 0.195 (SD = 0.0023) across populations (see also

| Outlier testing and genome alignment
Similar patterns of outlier loci were identified by both BayeScan and OutFLANK (Table 3). Across the coastal samples, 76 and 77 loci showed significant (FDR q < 0.05) evidence for selection by BayeScan ( Figure 4a) and OutFLANK ( Figure S2), respectively, and 65 loci were identified by both analyses (Table 3). BayeScan identified 48 of these loci as having "decisive" evidence for selection (log(PO)>2). Between coastal and Salish Sea samples (WC05 and HS04 vs. SS12/13), nine and 36 outlier loci were identified by BayeScan and OutFLANK, respectively, and seven loci were identified as outliers by both methods (Table 3). Four of the nine loci had decisive evidence for selection according to BayeScan  (Table S3).
A large number of loci aligned uniquely to the Atlantic cod genome.
Of the total 6,425 loci identified in this study, 5,007 aligned uniquely to the gadMor1 genome assembly and 866 of these aligned to coding regions. Forty-four of 65 outlier loci identified by both BayeScan and OutFLANK in coastal samples aligned uniquely to the Atlantic cod genome and 14 aligned uniquely to coding regions (Tables 3 and   S3). Five of seven outliers identified by both analyses in the comparison between the Salish Sea to coastal samples aligned to unique scaffolds and one aligned to a coding region (Tables 3 and S3). In total, four Atlantic cod scaffolds (GeneScaffold_1739, GeneScaffold_1755, GeneScaffold_1981, GeneScaffold_2531) had more than one outlier locus align to them in coastal comparisons, while one scaffold (GeneScaffold_2140) had more than one locus align in the Salish Sea versus coastal comparison. In total, 23 outlier loci were assigned to Atlantic cod linkage groups (Table S3). Four linkage groups (10, 14, 17, and 20) contained only one outlier, while linkage groups 1, 2, 8, 10, and 16 contained two, two, three, and four outlier loci, respectively (Table   S3). Lastly, locus 15623 was the only outlier to colocalize to one of the

| Population assignment
Overall, GeneClass2 assigned 88% of individuals back to the sample from which they were captured ( Among northern samples, assignments appeared to be pulled toward a center point between KOD03 and UP03. Likewise, samples from the south appeared to be drawn to a center point between HS04 and WC05 ( Figure 6).

| Individual ancestry
Structure results, with an optimal K value of two, suggested that JDF12 individuals were largely fish of coastal origin (Figure 7). All WC05 and HS04 fish clearly clustered together as a genetically homogenous group, while all SS13 and SS12 individuals clustered together as a distinct group in the Salish Sea. Of the JDF12 fish, 14 (87.5%) individuals had similar ancestries as coastal fish, and two JDF12 individuals appeared to cluster with the Salish Sea sample.
Only three fish (2.6% of total; one JDF12 and two SS12) may have had mixed ancestries between the coast and Salish Sea (Q ≥ 0.10; Vähä & Primmer, 2006).

| D ISCUSS I ON
Our results confirmed the geographic patterns in population structure using other genetic markers. The IBD gradient previously reported Numbers on the diagonal represent the total number of outliers identified by each test for each data set-for the two BayeScan analyses, the second number is the number of loci with decisive evidence for selection (log(PO)>2). Numbers below the diagonal represent the number of outliers in common between the two analyses. Numbers within parentheses represent how many outliers aligned uniquely to the Atlantic cod genome (gadMor1) and number of outliers that aligned uniquely to a coding region, respectively.
for North American coastal Pacific cod using mitochondrial DNA sequences (Canino et al., 2010) and microsatellites (Cunningham et al., 2009;Spies, 2012) was here reproduced with SNPs. In addition, the relatively high genetic differentiation between the Salish Sea and coastal Pacific cod reported for microsatellites (Cunningham et al., 2009) and mtDNA (Canino et al., 2010) was supported by this study, although the level of differentiation was lower than in these previous studies. Despite this similarity in overall genetic patterns between

Number of loci Assignment success (Proportion)
the three marker systems, only RAD data provided the power to accurately assign individuals to population and location of origin.

| Patterns of population differentiation
The slope of the IBD relationship was remarkably consistent with that found in two previous microsatellite studies (Cunningham et al., 2009;Spies, 2012), once estimates of differentiation were adjusted for differing levels of variability between marker systems (Meirmans & Hedrick, 2011 (Knutsen et al., 2010;Waples, 1998). However, the confirmation of the slope by a completely different marker system further supports the genetic pattern itself and the underlying basis of limited dispersal distances, despite a lengthy larval period and highly migratory adults. It also supports the increasingly acknowledged pattern of subtle population structure in species with no obvious barriers to dispersal, such as marine fishes (Hauser & Carvalho, 2008).
The tight isolation-by-distance pattern along the coast seems to be a contradiction to the DAPC that very clearly shows three clusters, with two clusters along the coast (Figure 3). Sampling was not evenly spaced (Figure 1) as the commercial fishery is concentrated west of Prince William Sound (Barbeaux et al., 2017), and samples from southeast Alaska were not available. The clustering in the DAPC is therefore most likely to be a result of the geographic distribution of samples, rather than true genetic differentiation among isolated populations. The high correlation coefficient of the IBD relationship (r = 0.9, Figure 2) supports such an interpretation, and the DAPC clearly indicates IBD within Alaska (Figure 3), which was also found in an independent study (Spies, 2012). Furthermore, the previous microsatellite study failed to identify genetic barriers in that region (Cunningham et al., 2009). Nevertheless, finer scale sampling would be required to determine the exact mechanisms of genetic differentiation, which could either be caused by limited dispersal within a continuous population or a stepping-stone pattern of migration along a series of isolated subpopulations (Meirmans, 2012).
The significant genetic differentiation documented previously between the Salish Sea and coastal Pacific cod (Canino et al., 2010;Cunningham et al., 2009) was also observed in this study. However, G ′′ ST between the two groups was much higher with microsatellites than with RAD tags (Figure 2). This difference was not due to selection at microsatellites. Eight of the 11 microsatellites supported the split between the two populations at roughly equal levels (Cunningham et al., 2009). Similarly, evidence for selective differentiation was found at only a few SNP loci and would result in the opposite pattern. However, microsatellite R ST (an F ST measure accounting for differences in allele sizes) was almost twice as high as F ST , suggesting mutation rather than just genetic drift as an important contributor to the observed differentiation (Canino et al., 2010).
Mutations accumulate much faster in microsatellites than SNPs (Weber & Wong, 1993) and in isolated populations would arguably also have a larger effect in multiallelic microsatellites than in biallelic SNPs, where mutations do not lead to new alleles but only marginally change allele frequencies.
However, despite the lower levels of differentiation observed from RAD-seq data in this study, ancestry estimates from Structure showed clearer separation between coastal and Salish Sea samples than for microsatellites (Canino et al., 2010), likely because of higher power conferred by the large number of loci. Indeed, there was little F I G U R E 6 Assignment locations for each sample using the leave-one-out technique and SCAT to assign to a location of origin based on genetic background. Plus signs represent assignment locations, while triangles are capture locations WC HS PWS KOD UP AD F I G U R E 7 Ancestry results for samples from the southeastern portion of the range. Ancestry proportions were determined using Structure with 50,000 and 100,000 burn-in and MCMC iterations and the previously determined optimal K-value of 2

SS12
SS13 JDF12 WC05 HS04 Ancestry evidence for hybridization between the two populations (Figure 3), despite spatial overlap in the Strait of Juan de Fuca. Pacific cod collected in the Strait of Juan de Fuca did not constitute a "population," but were instead made up mostly of individuals from the Washington coast and two Salish Sea cod, as indicated by GeneClass2, Structure, and DAPC results. Other than these two individuals, there was a sharp genetic break between the two populations near Victoria, BC, Canada. Similar genetic differentiation between Salish Sea and coastal populations has been found previously in hake (Iwamoto et al., 2015), brown rockfish (Buonaccorsi, Kimbrell, Lynn, & Vetter, 2005), and copper rockfish (Buonaccorsi, Kimbrell, Lynn, & Vetter, 2002).
This high differentiation with apparent lack of hybridization may seem surprising on such a small geographic scale and in an oceanographic system as dynamic as the Salish Sea. The genetic break coincides with the Victoria Sill, a shallow submarine feature that divides the Strait of Juan de Fuca from the interior Salish Sea (Ebbesmeyer & Barnes, 1980). This sill rises to depths of 55 m to 100 m and causes considerable vertical mixing during strong tidal currents in the Strait of Juan de Fuca (Ebbesmeyer & Barnes, 1980;Khangaonkar, Long, & Xu, 2017;Soontiens & Allen, 2017). This mixing creates spatial discontinuities in sediment load (Johannessen, Masson, & Macdonald, 2006), salinity (Masson & Cummins, 2000), nutrient load (Peña, Masson, & Callendar, 2016), and primary production (Masson & Peña, 2009). Sills are therefore seen as a mechanism limiting larval dispersal, and indeed, similar sills in Norwegian fjords seem to represent a barrier for the export of pelagic eggs of Atlantic cod, thus explaining the small-scale differentiation found there (Knutsen et al., 2007). Adults, however, seem to be able to transverse such barriers, as evidenced by the two Salish Sea individuals found in the western Strait of Juan de Fuca. It therefore seems possible, that in addition to such extrinsic mechanisms maintaining isolation, intrinsic genomic mechanisms may be at play. Indeed, chromosome inversions in Atlantic cod in Scandinavia seem to enhance adaptive differentiation between populations (Barth et al., 2017), including stationary and migratory ecotypes , and populations inhabiting fjords and outer coasts . Such inversions are likely to cause genome regions of high differentiation among populations ("islands of divergence," Sodeland et al., 2016).
Given the weak evidence for adaptive divergence in Salish Sea cod found in our study, there is, however, currently no evidence for such inversions in Pacific cod.
A potential secondary split was found within the SS12/13 sample (Figure 3 bottom). Although more similar to one another than to any other sample (Figure 3 top), Pacific cod collected in 2012 in U.S. waters may be different from those collected in 2013 in the northern Georgia Basin in Canada. Consistent with the north-south break, another shallow sill is located north of the San Juan Islands at Boundary Pass (Masson, 2006), and the Fraser River freshwater plume and deltaic plain extends from Vancouver, BC, Canada, across the Basin toward Vancouver Island. One or both of these physical features may help isolate Pacific cod in the deeper northern Georgia Basin from those to the south. Indeed, freshwater and sediment barriers have been implicated in driving population divergence and endemism in marine taxa elsewhere. The Amazon, Yangtze, and Mississippi River outflows, for example, limit cross-barrier dispersal in corals, limpets (Cellana toreuma), seahorses (Hippocampus spp.), various reef fishes, red snapper (Lutjanus campechanus), and other organisms (Boehm et al., 2013;Dong et al., 2012;Potts, 1983). Yet, it must be noted that temporal structuring, which was not examined in this study, cannot be ruled out.
The seemingly limited gene flow between Pacific cod in the Salish Sea and coastal populations may have important management implications. Pacific cod abundance in the Salish Sea has declined dramatically over recent decades (Johannessen & McCarter, 2010;Palsson, Hoeman, Bargmann, & Day, 1996). In 2000, the National Marine Fisheries Service published a status review of Pacific cod in Puget Sound and determined that their listing under the Endangered Species Act was not appropriate at the time (Gustafson et al., 2000).
This determination was partially based on the lack of genetic information and an unclear northern boundary delineation that they hypothesized likely extends as far north as the border between British Columbia and southeastern Alaska. Microsatellite (Cunningham et al., 2009) and mtDNA (Canino et al., 2010) Gustafson et al., 2000). Given these results, another status review of Salish Sea cod may be warranted.

| Evidence for selective differentiation
There was strong evidence from outlier locus testing for selective differentiation across the range. Sixty-five loci were putatively categorized as being under selection across the sampled coastal range in both outlier tests (BayeScan and OutFLANK), and seven were identified among HS04, WC05, and SS12/13. The identification of these loci by two different approaches supports the notion that they are truly under selection, although some of the observed patterns may still be the result of IBD (Lotterhos & Whitlock, 2014;Meirmans, 2012). The adaptive significance of these outlier loci is currently unknown, but candidate outlier genes identified by alignment to the Atlantic cod genome provide suggestive evidence that reproductive biology may play an important role in differentiation among coastal samples. First, locus 53780 had the largest overall F ST value of any locus (F ST = 0.71) and was identified as an outlier by both BayeScan and OutFLANK. This locus aligned to the coding region of a novel protein sequence that is orthologous to genes involved in reproduction and sexual development: steroidogenic factor 1 (SF1) and zona pellucida glycoprotein 3 (ZP3). Second, locus 8253 had the second largest overall F ST value (F ST = 0.46), was identified as an outlier in both coastal analyses, and aligned to an intron of the folliculogenesis-specific bHLH transcription factor, which may regulate the expression of the reproductively important zona pellucida genes (Liang, Soyal, & Dean, 1997). Zona pellucida has been previously identified to undergo rapid selection in other species (Aagaard, Yi, MacCoss, & Swanson, 2006;Heras, McClintock, Sunagawa, & Aguilar, 2015), and both changes to coding regions as well as differences in expression patterns may be important to diversification in Pacific cod.
Fewer outlier loci were identified in the comparison between Salish Sea and coastal cod, a surprising result, given the sharp genetic break between populations and the dearth of identifiable hybrids. However, outlier tests have notoriously low power with few sampled populations (Foll & Gaggiotti, 2008), and the limited number of populations may have led to many false negatives. In addition, the overall high F ST between coastal Pacific cod and Salish Sea cod may have reduced our power to detect outliers in that comparison. Nevertheless the strong genetic differentiation of Puget Sound cod (Cunningham et al., 2009) and their long isolation from coastal populations (Canino et al., 2010) provide the opportunity for adaptive genetic differentiation and local adaptation to warmer conditions in the southern limit of the distribution.
Bottom temperatures during the larval period (February, March, 7-8°C; Moore et al., 2008) routinely exceed optimal temperatures for cod larval rearing (4°C; Alderdice & Forrester, 1971), suggesting potential effects of long-term heat stress on growing larvae in Puget Sound. Similar patterns have been shown in Atlantic cod, where recruitment success of several stocks is tightly linked to sea surface temperatures during the larval period (Brander, 2010).
Further evidence for local adaptation in Pacific cod stems from steep latitudinal clines of increased growth and mortality rates of Puget Sound cod (Karp, 1982), although the underlying genetic basis of such differences is currently unknown. Local adaptation has been demonstrated in Atlantic cod at very small geographic scales from both common garden experiments (Hutchings et al., 2007) and molecular studies (Case, Hutchinson, Hauser, Van Oosterhout, & Carvalho, 2005). Local adaptation increases the biocomplexity of the species as a whole, which may help prevent extreme abundance fluctuations in a portfolio effect (Schindler et al., 2010). Fisheries management approaches that incorporate evolutionary and plastic biocomplexity have been suggested for Atlantic cod, for example, in Olsen et al. (2008). Moreover, populations at the southern edge of the species' distribution may be old and harbor genetic variation important for the survival and evolution of a species (i.e., rear edge effect; Hampe & Petit, 2005).
These stabilizing characteristics may become increasingly important as climate change progresses (Hauser & Carvalho, 2008).

| Population assignment
Poor population assignment rates are often expected for species exhibiting IBD (Pritchard et al., 2000) depending on the level of differentiation and the spatial scale of sampling. Nevertheless, RAD sequencing increased our ability to re-assign individuals to their population of origin relative to microsatellite loci used in previous studies. Excluding JDF12, median reassignment success using all RAD tags was ≥85% for each method (Table 4 and Supporting information), while rates for microsatellites using GeneClass2 averaged only about 32% (L. Hauser, unpublished data). The higher assignment success of RAD sequencing is likely due to both the greater number of loci in this study and the inclusion of outlier loci with high levels of differentiation.
Despite the recent and ongoing decline in sequencing costs, RAD sequencing is too costly for routine application on large-scale fisheries samples. However, simulations based on loci identified in a training data set and applied in a leave-one-out fashion on a holdout data set (Anderson, 2010) suggested that 100 to 200 loci may be sufficient to approach a similar assignment success as for all 6,425 loci. Such a set of 100 loci could be easily scored using reduced representation libraries (e.g., GT SEQ ; Campbell, Harmon, & Narum, 2015), and opens potential applications to source verify catch or fish products (Nielsen et al., 2012), estimate contributions to mixed stock fisheries (Ruzzante, Taggart, Lang, & Cook, 2000), directly estimate dispersal rates (Bradbury et al., 2011), and identify skiped spawning cod remaining on the feeding grounds (Skjaeraasen et al., 2012), among other questions. Whether ~80% assignment success is sufficiently high will depend on the specific question, but it does represent a significant increase in power over microsatellites, and, to our knowledge, is among the highest in marine fishes to date.
Increasing assignment capabilities in the marine environment may have implications for the future of Pacific cod fisheries management. Currently, Pacific cod established fishery management units may not reflect suspected stock structure (Spies, 2012). Although the Bering Sea/Aleutian Islands management area was separated into two units beginning in 2015, there is evidence that stocks may be subdivided over smaller spatial scales. For example, cod fisheries in the U.S. Exclusive Economic Zone (EEZ; 3-200 nm offshore) in the Gulf of Alaska Management Unit fall under federal authority, while the State of Alaska manages the fishery within state territorial (0-3 nm) waters. Typically, the Alaska Department of Fish and Game (ADF&G) creates parallel fishing seasons, allowing vessels to fish for cod at the same time as federally controlled fisheries. The significant differentiation between KOD03 and PWS12 found in this study suggests potential isolation of the Prince William Sound population, which is managed by ADF&G and may require independent assessment to annual federal stock assessment reports (Thompson & Palsson, 2013). When management units and associated harvest quotas are not based on realized stock structure, the spatial allocation of fishing may not be proportionate to stock biomass and could result in localized overfishing, particularly if two or more stocks are managed as one . Moreover, disproportionate exploitation of certain stocks may have the potential to destabilize the adaptive resilience of the stock complex as a whole (Hilborn, Quinn, Schindler, & Rogers, 2003).
The success of population assignment tests is not only limited by the power of genetic markers, but also by the availability of baseline samples for assignment. In marine fishes with relatively poorly known and potentially spread-out spawning locations  as well as IBD patterns, such population assignment can become very problematic. However, by interpolating allele frequencies in a very defined IBD relationship, SCAT (Wasser et al., 2004) provided some very promising identification of location of origin. To the best of our knowledge, this study is the first to use this method in a low F ST marine fish species. Misassignment distances (i.e., the distance between capture location and putative location of birth) in this study (median = 220 km) were greater than generational dispersal from IBD patterns (<100 km; Cunningham et al., 2009). Numerous factors may contribute to this discrepancy. First, there is undoubtedly a statistical error in the estimation of location of origin, although further work is required to quantify that error. A sign of this statistical error may be the tendency of assigning cod to locations at the center of the distribution, both in the Alaskan peninsula/Gulf of Alaska and in the Washington/BC area ( Figure 6). In an IBD pattern, misassigments to adjacent locations are expected (and indeed observed with traditional population assignment, Figures S2-S8), but this tendency toward the center of the distribution may have more to do with geographic location of samples and edge effects. Second, SCAT assigns individuals to a population of origin based on their genetic background. As such, progeny of immigrants may have an inflated migration distance that is influenced by their pedigree. Third, although every effort was made to collect spawning fish, samples may have unintentionally included some skip spawners from external populations, thus inflating misassignment and potentially biasing baseline allele frequencies. Skipped spawning has been documented in a wide range of marine species and likely occurs also in Pacific cod (Rideout & Tomkiewicz, 2011).
Despite the challenges of fine-scale geographic assignment, the results observed here can serve as a starting point for using genetics to investigate movement patterns of Pacific cod and other marine species. Tagging studies (Nichol, Kotwicki, & Zimmermann, 2013;Shimada & Kimura, 1994) have indicated seasonal movements "related to, but not limited to, spawning and feeding" (Rand et al., 2014), but coherent spatiotemporal patterns have not emerged.
Although no range-wide studies have been performed, Shimada and Kimura (1994) estimated that 15%-17% of Pacific cod from the eastern Bering Sea may migrate to the Gulf of Alaska to spawn, and a significant proportion of the stock may undertake large migrations (>100 km; Rand et al., 2014;Shimada & Kimura, 1994). Currently, tagging studies are the standard method of studying movement patterns of Pacific cod, but such studies are expensive in terms of ship time, labor intensive (requires capturing the same fish twice), potentially affected by tagging mortality and atypical behavior of tagged fish, and often limited to few fish being recovered. In contrast, with genetic tools, all individuals in a population are effectively "tagged," and so all captured individuals can be assigned to a geographic location of origin, resulting in considerably larger sample sizes and more effective coverage of the population. Results from this study suggest that assignment methods may provide an informed approach toward resolving seasonal movement; fish captured during the nonspawning season (April-December) could be assigned geographically using SNP data to putative spawning source populations collected from January to March. This approach may also allow distinction of stationary (i.e., nonmigratory) and migratory ecotypes that are suggested in Pacific cod (Rand et al., 2014) and documented in Atlantic cod (Karlsen et al., 2013).

| CON CLUS IONS
Our study provided unprecedented resolution of the stock structure of Pacific cod along the North American coast. Although previously identified isolation-by-distance patterns were confirmed, the higher power of RAD sequencing provided highly significant differentiation among almost all samples. Population structure relevant to management was evident on large (Washington/British Columbia vs. Alaska) as well a small geographic scales (e.g., Prince William Sound vs Kodiak Island), thus providing important information for the definition of management units. Potentially, more importantly, we established the utility of assignment to both population and location of origin and applied that latter approach for the first time in a marine species. SNP-based genetic stock identification has the potential to provide deeper resolution of population structure and understanding of seasonal migrations in a species where spatial dynamics are not well characterized . Knowledge of seasonal movement and mixing may provide support for the development and implementation of management practices, such as determination of stock-specific total allowable catches designed to preserve spatiotemporal stock structure and avoid local depletion.
More generally, the application of next-generation sequencing techniques in combination with sophisticated statistical approaches may provide a significant advance in the management of exploited marine species such as Pacific cod.

ACK N OWLED G EM ENTS
This work was funded in part by a grant from Washington Sea Grant, University of Washington, pursuant to National Oceanic and Atmospheric Administration award no. NA10OAR4170057.

CO N FLI C T O F I NTE R E S T
None declared.