Genetic evidence of a northward range expansion in the eastern Bering Sea stock of Pacific cod

Abstract Poleward species range shifts have been predicted to result from climate change, and many observations have confirmed such movement. Poleward shifts may represent a homogeneous shift in distribution, seasonal northward movement of specific populations, or colonization processes at the poleward edge of the distribution. The ecosystem of the Bering Sea has been changing along with the climate, moving from an arctic to a subarctic system. Several fish species have been observed farther north than previously reported and in increasing abundances. We examined one of these fish species, Pacific cod, in the northern Bering Sea (NBS) to assess whether they migrated from another stock in the eastern Bering Sea (EBS), Gulf of Alaska, or Aleutian Islands, or whether they represent a separate population. Genetic analyses using 3,599 single nucleotide polymorphism markers indicated that nonspawning cod collected in August 2017 in the NBS were similar to spawning stocks of cod in the EBS. This result suggests escalating northward movement of the large EBS stock during summer months. Whether the cod observed in the NBS migrate south during winter to spawn or remain in the NBS as a sink population is unknown.


| INTRODUC TI ON
Several recent studies have confirmed the prediction that increases in water temperature driven by climate change can cause range shifts of marine species toward higher latitudes and contraction at lower latitudes (Booth, Bond, & Macreadie, 2011;Pinsky, Worm, Fogarty, Sarmiento, & Levin, 2013;Sunday, Bates, & Dulvy, 2012).
As productive, commercially fished species move northward, understanding the effect of changes in climate on species range shifts into higher latitudes and movement into marginal habitat will be particularly challenging. Even understanding whether habitat is marginal is difficult as conditions shift under climate change. Nevertheless, understanding changes in distribution is important for their conservation through dynamic fisheries management, to engage in decision-making that is based on the best estimates of stock dynamics and spatial distribution.
Populations in marginal habitats are expected to be less abundant than those in core habitats, and survival and reproductive rates may be lower (Kawecki, 2008;MacCall, 1990;Sexton et al., 2016). Newly established marginal populations are less well adapted, with low reproduction and survival, and must be maintained by migration from the core (Kawecki, 2008). Marginal populations may then adapt to local conditions but are still prone to local extinctions and can act as demographic sinks (Kawecki, (Alabia et al., 2018;Stabeno et al., 2012). Since 2014, the eastern Bering Sea (EBS) has exhibited anomalously warm temperatures (Alabia et al., 2018). As a result, the "cold pool," water below 2°C that remains along the EBS shelf during the summer following sea ice retreat, has been restricted to the northern parts of the EBS. Most recently, in 2018, the National Marine Fisheries Service (NMFS) survey of the EBS observed the smallest cold pool in the survey history , with only 1% of the total area of the EBS shelf bottom at less than 2°C (Siddon & Zador, 2018;Stabeno et al., 2018). The cold pool is an important factor affecting species distributions; for example, walleye pollock (Gadus chalcogrammus), Pacific cod (Gadus macrocephalus), and most flatfishes avoid it (Hollowed, Planque, & Loeng, 2013;Sigler et al., 2016;Stevenson & Lauth, 2019). Reductions in marine mammal and benthic prey populations have been observed in the northern Bering Sea (NBS) concurrent with increases in pelagic fish, as the NBS is shifting from arctic to subarctic conditions (Grebmeier et al., 2006). tons in 2018 (Siddon & Zador, 2018). In contrast, the first of the two most recent and comprehensive surveys, conducted in 2010, found very few Pacific cod anywhere north of 60°N latitude, estimating 29,091 tons for the entire northern region, or approximately 3.3% of the biomass for the EBS shelf for that year (Stevenson & Lauth, 2019). Data from previous surveys indicate that the historical range of Pacific cod did not include the northern portion of the Bering Sea; Pacific cod accounted for less than 1% of the total gadid biomass and were encountered "only in trace amounts" in the 1976 and 1979 surveys of Norton Sound and Chirikov Basin (Sample & Wolotira, 1985;Wolotira, Sample, & Morin, 1977). The EBS is considered part of the core habitat for the species, and longterm average   800,000 t, ten times higher than the Aleutian Islands and three times larger than the GOA (Barbeaux et al., 2017;Schmidt et, 2019).
However, the cod biomass measured by the NMFS EBS shelf survey declined 37% between 2016 and 2017, representing the largest decline since the survey began in 1982. Strikingly, the 2017 summer cod biomass in the NBS was equal to 83% of the reduction in biomass in the EBS in the same year and by 2018, summer cod biomass was higher in the NBS than the EBS (Figure 2, Table S2; Thompson, 2018).
The life history of Pacific cod is characterized by seasonal movement between summer feeding and winter spawning locations (Shimada & Kimura, 1994). This seasonal movement coupled with genetic evidence for population structure among spawning groups indicates that natal homing and spawning site fidelity are part of Pacific cod life history, and thus, collections of spawning samples are representative of a population while samples collected during summer migrations may be mixed. Spawning areas have been identified along the Bering Sea shelf, Aleutian Islands, and Gulf of Alaska ( Figure 1; Neidetcher, Hurst, Ciannelli, & Logerwell, 2014). It is doubtful that cod have historically spawned in the NBS, since seasonal sea ice typically covers the northern part of the Bering Sea shelf during January-May, even in warm years, and historical surveys observed few cod in the NBS (Sample & Wolotira, 1985;Stabeno et al., 2017;Wolotira et al., 1977).
Shifts in the spatial distribution of commercial species have the potential to disrupt fisheries with challenging implications for sustainable management (Karp et al., 2019;Pinsky et al., 2018). In the Bering Sea, Pacific cod fisheries are rationalized and managed through a complex catch share program involving a diverse suite of fishing sectors (Fina, 2011;Ono et al., 2017;Stram & Evans, 2009).
Of particular importance to our study, the northern region of the Bering Sea shelf is closed to commercial trawling, limiting access of one sector of the fleet when the population shifts north. The trawl fishery primarily operates during the fall and winter when fish have returned to the south to spawn. If climate change disrupts or retards this southward migration, fishing opportunities would diminish the catch potential for this portion of the fleet.  Cadrin et al., 2013). However, as in many other marine species, the mechanisms of northward shifts and the status of new northern populations are currently unknown.
The goal of this study was to examine the mechanisms and genetic effects of shifts in the abundance of Pacific cod in the northern portion of its range within the context of our current understanding of range expansions into marginal habitat. This work built upon baseline information (Drinan et al., 2018), which showed that spawning populations of Pacific cod were sufficiently discriminated by SNP loci for successful assignment to population of origin. We used 3,599 single nucleotide polymorphism (SNP) markers obtained through restriction-site associated DNA sequencing (RADseq) to compare a sample of Pacific cod taken from the NBS during a 2017 research survey to spawning fish from three locations in the EBS and three additional populations throughout their range in Alaska ( Figure 1). We tested several hypotheses to explain the origin of increasing abundances of Pacific cod in the NBS, specifically whether the NBS cod stock represents: 1. an isolated preexisting population, 2. a new population established via a founder event, 3. a range shift from core habitat elsewhere in the range of Pacific cod, 4. a sink population that is maintained by immigration from core habitat elsewhere in the range of Pacific cod, or 5. an extension of summer feeding migrations. Hypothesis 1 would be consistent with significant genetic differentiation among NBS samples and other known Pacific cod populations, hypothesis 2 would suggest a reduction in relative genetic diversity in NBS cod, while hypotheses 3, 4, and 5 would indicate genetic similarity with the stock of origin. Distinguishing among hypotheses 3, 4, and 5 will require additional research, but the management implications among them are significant. Therefore, results are discussed in the context of other relevant biological information that may provide some further insight and context for the current and future status of the population.  (Table 1). Previously published sequences (Drinan et al., 2018) were combined with new data for a total of 360 individuals (Table 1).

| RAD sequencing
Fin clips were preserved in 100% nondenatured ethanol. For the NBS collection, Pervenets, Pribilof, and Unimak, RAD libraries were prepared as described in Drinan et al. (2018) and sequenced in 100 bp single-end reads on three lanes on a HiSeq4000 (Illumina,

Inc.) at the University of Oregon Genomics and Cell Characterization
Core Facility (GC3F).  (1) Note: Samples marked with an asterisk * were sequenced previously (Drinan et al., 2018). N represents the number of samples that passed data quality checks, and numbers in parentheses represent the number removed.  A single SNP at each locus with the highest MAF was retained for further analysis.
Conformation of genotype proportions to Hardy-Weinberg expectations (HWE) was tested by chi-squared tests in pegas (Paradis, 2010) for each collection separately. Loci were removed if they were not present in at least six of the seven populations, significantly out of HWE following correction for multiple testing by the false discovery rate (FDR) method (Benjamini & Hochberg, 1995), or out of HWE in two or more populations.
BayeScan uses differences in allele frequencies between populations to identify loci under selection via a q-value, a stringent FDR analog of the p-value. OutFLANK infers the distribution of neutral markers from a trimmed distribution of F ST values and uses that distribution to detect outliers. BayeScan was run for a total of 100,000 iterations, a burnin of 50,000, and thinning interval 10, for a total sample size of 5,000. The prior odds for the neutral model was set to 100 and α = .05 (FDR). OutFLANK was run with left and right trim fractions 0.05, minimum heterozygosity required to include a locus 0.1, and a q-threshold 0.05. Linkage group and position were taken from the stacks output file sumstats.tsv and aligned with all four GadMor2 annotation files (https ://osf.io/4qsdw/ ) (Tørresen et al., 2017) using the Bioconductor package GenomicRanges (Lawrence et al., 2013).
F IS and Fisher's exact tests for differentiation were calculated using genepop with dememorization of 10,000; 1,000 batches; and 5,000 iterations per batch (Rousset, 2008). The number of alleles per population and expected heterozygosity (H e ) were calculated using adegenet (Jombart & Ahmed, 2011). Effective population size was estimated for each population using NeEstimator v2.01 with the random mating model, linkage disequilibrium method, parametric 95% confidence intervals, and minimum allele frequency cutoff of 0.05 (Do et al., 2014). Hierfstat (Goudet, 2005) was used to calculate rarefied allelic counts and pairwise F ST (Weir & Cockerham, 1984).

| RE SULTS
Data processed on the Illumina HiSeq 4000 averaged 3.8 × 10 8 total sequences per library, 5.4% ambiguous barcode drops, 0.6% lowquality read drops, and 88.5% retained read rate. Data processed on an Illumina HiSeq 2000 presented lower average scores: 1.9 × 10 8 total sequences per library, 13.0% ambiguous barcode drops, 8.0% low-quality read drops, and 53.1% retained read rate. A total of 6,133 loci were retained following initial STACKS v1.44 data processing. Five individuals were removed during the STACKS pipeline due to insufficient data (Table S1).
There were 3,731 loci remaining after selection for loci present in at least six out of seven populations and in HWE by population.
Sixty-six loci were removed due to linkage disequilibrium, and the largest linked cluster consisted of three loci. Seven individuals were removed with more than 30% missing data (Table S1).  Figure 3).
F IS was generally low and positive, indicating a heterozygote deficiency, which can be due to inbreeding or allelic dropout (  F ST ranged from 0.0002 to 0.0142 across all comparisons (Table 3).  (Table S3). OutFLANK identified 218 outlier loci in the full dataset and four in the Bering Sea dataset, while BayeScan identified 90 outlier loci in the full dataset and one in the Bering Sea dataset. All but one of the candidate loci for selection identified by BayeScan in the full dataset was also identified with OutFLANK. These 89 loci mapped to a limited number of linkage groups: LG2,6,8,9,11,14,and 15,suggesting that some linkage groups may contain more genes under selection than others ( SNPs identified as outliers, indicating diversifying selection rather than balancing. BayeScan results are shown visually for the full and Bering Sea datasets, with log 10 (q-value) plotted versus F ST , and a vertical line depicting the threshold level for significance ( Figure   S1).
Effective population size estimates were generally large and ranged from 1,691 in the Adak population to over 5,000 in samples from the Bering Sea, although results are likely imprecise as some values could not be estimated ( The estimated posterior density for the RUBIAS mixture proportion of the NBS individuals was over 98% (0.957, 1.000) EBS when individuals from Pervenets, Pribilof, and Unimak were combined into a single EBS reference population (Figure 5a, Table 4a). When Pervenets, Pribilof, and Unimak were treated as separate reference populations, the mixture proportions were 30% Pervenets, 19% Pribilof, and 50% Unimak (Figure 5b), with <1% for populations outside of the EBS (Table 4b).
Comparison of the Z-scores to the normal distribution provided information on whether the NBS collection came from a population not represented by reference populations. The Kolmogorov-Smirnov test rejected the null hypothesis that the log-likelihood Z-score came from a normal distribution when all EBS samples were excluded from the reference but not otherwise ( Table 5). Removal of Pervenets, Pribilof, or Unimak from the reference dataset did not significantly change the Z-score, even though each had a high probability of assignment from the NBS (Table 5). Simulations indicated that mixture proportions were highly accurate for all reference groups when EBS collections were combined (Figure 6a). When considered separately, there was a definitive loss of accuracy in the mixture proportions of Pervenets, Pribilof, and Unimak populations (Figure 6b), which was likely associated with low levels of differentiation among those populations.

| D ISCUSS I ON
We found evidence for large-scale summer redistribution of Pacific cod from their historical range in the southern EBS to ~1,000 km north during anomalously warm conditions between 2010 and F I G U R E 5 Posterior density of the probability of assignment of the northern Bering Sea sample to reference spawning populations, labeled as "Population" in the legend. In panel ( 2017. Movement appeared to be most likely from spawning populations in the EBS core habitat into the NBS, which has historically been marginal habitat as is evidenced by low encounter rates for surveys conducted in the region before 2017 (Bakkala et al., 1992;Bakkala, Traynor, Teshima, Shimada, & Yamaguchi, 1985;Goddard & Zimmermann, 1993;Sample & Wolotira, 1985;Stevenson & Lauth, 2012Walters et al., 1988;Wolotira et al., 1977;Wolotira, Sample, Noel, & Iten, 1993). Individual cod collected from the NBS in August of 2017 assigned to EBS spawning populations with high levels of certainty (>98%) and not to the Aleutian Islands or the GOA ( Figure 5; Table 4). This indicated support for hypotheses 3, 4, and 5; the movement into the NBS represented either a large-scale redistribution (hypothesis 3), a sink population (hypothesis 4) from EBS cod spawning populations, or an extension of the northward feeding migration of the EBS stock (hypothesis 5). It is difficult to distinguish among hypotheses 3, 4, and 5 without further analyses, but a feeding migration is the simplest and most consistent with known cod life history. Pacific cod sampled in the NBS had above average condition in 2018, and anecdotal observations indicated their stomachs were full of snow crab, Chionoecetes opilio (Siddon & Zador, 2018). The northern Bering Sea may represent a new location for summer feeding migrations of EBS cod in years with a diminished cold pool. If cod migrate northward to feed during the summer, they may either undergo long migrations back to their spawning location of origin (hypothesis 5) or perhaps overwinter in the NBS if conditions allow (hypotheses 3 or 4).
It is possible that the NBS collection, or some portion thereof, could have originated from Russian waters or another unsampled location, if that population were genetically similar to the reference populations from the EBS. Pacific cod spawn along the Bering Sea shelf, and fishing takes place along the shelf on both sides of the Russia-U.S. maritime boundary during spawning season. It is unlikely that cod in Russian water adjacent to the Pervenets sampling location are strongly genetically distinct from other EBS cod, but this has not been tested. Nonetheless, Pervenets, Pribilof, and Unimak represent the largest and most geographically proximate known spawning areas in the EBS and are therefore the most likely sources of the NBS collection (Neidetcher et al., 2014).
Results indicate against hypothesis 1, a founder event, and hypothesis 2, that large numbers of cod in the NBS originated from a pre-existing population in the region. Genetic results did not provide evidence that cod observed in the NBS were genetically distinct from populations of cod in the EBS (Tables 2 and 3). The level of allelic richness was not lower in the NBS collection than any of the EBS populations, which would be expected in the case of a founder effect. Rather, slightly higher allelic richness in the NBS may indicate that it represents a mixture of other stocks, a reasonable assumption considering that stocks of cod likely mix during the summer feeding season (Table 2). Further, surveys indicated very low abundance of Pacific cod in the NBS prior to 2010.
The length-frequency distribution of cod collected in the NBS was also similar to that of the EBS population, although slightly larger (Stevenson & Lauth, 2019), which could be due to better feeding conditions or if cod preferentially undertake more distant summer migrations as they grow. Oceanographic evidence provides a mechanism for northward transport of juveniles. It is unknown why younger fish would move northward, but it has been shown that young fish follow older fish to learn migration routes (Dodson, 1988). The smallest fish collected from spawning populations was 46 cm from Pervenets and could have been a maturing 3 or 4 years old, as 50% maturity is estimated at 4.8 years in EBS cod (Figure 3, Thompson, 2018). The smallest fish sampled in the Northern Bering Sea were 33 cm and 2 years old, likely too small to swim 600-1,000 km if they originated along the Bering Sea shelf (Figure 3). There is evidence to suggest that larval fish can position themselves in the water column to take advantage  (Thompson, 2018). Hyperaggregation occurred in Atlantic cod in the Bonavista Corridor off Newfoundland, Canada, where the density of cod from 1990 to 1993 was fourfold that of the 1980s even though abundance had declined fivefold (Rose & Kulka, 1999). While hyperaggregation has been documented in Atlantic cod and other fish species (Erisman et al., 2011;Neuenhoff et al., 2018), it has not yet been evaluated in Pacific cod. Length-weight residuals, used as a proxy for fish condition, have declined in Pacific cod in the EBS since 2003 but increased in the NBS in 2018 relative to the EBS population (Siddon & Zador, 2018), indicating that NBS summer feeding is favorable for cod. Although the spawning stock was large in 2018, estimated annual recruitment has been below average since 2013 (Thompson, 2018). Anomalously low recruitment may be an indication of a disruption in typical spawning behavior or a reduced spawning stock, although estimates of total spawning biomass have appeared relatively stable since 1990s ( Figure 1.21, Thompson, 2018). Together, this evidence points to reduced recruitment or reduced spawning activity in the EBS, but high-quality foraging in the NBS.
Significant stock status changes, likely coupled to environmental change, have occurred recently in other parts of the range of Pacific cod, as well. In the GOA, Pacific cod encountered a different fate than Bering Sea cod during the recent warm period. During 2016-2017, warmer temperatures led to increased metabolic demands in Pacific cod, which may have exceeded energetic consumption and resulted in lower body condition (Barbeaux et al., 2017). Indeed, this elevated mortality was likely the cause of a 58% decline in Pacific cod biomass in the GOA estimated from NMFS biomass trawl survey data in 2017 relative to 2016, the lowest estimate from the time series  of standardized surveys by more than half (Barbeaux et al., 2017). In addition, at the southern end of their range in Puget Sound, WA, abundance declined as predicted in a range contraction typical for ectotherms in general (Sunday et al., 2012). Catches of cod in Puget Sound averaged approximately 2 million pounds from 1958 to 1967 (Alderdice & Forrester, 1971). Today, however, few cod currently reside in Puget Sound, and abundance estimates from WDFW trawl surveys are considered unreliable due to the infrequency of cod encounters within most Puget Sound sub-basins (Pacunski, R., WDFW, pers. comm.).
Several tangential results are worth mentioning that provide potentially useful management information. Significant F ST values among EBS shelf populations were present even though the magnitude of F ST was too small to distinguish among those populations using assignment testing. Fisher's exact test has been shown to have low statistical power for SNPs; therefore, the significant levels determined using pairwise F ST are considered more reliable (Ryman et al., 2006). The significant pairwise F ST estimates among EBS populations (Pervenets, Pribilof, and Unimak) warrant further research on the level of genetic differentiation among those populations, perhaps using a more powerful dataset, such as whole-genome sequencing.
A lack of differentiation among the NBS collection and Bering Sea spawning populations provides support that those fish represented a mixture of Bering Sea populations. Overall, effective population size estimates are likely inaccurate, and larger sample sizes are likely required to obtain more precise estimates (Marandel et al., 2019).
Even so, effective population sizes support independent population size estimates that indicate the Bering Sea stock is larger than either the Gulf of Alaska or the Aleutian Islands. The data used in this study consisted of primarily selectively neutral loci, with roughly 89 loci in the full dataset identified as candidates for diversifying selection.
Assignment testing is not commonly applied to marine populations, which are typified by low levels of F ST . Incorporating a combination of neutral SNPs and SNPs under diversifying selection is a useful technique to improve assignment success, as in Drinan et al. (2018).
The dataset used here included many of the same loci in that study, including one SNP found within the zona pellucida subunit 3 coding sequence.
The results of this study are directly relevant to the management of this valuable resource underscoring the growing importance of genetics as a tool for sustainable management in a changing climate. Our study supports the hypothesis that climate change will extend the range for many subarctic species including Pacific cod (Cheung, Reygondeau, & Frölicher, 2016). However, our study does not support the hypothesis that climate-induced shifts in suitable habitat will also result in increased catch potential. If the low recruitment trend is not abated, the stock and the future catch potential will decline despite the expansion of suitable habitat. The findings supported continued single stock management in the eastern Bering Sea, although observations of cod along border of the United States and Russia suggests that transboundary management agreements may be necessary for the future (Pinsky et al., 2018). Results also emphasize the need for ecosystem-based fisheries management as a holistic approach that encompasses all interactions within the ecosystem rather than a single-species approach. In addition to genetics, tagging studies and continued research surveys are essential for understanding how a shifting climate will influence the extent of postspawning migrations in existing populations (e.g. from the EBS) 1