Strontium isotopes delineate fine-scale natal origins and migration histories of Pacific salmon.

Highly migratory organisms present major challenges to conservation efforts. This is especially true for exploited anadromous fish species, which exhibit long-range dispersals from natal sites, complex population structures, and extensive mixing of distinct populations during exploitation. By tracing the migratory histories of individual Chinook salmon caught in fisheries using strontium isotopes, we determined the relative production of natal habitats at fine spatial scales and different life histories. Although strontium isotopes have been widely used in provenance research, we present a new robust framework to simultaneously assess natal sources and migrations of individuals within fishery harvests through time. Our results pave the way for investigating how fine-scale habitat production and life histories of salmon respond to perturbations-providing crucial insights for conservation.

INTRODUCTION

Effective conservation of highly mobile species (for example, birds and anadromous fish) is very difficult (), primarily because of limitations in our ability to trace individuals throughout their lives (). Migratory organisms traverse large geographic areas, which encompass not only breeding and nonbreeding habitats but also migratory corridors, bottlenecks, and refugia—all of which can be important to survival at some point during an individual’s life cycle (). Tracing lifelong habitat use and habitat productivity at the level of individuals and populations will enhance our ability to protect critical habitats of highly mobile species. This is especially important given that environmental changes, including habitat loss and climate change, affect the population dynamics of mobile species in unpredictable ways (, ).

Conservation of Pacific salmon (Oncorhynchus spp.) is especially problematic because of their high mobility and complex population structure. The Bristol Bay region of Alaska produces some of the largest remaining wild salmon runs in the world. Three relevant perturbations to this region include climate change (), industrial development (), and commercial fisheries (). These generally come in the form of changes to oceanic and freshwater habitat conditions or overfishing of distinct populations. Precise natal homing of adults coupled with heterogeneity in spawning environments gives rise to locally adapted populations. This biodiversity imparts resiliency to environmental change, lending temporal stability to their overall productivity and dependent human systems (for example, fisheries) (–). Broad-scale (watershed-specific) and fine-scale (stream-specific) levels of population structure, plus variation in life history, contribute to this phenomenon (). In ecology, this phenomenon has been likened to the portfolio effect (, ), whereby a diversification of contributing populations imparts stability to overall productivity. If Bristol Bay sockeye salmon (Oncorhynchus nerka) were instead one homogeneous population (compared to hundreds), interannual variability in runs would be 2.2 times greater and would lead to complete fishery closures 10 times the current frequency (one closure per 100 years) (). Because of the ecosystem services and economic value provided by salmon populations, as well as their high exploitation rates, there has been considerable effort invested into tracking the relative productivity of populations.

One accurate way to determine the relative productivity of populations year to year is by using genetic differentiation among populations to apportion mixed stock fishery harvests [that is, mixed stock analysis (MSA)] (). Genetically based MSAs, however, do not yield migratory information of individuals and are often limited to apportioning harvests to coarse spatial scales. For example, harvests of Bristol Bay sockeye salmon are apportioned to this region’s nine major watersheds (). Even after large-scale SNP (single-nucleotide polymorphism) discovery efforts (), Bristol Bay Chinook salmon (Oncorhynchus tshawytscha) populations are indistinguishable from one another and from other large adjacent watersheds (the lower Kuskokwim River). Nonetheless, the fine-scale level of biodiversity has been recognized as an important aspect of sustainable fisheries (), and questions remain regarding how it responds to perturbation. Tools are needed to determine finer-scale levels of population structure in fishery harvests to develop time series of fine-scale productivity.

Trace element–to–calcium ratios recorded in otoliths, which grow via metabolically inert concentric rings of CaCO3, are another effective MSA tool, especially in regions with shallow genetic structure (). Otoliths also represent lifelong chemical records of fish. Using this approach to develop time series, however, has been limited by a fundamental assumption of all provenance studies, regardless of tracer: that the tracers used remain temporally stable (). This assumption is rarely met with trace elements. Confounding environmental and physiological effects on how trace elements get incorporated into otoliths () and their typically high temporal variability render evaluation of temporal variability throughout habitats difficult. This is especially true for species with large migration and dispersal ranges (for example, anadromous fish). Thus, the trace element approach must be cohort-specific, requiring researchers to match years of baselines to subsequent years of mixed stock harvests ().

Unlike trace elements, strontium (Sr) isotope (87Sr/86Sr) ratios in otoliths (regardless of species, environment, or physiology) directly reflect ambient environments (). This attribute simplifies interpretations of 87Sr/86Sr-based provenance studies and allows evaluation of 87Sr/86Sr temporal variability throughout study regions via sedentary organisms (). When included in provenance studies, 87Sr/86Sr ratios generally provide the most discriminatory power to determine natal sources of individuals (). Last, 87Sr/86Sr variability scales with geologic heterogeneity (), existing at multiple spatial scales, include broad (among watersheds) and fine (within watersheds) geographic scales.

RESULTS AND DISCUSSION

We developed a 87Sr/86Sr-based MSA model for the Nushagak River [the third largest river in Western Alaska (34,700 km2), the largest of the nine major watersheds flowing into Bristol Bay, and a major producer of wild Pacific salmon] to apportion a mixed stock harvest of Chinook salmon with respect to two ecological dimensions: (i) natal origins and (ii) freshwater life histories. Geologic heterogeneity within the Nushagak River gives rise to geographic patterns in 87Sr/86Sr ratios of river waters, which are temporally stable at sub- and interannual time scales (). Our MSA model was built using water data from throughout the Nushagak River (n = 95 waters) (), producing seven strontium isotopic groups (SIGs) (Fig. 1). To independently validate the performance of the model to accurately classify natal origins of salmon, we used otolith samples of known origin [juvenile Chinook salmon, n = 153; slimy sculpin (Cottus cognatus), n = 33] (). These two tests yielded overall classification accuracies of 90 and 88%, respectively (table S1). To evaluate if juvenile Chinook salmon otoliths reflected the 87Sr/86Sr composition of waters at capture sites and to assess juvenile movements before capture, we conducted (i) a 87Sr/86Sr water-otolith regression, which yielded a slope = 0.983 ± 0.017 (2 SE) and r2 = 0.99 (fig. S1), and (ii) full life history profiles of a subsample at each capture site. These analyses indicated that some juveniles moved into capture sites before collection (). Thus, using water data to build the MSA was more appropriate than using otolith data, the more common approach ().

Fig. 1

Map of the Nushagak River, SIGs, and sampling sites for juvenile Chinook salmon and waters (black filled circles) ().

Our MSA model apportioned a fishery harvest to seven isotopically different sub-basins within the Nushagak River (Fig. 2A); this equates to a sevenfold increase in the current resolution for determining the relative productivity of populations from this major watershed (). Of the 255 harvested adult Chinook salmon analyzed, our MSA model predicted a total of 71% originated from three of the SIGs (SIG1, 22%; SIG2, 22%; SIG6, 27%) (see Fig. 1 for SIG locations). Our MSA results for the overall proportions of SIG membership changed only slightly (table S2) when we set a probability threshold for group membership of >80 or >70%. Thus, our overall results were 90% accurate, as indicated by our validation tests and posterior probabilities of SIG membership.

Fig. 2

Chinook salmon 87Sr/86Sr-based MSA.

(A) Proportion of each SIG in the incidental catch of adult Chinook salmon in Nushagak Bay in 2011. (B) Day-specific proportions of each SIG during the 3-day fishing period.

The proportion of harvests of Chinook salmon in each SIG scaled strongly (correlation coefficients >0.76) with estimates of potential habitat amount (stream length and basin area within each SIG) (fig. S2 and table S3). Thus, in 2011, the number of fish produced from each SIG was largely related to differences in habitat quantity ().

The proportion of harvested individuals from respective SIGs changed over a 3-day period and exhibited a geographic pattern in terms of day-specific MSAs (Fig. 2B). For example, the proportion of fish from SIG1 and SIG2 increased steadily over the 3-day period (by ~13%). Correspondingly, SIG5, SIG6, and SIG7 decreased (for example, SIG6 by ~8%). Geographically, these groups of SIGs represent the two primary upper regions of the Nushagak River (the Mulchatna and the Upper Nushagak Rivers).

Four distinct 87Sr/86Sr life history patterns during freshwater residence were evident in the fishery harvest (Fig. 3 and table S4A). First, individuals indicating site fidelity (no apparent change in 87Sr/86Sr) during freshwater residence made up 72% of the catch. Second, 7% reared in a different 87Sr/86Sr environment than their natal origin. Third, 17% originated and reared in the same 87Sr/86Sr environment, but showed short forays (occurring over an otolith distance of <100 μm just before marine migration, Fig. 3) toward lower river ratios (the migration corridor through which all upstream fish must swim). Last, 4% of the individuals reared in a nonnatal environment showed a lower river foray just before their migration to the ocean. Using the MSA model, we were also able to infer actual migration routes, refine natal origins, and determine habitat use of the individuals that reared in a different SIG than their natal origin (table S4B) ().

Fig. 3

Representative 87Sr/86Sr life history profiles from each SIG of adult Chinook salmon.

Colors of profiles correspond to Fig. 1. Also shown is an adult, which originated from SIG1 but was reared in SIG2.

Surprisingly, since it was first shown that calcified structures of fish reflect ambient water 87Sr/86Sr ratios (, ), 87Sr/86Sr-based MSAs of fishery harvests of anadromous species have been uncommon (, ). Previous attempts were unable to discern group membership at the sub-basin scale (), assess freshwater life history variation (, ), evaluate 87Sr/86Sr temporal variability throughout all natal sources (), or encompass all potential natal sources (). These limitations are fundamental assumptions of provenance studies and imperative to building MSA models, which can be used to develop time series of the relative production of populations at the fine-scale level.

Here, we simultaneously determined the relative production of sub-basins and life histories by conducting an 87Sr/86Sr-based MSA of fishery harvests. Most importantly, this MSA model can be used to develop time series of production at the scale of individual populations (). In 2011, we determined that (i) three SIGs produced >70% of the catch; (ii) the majority of recruitment (>70%) exhibited site fidelity during freshwater residence before marine migration; (iii) the number of fish produced by an SIG scaled with the amount of habitat; and (iv) eastern SIGs contributed more fish earlier during the 3-day period for which fishery samples were collected, whereas western SIGs contributed more fish later.

Additional years of data, which assess these same parameters, will provide valuable insights into the role of individual sub-basins and life history patterns on the overall Nushagak River stock productivity and stability, and how these populations respond to changes in the environment. Comparing the MSAs of out-migrating salmon (smolts) to the MSAs of adult recruitment will provide a new way to assess how fine-scale population structure and freshwater life history patterns affect marine survival and overall recruitment. Although we focused on Chinook salmon (), our MSA model is also applicable to other fish species (), including sockeye and coho (O. nerka) salmon that rear in freshwater before going to sea, because of the conservative nature of 87Sr/86Sr ratios.

Wild salmon and salmon-based cultures have been decimated worldwide. Alaska is home to some of the world’s last thriving wild salmon runs, productive wild salmon fisheries (), and indigenous cultures where salmon have remained culturally, nutritionally, and economically central for at least 4000 years (for example, the Yup’ik and Dena’ina of Bristol Bay, Alaska) (). Developing metrics able to assess the relative production of natal habitats and of different life histories through time, and how these ecological dimensions respond to perturbation (for example, commercial fisheries, industrial development, and climate change), will be integral to the conservation of wild salmon. Strontium isotopic heterogeneity within watersheds provides a viable framework to develop such metrics. The issues of provenance and migration pose vexing conservation problems for all mobile species, not just anadromous fishes. Coupling spatially and temporally robust 87Sr/86Sr baselines with sequentially growing animal tissues (for example, mammalian teeth or bird feathers) offers a sound framework for identifying crucial habitats to inform conservation planning.

MATERIALS AND METHODS

The fishery

Here, we focused on the incidental catch of Chinook salmon during the Bristol Bay sockeye salmon commercial fishery in the Nushagak District because (i) no methods for discerning this species stock composition within this watershed currently exist; (ii) its annual bycatch during the sockeye fishery is often >10,000 individuals (, ); (iii) juveniles from Subarctic populations spend one full year in the freshwater environment before out-migration to sea () and, thus, adopt accurate 87Sr/86Sr ratios from their natal habitats (); and (iv) Chinook salmon populations have shown dramatic changes across Alaska in the past decade, and tools are needed to accurately discern their fine-scale population structure (, ). Current management practices deliberately restrict the commercial sockeye salmon fishery to take place within the natural boundaries of Nushagak Bay (fig. S3), such that the fishery targets primarily fish ultimately bound for Nushagak Bay watersheds. Because the Nushagak River Chinook salmon run is so large (>200,000, 20-year average) compared to the small runs of proximal watersheds within Nushagak Bay (for example, the largest of which is likely the Muklung River within the upper Wood River system, ~400 total) (), incidentally caught Chinook salmon within Nushagak Bay can be assumed to be bound for the Nushagak River. Further, Nushagak River Chinook salmon make up about 90% of Chinook salmon in the Bristol Bay region (). By far, the largest nearby runs of Chinook salmon are bound for rivers in the Togiak District (fig. S3), where the total district-wide run is about 20,000 fish (20-year average) (). Because of the migratory trajectory of returning Chinook salmon to Bristol Bay rivers, the Togiak District is the only district “downstream” of the Nushagak River. Thus, we think that strays from other systems are likely to have minimal impact on our conclusions regarding the provenance of fish captured in Nushagak Bay because: (i) of the size of the Nushagak River run relative to nearby systems, and (ii) harvests are directed in the terminal area of the river mouth. These are important considerations because 87Sr/86Sr ratios of the other Bristol Bay watersheds overlap the lower end of the isotopic range (0.704 to 0.706) within the Nushagak River ().

Otolith collections

Juvenile salmon (n = 153 total) were captured using minnow traps in the autumns of 2011, 2012, and 2013 at 24 locations throughout the Nushagak River watershed (Fig. 1). Sagittal otoliths were dissected in the field and stored dry in polypropylene tubes until sectioning and analysis. Trapping locations included all of the major tributaries and main river channels of the primary branches of the Nushagak River (for example, Mulchatna, Nuyakuk, Upper Nushagak, and Lower Nushagak Rivers) (Fig. 1). We previously generated water 87Sr/86Sr ratios corresponding to these trapping sites, in addition to many more sites (n = 95 total) (black filled circles, Fig. 1) (). Otoliths were collected from adult Chinook salmon (n = 255 total) incidentally caught during the 2011 commercial sockeye salmon fishery conducted in the Nushagak Fishing District. Collections occurred during a 3-day period in late June (27 to 29 June 2011) during the peak of the Chinook salmon run () and included catches from different gear types (that is, both set and drift gillnets) and multiple boats. Fish were collected under Alaska Department of Fish and Game Fish Resource Permit numbers SF2011-236, SF2012-231, and SF2013-255 and Institutional Animal Care and Use Committee protocol number 178401-11.

87Sr/86Sr ratio analyses of otoliths

87Sr/86Sr ratios of otoliths were measured using laser ablation (LA) (193-nm excimer laser, Photo Machines) multicollector inductively coupled plasma mass spectrometry (MC-ICPMS) (Thermo Scientific, high resolution Neptune) at the University of Utah, Department of Geology and Geophysics ICPMS laboratory. Sectioning and mounting methods followed those outlined by Donohoe and Zimmerman (). Juveniles were sectioned in the sagittal plane, whereas adults were sectioned in the transverse plane. Before isotopic analyses, otoliths were sonicated for 5 min in MilliQ water, rinsed, and dried in a laminar flow hood. LA transects were implemented using a 31.4-μm-diameter circle for juveniles and a 53.1-μm circle for adults with a pulse rate of 10 Hz, a scan rate of 2 μm/s, and a laser energy of 45% (that is, 3.23 J/cm2). Counts per second of 88Sr, 87Sr, 86Sr, 85Rb, and 83Kr were measured at the Faraday cups with an integration time of 0.524 s, corresponding to one cycle. Before each ablation, transect background intensities (V) of each isotope were measured for 120 cycles, and the mean was used as a blank correction during sample analyses. Analytical accuracy of LA data was evaluated by measuring the 87Sr/86Sr ratio of a modern marine shell during and after each LA run (n = 10 runs; 9 to 12 shell analyses per run). LA run-means of shell ratios ranged from 0.70922 to 0.70925 with an average weighted 2SE of ±0.00002, which is consistent with our MC-ICPMS solution analyses of the same shell () and the global marine value (0.70918 ± 0.00006 2SD) (). We therefore applied no standard correction to otolith data. 87Sr/86Sr ratios of otolith samples were corrected for mass bias using an exponential law and for isobaric interference.

Adult Chinook salmon otoliths were measured perpendicular to growth axis, and each transect encompassed the core, entire freshwater residence, and migration into the marine environment. We determined the natal origin, rearing area (if different from natal origin), and overall freshwater movement pattern. We considered freshwater residence within each otolith to be the region between the distal extent of the core (~250 μm from primordia) and the distal extent of the first annulus (that is, before marine migration). Specifically, we determined the freshwater residence portion of each transect by inspecting the 87Sr/86Sr ratio profile and the corresponding 88Sr intensity (V) profile and superimposing transects on respective otolith images (taken in reflected light). For each otolith 87Sr/86Sr life history profile, we fitted a generalized additive model (GAM) using the MGCV package in R (http://cran.r-project.org/), which uses penalized iteratively reweighted least squares (P-IRLS) to maximize goodness-of-fit and general cross-validation (GCV) methods to minimize overfitting of each 87Sr/86Sr profile (). We also calculated Bayesian 95% confidence intervals along each transect, which allowed us to inspect if and where changes in the 87Sr/86Sr ratio occurred (). Details on the use of GAMs for otolith LA data are provided elsewhere (). If there was no change in the 87Sr/86Sr ratio during freshwater residence, we determined the natal origin to be represented by the entire freshwater residence period. If, however, there was a change during freshwater residence, we determined the natal origin to be between the distal extent of the core and the inflection point of the first 87Sr/86Sr change. The rearing signal was determined to be between the end of the natal origin signal and the migration to ocean.

The 87Sr/86Sr ratios of all juvenile Chinook salmon otoliths were measured ~50 to 80 μm from the edge of the otolith and parallel to growth axis for 200 μm. This region was targeted because we assumed it to be representative of the ambient environment experienced by each individual just before being captured. It was also sufficiently distal to the maternally influenced otolith core region (, ). Additionally, to investigate if juvenile salmon captured at a particular site also originated at that site (that is, the freshwater movements before being captured), we performed additional LA transects perpendicular to growth axis from the otolith core to the edge on a random subsample from each trapping site (n = 3 individuals per site). Because we measured from the core to the edge of each otolith of these subsampled juveniles, we also computed the 87Sr/86Sr ratios for the core (the middle ~100 μm), a natal origin signal, and a rearing signal (if indicative of movement). Using the former two ratios, we calculated the 87Sr/86Sr mass balance for Sr within the otolith core between marine-derived and freshwater-derived sources of Sr (fig. S4) (), assuming a two end-member mixing model defined by the equation below:where P is the proportion of Sr derived from marine sources via the egg, 87Sr/86Srcore is the measured ratio in the otolith core, 87Sr/86Srnatal is the measured ratio from the natal origin otolith region, and 87Sr/86Srmarine is the global marine value (0.70918).

Classification statistics

To conduct the MSA, we used discriminant function analysis (DFA) to predict the natal origins and freshwater movement patterns of adult Chinook salmon caught in Nushagak Bay using the MASS package in R (http://cran.r-project.org/). We trained the DFA using a river water data set (n = 95) of the Nushagak River watershed (). To demonstrate that using water data to train the DFA model was appropriate, we conducted a regression of ambient water versus otolith 87Sr/86Sr ratios. Because water 87Sr/86Sr ratios exhibited a strong geographic pattern throughout the Nushagak River (), we determined the SIGs for the DFA on the basis of similarity in both 87Sr/86Sr ratios and geography (Fig. 1 and fig. S1A). In the DFA, each SIG was given equal a priori probabilities (that is, 1/7). To validate the DFA model, we conducted two independent tests of model accuracy using measurements of otoliths from known origin: (i) juvenile Chinook salmon (reported here) and (ii) slimy sculpin indicative of site-specific temporal variation due to their low dispersal (). We then used the DFA model to predict the natal origins and movement patterns of adult Chinook salmon (n = 255) incidentally caught during the 2011 commercial sockeye salmon fishery. For each unknown individual, the posterior probabilities of SIG membership computed by DFA (that is, SIGs 1 to 7) summed to 1, and individuals were assigned to the SIG with the highest probability of SIG membership. To calculate the uncertainty for each SIG proportion computed by DFA (Fig. 2A and table S2), we used the mean maximum posterior probability of all individuals within a said SIG, where the uncertainty of a SIG was its error rate (that is, 1 − Pmax, where Pmax is the mean maximum posterior probability of SIG membership).