Non-stationary and interactive effects of climate and competition on pink salmon productivity.

Pacific salmon (Oncorhynchus spp.) are exposed to increased environmental change and multiple human stressors. To anticipate future impacts of global change and to improve sustainable resource management, it is critical to understand how wild salmon populations respond to stressors associated with human-caused changes such as climate warming and ocean acidification, as well as competition in the ocean, which is intensified by the large-scale production and release of hatchery reared salmon. Pink salmon (O. gorbuscha) are a keystone species in the North Pacific Ocean and support highly valuable commercial fisheries. We investigated the joint effects of changes in ocean conditions and salmon abundances on the productivity of wild pink salmon. Our analysis focused on Prince William Sound in Alaska, because the region accounts for ~50% of the global production of hatchery pink salmon with local hatcheries releasing 600-700 million pink salmon fry annually. Using 60 years of data on wild pink salmon abundances, hatchery releases, and ecological conditions in the ocean, we find evidence that hatchery pink salmon releases negatively affect wild pink salmon productivity, likely through competition between wild and hatchery juveniles in nearshore marine habitats. We find no evidence for effects of ocean acidification on pink salmon productivity. However, a change in the leading mode of North Pacific climate in 1988-1989 weakened the temperature-productivity relationship and altered the strength of intraspecific density dependence. Therefore, our results suggest non-stationary (i.e., time varying) and interactive effects of ocean climate and competition on pink salmon productivity. Our findings further highlight the need for salmon management to consider potential adverse effects of large-scale hatchery production within the context of ocean change.

INTRODUCTION

Pacific salmon (Oncorhynchus spp.) inhabiting the North Pacific Ocean and adjacent freshwater habitats are exposed to natural environmental changes and anthropogenic stressors. The joint impacts of stressors such as changes in climate, ocean acidification, and industrial scale hatchery production (the release of artificially propagated juvenile salmon into the ocean) on the productivity of wild salmon populations have received increasing attention in recent years (Cline et al., 2019; Connors et al., 2020; Cunningham et al., 2018). Understanding how multiple stressors interact to influence wild salmon populations is critical for anticipating future responses to global change and for improved resource management.

Pink salmon (O. gorbuscha) are the most abundant salmon species in the North Pacific Ocean, and their numbers have more than doubled over the past 50 years due to strong production of wild and hatchery populations (Ruggerone & Irvine, 2018). The majority of hatchery produced pink salmon originate in the state of Alaska, where they were used in the 1970s as a management tool to supplement low abundances of wild populations. Pink salmon have trophic impacts on marine ecosystem components, including phytoplankton and zooplankton (Batten et al., 2018; Shiomoto et al., 1997), other salmonids (Cline et al., 2019; Kendall et al., 2020; Ruggerone & Connors, 2015; Ruggerone & Nielsen, 2004), and predators such as seabirds (Springer & van Vliet, 2014; Springer et al., 2018; Toge et al., 2011). While pink salmon are considered a keystone species in the North Pacific due to their broad impacts on other organisms, less is known about factors that limit population productivity and abundance of pink salmon. Understanding these factors may help inform management decisions in Alaska and beyond, including other parts of the world where pink salmon are increasing in abundance and considered invasive (Nielsen et al., 2020; Sandlund et al., 2019).

About 5 billion hatchery salmon, primarily chum (O. keta) and pink salmon, are released into the North Pacific Ocean each year (Ruggerone & Irvine, 2018). Hatchery fish account for roughly 40% of the total salmon biomass in the North Pacific Ocean (Ruggerone & Irvine, 2018). Supplementation from hatcheries can affect wild populations in multiple and complex ways. Hatchery programs were developed primarily to mitigate declines in wild populations; however, concerns over adverse genetic or competitive effects of hatcheries on wild salmon have been raised repeatedly (Araki et al., 2007; Hilborn, 1992; Jasper et al., 2013; Naish & Hard, 2008; Waples, 1991). Although the majority of hatchery reared pink salmon returning as adults are harvested in ocean fisheries, several million fish that spawn in regional streams may interbreed with wild‐origin fish (Knudsen et al., 2021). Some hatchery programs might replace rather than augment wild production due to increased competition for resources, as has previously been suggested for pink salmon in Alaska (Hilborn & Eggers, 2000; but see Wertheimer et al., 2001). The potential for negative impacts of increasing abundances of hatchery reared salmon in the ocean has led to calls for an open dialogue on the number of hatchery fish being released each year (Connors et al., 2020; Holt et al., 2008).

Long‐term trends in environmental conditions may pose substantial threats to salmon survival. Ocean acidification has been identified as a potential stressor for aquatic species globally, including in the Northeast Pacific (Feely et al., 2016; Orr et al., 2005), for example due to the dissolution of pteropods that are a primary food source for pink salmon (Bednarsek et al., 2021). Experimental studies have also found negative effects of ocean acidification on the growth and behavior of coho (O. kisutch) and pink salmon (Ou et al., 2015; Williams et al., 2019).

Furthermore, climate regime shifts have occurred in the Northeast Pacific Ocean in 1976–1977 and 1988–1989 (Hare & Mantua, 2000; Irvine & Fukuwaka, 2011). While the 1976–1977 event was characterized by a change in sign of the Pacific Decadal Oscillation Index (PDO), the 1988–1989 shift involved a change in the leading mode of North Pacific climate variability, from PDO‐like to more North Pacific Gyre Oscillation (NPGO)‐like (Yeh et al., 2011). Recent work has shown that the effects of ocean temperature on salmon productivity, including on pink salmon, have changed around 1988–1989, with generally weaker links after the 1988–1989 event (Litzow et al., 2019, 2020). Non‐stationary relationships among physical and biological variables indicate that the effects of environmental conditions, such as temperature, on population processes, such as recruitment, can vary over time in intensity and/or direction.

In contrast to other species of Pacific salmon, pink salmon have a fixed 2‐year life history such that all fish return to reproduce 2 years after eggs were deposited in the gravel (Ricker, 1962). Odd‐ and even‐year lineages that spawn in the same rivers are thus reproductively isolated and genetically distinct (Aspinwall, 1974; Beacham et al., 2012). One of the broodlines is typically more abundant than the other, causing 2‐year cycles in pink salmon returns. This pattern of broodline dominance can shift over time and varies among regions along the west coast of North America (Irvine et al., 2014; Krkosek et al., 2011). Cycles in return abundances likely result from a combination of stochastic recruitment and density‐dependent interactions between broodlines, such as competition, cannibalism, or disease transmission (Krkosek et al., 2011). Competition for resources in pink salmon may arise from interactions within broodlines, between broodlines, or with other salmon species, including both wild and hatchery fish.

We investigated the combined effects of changes in ocean conditions and salmon abundances (pink, chum, and sockeye salmon) on the productivity of wild pink salmon in Prince William Sound (PWS), Alaska. Each year about 600–700 million hatchery pink salmon fry are released into PWS, which constitutes roughly 50% of the global pink salmon hatchery production in recent years (NPAFC, 2020; Stopha, 2019). We used 60 years of available time‐series data (1960–2019) on wild pink salmon abundances based on catch and escapement records from PWS to study how wild pink salmon productivity is affected by changes in ocean temperature and acidification, as well as competition among pink salmon, including hatchery reared fish, and interspecific competition with other species of Pacific salmon.

MATERIALS AND METHODS

Pink salmon ecology

Pink salmon eggs are laid in late summer to fall in coastal rivers, streams, and brackish estuaries (year of spawning is referred to as “brood year,” BY). They are semelparous and die after spawning (Quinn, 2005). Juveniles do not reside in fresh water for an extended period and migrate to sea during their first spring (BY+1). Prince William Sound pink salmon spend the summer in nearshore habitats of PWS before entering the Gulf of Alaska (GoA) by late summer or early fall and migrate southwest in association with the Alaska Current and the Alaska Coastal Current (Armstrong et al., 2008). The typical migration pattern of PWS pink salmon and the general ocean distribution of other Pacific salmon is illustrated in Figure 1. When adults return to PWS the following summer (BY+2) they may overlap with juveniles of the other broodline in coastal environments before entering their natal streams to spawn. Pink salmon primarily feed on zooplankton, pteropods, squid, and other fishes, and diet overlap is largest with sockeye salmon (O. nerka) and chum salmon (O. keta) (Johnson & Schindler, 2009; Kaeriyama et al., 2000).

FIGURE 1

Map of the Northeast Pacific Ocean. Shown are the typical migration routes of pink salmon originating from Prince William Sound (PWS, black dotted line) and other west coast Pacific salmon (thick gray line). Filled blue circles indicate hypothesized effects of competition during the pink salmon life cycle. Thick blue line shows the maximum distributions of North American pink salmon

Hatchery pink salmon

The State of Alaska hatchery program was developed in the early 1970s in response to declining salmon abundances, with the goal of supplementing wild stock abundance for the public benefit (Wilson, 2020). Hatchery reared pink salmon are currently released into PWS as fry by four major hatcheries operated by the Prince William Sound Aquaculture Corporation (PWSAC: Armin F. Koernig, Cannery Creek, Wally Noerenberg) and the Valdez Fisheries Development Association (VFDA: Solomon Gulch). Another facility, Main Bay Hatchery (PWSAC), released pink salmon until 1989. The PWSAC operates the largest hatchery program in North America (Hilborn & Eggers, 2000). These hatcheries are located in close geographic proximity to wild spawning areas, and hatchery reared juveniles occupy similar nearshore marine habitats as juveniles from wild populations (Boldt & Haldorson, 2004; Moss et al., 2005). Like wild pink salmon, hatchery fish attain most of their growth in the ocean, and return during the year following their ocean entry.

Data sources

We obtained estimates of escapement, harvest, and total run size of PWS wild pink salmon, as well as release and harvest statistics of PWS hatchery pink salmon from the Alaska Department of Fish and Game (ADF&G) (Figure 2a, d). Pink salmon harvest statistics were downloaded from the ADF&G electronic fish ticket database (ADF&G, 2019), where hatchery contributions are assessed by examining otoliths for hatchery thermal marks and expanding the estimates to the entire catch (Haught et al., 2019). Hatchery release numbers were obtained from the ADF&G Mark, Tag, and Age Lab (https://mtalab.adfg.alaska.gov/cwt/reports/hatcheryrelease.aspx). Estimates of spawning escapement that produce PWS wild pink salmon were derived from annual aerial surveys of index streams (Russell & Haught, 2020). Specifically, area‐under‐the‐curve (AUC) estimates of daily stream counts were divided by stream‐specific average stream life of the fish and expanded to account for mean observer efficiency (0.436) (Bue et al., 1998) and the proportion of escapement in non‐aerial index streams (0.20) to estimate total spawner escapement. Our escapement estimates are identical to the methodology used by the Alaska Department of Fish and Game, but these approaches do not account for time‐varying observation errors or changes in methodology. Observation errors may depend on conditions during surveys such as tides, visibility, survey frequency, and variation in observer efficiency (Fried et al., 1998). Methodological changes include changes in the number of streams assessed over time and changes to escapement goals (Haught et al., 2017). We assume that effects of observation errors and changes in methodology over time are relatively small compared with the year‐to‐year variation in the numbers of pink salmon spawning in PWS. Similarly, aerial escapement surveys do not discern between hatchery and wild fish, and the fraction of hatchery‐origin fish on the spawning grounds varies over time and space (Brenner et al., 2012; Joyce & Evans, 1999; Knudsen et al., 2021). A recent study estimated the contribution of hatchery fish to spawning streams at 5%–15%, based on a subset of streams sampled over a period of 3 years (Knudsen et al., 2021); however, straying proportions for individual streams can be substantially higher (Brenner et al., 2012; Joyce & Evans, 1999).

FIGURE 2

Time series of response and predictor variables. Shown are time series of (a) wild run size as escapement (dark) plus harvest (light) of even‐year (purple) and odd‐year (yellow) pink salmon in PWS, (b) recruits per spawner and ln(recruits/spawner) for wild pink salmon in PWS (gray line indicates replacement at one recruit/spawner), (c) spring sea surface temperature (SST) anomaly in the eastern GoA, and (d) hatchery releases of juvenile pink salmon into PWS (purple, solid line), harvests of adult hatchery pink salmon returning to PWS (blue line and circles), and competitor abundance (sockeye and chum salmon) in the GoA (gray, dashed line)

Temperature time series were obtained from the Extended Reconstructed Sea Surface Temperature dataset (ERSSTv4, Huang et al., 2015) for the years 1960–2019. Winter and spring sea surface temperature (SST) anomalies in the eastern and western GoA were calculated as averages for the months November to March and April to June, respectively (Litzow et al., 2020). Reconstructions of pCO2 and pH in the GoA and PWS, which are used as proxies for ocean acidification, were taken from a recently developed regional ocean biogeochemical model of the GoA that were available for the time period 1980–2013 (Hauri et al., 2020). We averaged over three spatial areas representing ocean conditions in PWS, the Gulf of Alaska south of PWS, and the central subpolar gyre (Hauri et al., 2021), defined by latitudes and longitudes 60°–61°N and 146°–148°W, 57.5°–59°N and 148°–151°W, and 55°–57°N and 147°–150°W, respectively. In addition, we used reconstructions of the Northern Gulf of Alaska Oscillation (NGAO), an index of sea surface height variability that is directly linked to changes in upwelling of nitrate and dissolved inorganic carbon and thus ocean acidification in the northern Gulf of Alaska (Hauri et al., 2021). Finally, abundances of other Pacific salmon, including sockeye (O. nerka) and chum salmon (O. keta), were taken from Ruggerone and Irvine (2018) for the return years 1960–2015, and Ruggerone et al. (2021) for the return years 2016–2019.

Statistical analyses

Our analysis is based on the Ricker model (Ricker, 1954) that describes the relationship between spawner abundance, that is, the number of fish that escape the fisheries and are assumed to reach the spawning grounds, and subsequent recruitment. The linearized version of the Ricker function is:

Recruitment can be measured at juvenile or adult stages and here we use the latter (recruitment defined as catch plus escapement in return year ), is spawner abundance 2 years prior, is the intercept, is the rate at which population productivity declines with spawner abundance (), and are normally distributed errors. Here is the initial slope of the stock–recruit curve which reflects the number of recruits per spawner at low spawner abundance. This linearized Ricker model can be extended by incorporating additional explanatory variables to model changes in productivity over time: where are the regression coefficients and are the covariate time series lagged relative to return year when the covariate is hypothesized to affect pink salmon productivity. We thus used the natural logarithm of recruits per spawner as our metric of population productivity. An alternative and commonly used metric of population productivity is the residuals of a Ricker stock–recruitment fit. As a sensitivity analysis, we performed a second round of modeling, fitting a Ricker model to the data with different stock–recruit relationships in each broodline. We subsequently compared the model residuals to the ln(recruits/spawner) time series. These two time series showed a high correlation (Pearson correlation of 0.92, Figure S1), and using Ricker residuals as an alternative metric of productivity would yield similar results. Because correlative statistical analyses carry the risk of identifying spurious correlations, we limited our full model to only incorporate covariates with hypothesized mechanisms that may affect pink salmon growth or survival (Table 1). All covariate time series were standardized to a mean of zero and one standard deviation for the analysis. Centering ensures that main effects are biologically interpretable in the presence of interactions or polynomials, and standardization of covariates results in standardized slopes that are comparable within and between models (Schielzeth, 2010).

TABLE 1

Explanatory variables considered in the model. Shown are names, species or season, rearing type or location, spatial scale, unit, and hypothesized mechanism for each covariate (including the lag in years relative to the year of wild return)

Covariate	Species/Season	Type/Location	Spatial scale	Unit	Mechanism
Spawner abundance1	Pink	Wild	Local/Regional	Millions of fish	Competition for spawning sites or juvenile resources (spawners two years prior)
Hatchery releases same broodline2	Pink	Hatchery	Regional	Millions of fish	Competition among juveniles during first summer in PWS (previous year)
Total run size other broodline1,3	Pink	Wild and hatchery	Regional	Millions of fish	Competition between juveniles and returning adults in PWS (previous year)
Abundance of competitors4	Sockeye + chum	Wild and hatchery	Continental or global	Millions of fish	Interspecific competition in the GoA or North Pacific (return year)
Sea surface temperature5	Winter or spring	GoA	Regional	°C	Changes in growth or survival in the ocean (previous year)
Northern Gulf of Alaska Oscillation6	Annual	Northern GoA	Regional	–	Changes in growth or survival in the ocean (previous year)
Carbon dioxide partial pressure or pH6	Annual	PWS or GoA	Regional	μatm (−)	Changes in growth or survival in the ocean (previous year)

Data sources: 1Russell & Haught, 2020, 2ADF&G Mark, Tag, and Age Lab, 3ADF&G electronic fish ticket database, 4Ruggerone & Irvine, 2018, Ruggerone et al., 2021, 5ERSSTv4, Huang et al., 2015, 6Hauri et al., 2020, 2021.

We considered the following continuous predictor variables in addition to spawner abundance: total run size of the other broodline in the previous year (wild plus hatchery adults returning to PWS during the year of wild juvenile outmigration), hatchery releases of the same broodline (released in the year of wild juvenile outmigration), competitor abundance in the GoA (sockeye and chum salmon from Kodiak to Washington State), competitor abundance in the North Pacific Ocean (sockeye and chum salmon from Asia and North America), proxies of ocean acidification in PWS, the GoA, and the subpolar gyre (pCO2, pH, NGAO), and spring or winter sea surface temperature (SST) in the eastern or western GoA. We also considered nonlinear effects of the biotic covariates (total run, hatchery releases, competitor abundance) and temperature by including quadratic terms in addition to the linear effects. Time series used in the selected covariate model are presented in Figure 2, and other time series considered in the model selection are presented in Figure S2.

We assumed that the even and odd pink salmon broodlines were separate populations as they are reproductively isolated and genetically distinct (Aspinwall, 1974; Beacham et al., 2012), with different relationships between population productivity and spawner abundance estimated for each broodline. We therefore included spawner abundance, broodline (odd/even), and a spawner abundance by broodline interaction as fixed terms. In addition, we considered a binary predictor representing the period before and after the 1988–1989 climate regime shift as a factor in interaction with SST. This interaction allows for different relationships between temperature and salmon productivity during each regime, as identified in previous studies on Pacific salmon (Litzow et al., 2019, 2020). We considered an interaction of regime with spawner abundance by broodline to allow for different relationships between productivity and spawners during each regime. We did not include a main effect for regime, which was strongly confounded with hatchery releases (hatchery production was limited in the 1970s and early 1980s). Confounded predictors can cause inaccurate parameterization and exclusion of important predictor variables during model selection (Graham, 2003). We used variance inflation factor analysis to assess multicollinearity among explanatory variables and to determine which variables could be included in the same model (threshold of 5, Zuur et al., 2010).

Because the initial model selection did not support the inclusion of ocean acidification proxies, which had relatively short time series (1980–2013), model selection was repeated using the full time series (1960–2019) of other explanatory variables. We used combined sockeye and chum salmon abundances as an index of competition, because time series of non‐PWS pink salmon in the GoA and North Pacific strongly covaried with PWS hatchery releases and total returns. The different spatial aggregates of competitor abundances or SST could not be included in the same model such that we ran the model selection with one of the time series at a time and compared models based on their AICc values. The full model included wild spawner abundance in interaction with broodline and regime, total return of the other broodline in the previous year, hatchery releases of the same broodline, an index of competitor abundance in the GoA or North Pacific, and spring or winter SST in the eastern or western GoA in interaction with regime. All sub‐models of the full model were tested as part of the model selection.

Diagnostics and interpretation

We performed AICc‐based model selection (Burnham & Anderson, 2002), where the model with the lowest AICc value was selected using the dredge function of the package MuMIn (v. 1.43.15) in R (R Core Team, 2020). We present results for the top model and the model selection, instead of using model averaging which may not be reliable when models contain interactions between explanatory variables (Cade, 2015). To evaluate whether these regression results were specific to the choice of model framework used, we also constructed models using random forests (using the R package randomForest, version 4.6–14; Breiman et al., 2006). Because random forest models are nonparametric, coefficients cannot be extracted for each predictor, but these models can be used to calculate the relative importance of each variable. Specifically, we calculated the percent increase in mean squared errors of out‐of‐sample predictions when excluding a predictor as our metric of variable importance (Liaw & Wiener, 2002).

Finally, we performed a cross‐validation to assess the ability of covariate models of different complexity to make out‐of‐sample predictions. The cross‐validation was performed by splitting the data into a training and a test dataset. Models were run on the training data (44 years, 75% of data) to get parameter estimates that were used to predict the remaining test data (14 years, 25% of data). This procedure was repeated 1000 times by randomly drawing the training and test datasets to achieve reliable estimates of the root mean squared prediction error (RMSE). The procedure was applied to all sub‐models of the selected model to assess whether any of the simpler models would produce better out‐of‐sample predictions compared to the selected model.

Model code is available at: https://github.com/janohlberger/PinkSalmon2021Ohlberger.git, and data are openly available in Zenodo at: https://doi.org/10.5281/zenodo.5780246.

RESULTS

Run sizes and escapements of wild pink salmon in PWS have varied considerably over the past 60 years (Figure 2). Return abundances of the even‐year broodline have varied between 1.9 and 27.9 million with a mean of 8.0 million fish, and abundances of the odd‐year broodline have varied between 3.5 and 43.2 million with a mean of 14.3 million fish. Population productivity measured as ln(recruits/spawner) has also varied substantially over time, between −0.45 and 2.2 with a mean of 0.74 (Figure 2).

We find evidence that the productivity of wild pink salmon is affected by multiple biotic and abiotic factors, including changes in ocean climate and competition, but not ocean acidification. The selected model included several explanatory variables in addition to the fixed model terms spawner abundance (), broodline (, odd/even), and their interaction (Table S1): where is the sea surface temperature in the GoA during spring of wild juvenile outmigration, is the total pink salmon return (wild and hatchery) in the year of wild juvenile outmigration, are hatchery releases of pink salmon into PWS in the year of wild juvenile outmigration,  is the competitor abundance in the GoA (sockeye and chum salmon) in the year of wild return, and is the ocean regime (before or after 1988–1989).

Negative density dependence, that is, lower productivity at higher spawner abundance, was found in both broodlines; however, the even‐year broodline showed stronger density dependence after the 1988–1989 regime shift (Figure 3a), whereas the odd‐year broodline showed weaker density dependence after the regime shift (Figure 3b). Overall, the odd‐year broodline had a slightly higher productivity than the even‐year broodline (Figure S3). The relationship between pink salmon productivity and temperature has shifted from positive and significant before to non‐significant after the 1988–1989 regime shift (Figure 3c). The total run of wild and hatchery pink salmon to PWS had a nonlinear effect on wild pink salmon productivity from the following brood year (Figure 3d), where the highest total runs were associated with low wild pink salmon productivity, indicating competitive interactions between returning adults and juveniles in nearshore habitats. Wild pink salmon productivity was negatively associated with pink salmon hatchery releases, suggesting adverse effects of hatchery production (Figure 3e). The relationship between productivity and the abundance of sockeye and chum salmon in the GoA was nonlinear (f), indicating a threshold effect where increasing interspecific competitor abundances show a positive association up to a threshold, above which additional competitors are associated with reduced productivity of wild pink salmon. Model parameter estimates and confidence intervals are provided in Table S1 and Figure 4.

FIGURE 3

Predicted partial effects on wild pink salmon productivity. Shown are effects of (a/b) spawner abundance for the even‐year and odd‐year broodlines during the previous (blue, dashed line) and current (red, solid line) ocean production regime, (c) spring SST in the eastern GoA in the year of juvenile outmigration for the previous (red) and current (blue) ocean regime, (d) total run size of wild plus hatchery pink salmon of the other (previous year) broodline to PWS, (e) hatchery releases into PWS in the year of juvenile migration, and (f) abundance of competitors (chum and sockeye salmon) in the GoA in the year of adult return. Polygons are 95% confidence intervals

FIGURE 4

Model parameter estimates. Shown are parameter estimates with 95% confidence intervals for all terms included in the selected covariate model

The model explained about 62% of the variance in ln(recruits/spawner) over time (Figure S4). For comparison, a model without any covariates (but allowing for independent stock–recruitment relationships by broodline) explained only 10% of the variance (Figure S5). Model diagnostics indicated normality and homogeneity of residuals (Figure S6). The percent increase in mean squared errors when excluding predictors from the model suggested that the most influential variable, besides the fixed terms spawner abundance and broodline, was SST anomaly, followed by hatchery releases, regime, total returns, and competitor abundance (Figure S7). The selected covariate model had the highest out‐of‐sample prediction ability, as indicated by the lowest root mean squared error (RMSE), compared to any of the simpler models that included fewer predictors (Figure S8). A sensitivity analysis of the model selection and parameter estimates that simulated varying degrees of observation error in recruits/spawner estimates further indicated that our findings are robust to low and moderate levels of observation error (Figure S9).

Scaling up hatchery releases of pink salmon from 0 to 700 million fish each year decreased the expected wild productivity under average conditions of the other predictors by ~55% for both broodlines (Figure 5). The absolute change in expected recruits per spawner depends on other factors such as the ocean regime and wild escapement (Figure S10). These estimates account for negative effects of hatchery juveniles released into PWS (competition within the same brood year), but do not account for potential impacts of returning adults (competition from the previous brood year), because the effect of total return was only apparent at very high abundances, and the effects of hatchery versus wild adults could not be distinguished. The resulting effects of hatchery releases on the expected wild return are presented in Figure 6. While the percent decline in the expected wild return when increasing hatchery releases is independent of the level of wild escapements, negative impacts of hatchery releases in terms of forgone wild returns (i.e. the expected additional wild return in the absence of hatchery releases) are greater when wild escapement is intermediate to high (Figure 6).

FIGURE 5

Predicted wild recruits per spawner as a function of hatchery releases. Shown are predictions for the even‐year (purple) and odd‐year (yellow) broodlines during the current regime while setting all other predictors (wild escapement, temperature, total return, competitor abundance) to median values. Thin lines show predictions beyond the highest number of hatchery fish that has been released into PWS (~700 million, gray vertical line). The gray horizontal line indicates population replacement at one recruit per spawner and the gray vertical line indicates the highest historical hatchery releases. Polygons are 95% confidence intervals

FIGURE 6

Predicted return of wild pink salmon as a function of wild spawner escapements and hatchery releases. Shown are predicted wild pink salmon returns (contour, millions of fish) for the (a) odd‐year broodline and (b) even‐year broodline during the current ocean regime. Other predictors were set to median values for that regime. Note that the even‐year broodline shows stronger density dependence compared to the odd‐year broodline during this regime. Predictions do not account for uncertainty in model parameter estimates

DISCUSSION

We find evidence that density dependence in pink salmon is caused by intraspecific interactions at the regional scale and that hatchery releases can negatively affect the productivity of wild populations. We also find that the 1988–1989 ocean regime shift was associated with changes in the temperature–productivity relationship and affected the strength of density dependence within pink salmon broodlines. Our results suggest that climate and competition can have non‐stationary and interactive effects on salmon and that potential adverse effects of large numbers of hatchery reared fish on wild populations should be considered in salmon management.

Hatchery releases of pink salmon in PWS have increased over time, especially during the 1980s, and currently amount to ~700 million fry annually such that returning hatchery fish now greatly outnumber returning wild pink salmon. While there is no evidence for a significant increase or decrease of pink salmon hatchery production in the near future, we predicted the effects of higher or lower hatchery releases on wild pink salmon productivity based on our model results. We estimate that the productivity of wild fish would be reduced to around one recruit per spawner (replacement) under average wild escapements and environmental conditions, if current hatchery releases were doubled. Conversely, the expected recruits per spawner of wild fish would increase by about 50% if hatchery releases were halved to ~350 million annually (from 2.33 to 3.50 and from 2.63 to 3.96 recruits per spawner for the odd and even broodlines, respectively, Figure 5).

Previous work suggested that hatchery production contributed to declining run sizes of wild pink salmon in PWS (Hilborn & Eggers, 2000, 2001; but see Wertheimer et al., 2001; Shuntov et al., 2017). Hatchery fish most likely reduce the productivity of wild populations due to competition for resources between hatchery and wild juveniles in nearshore marine habitats, which can negatively affect body growth and reduce survival of wild pink salmon, because fast‐growing juveniles experience higher survival rates (Cross et al., 2008, 2009; Moss et al., 2005). Hatchery production may have little or even no detectable effects on wild pink salmon populations in regions where hatchery programs are spatially segregated by larger distances from nearshore rearing habitats of wild juveniles. Accordingly, productivity trends of pink salmon in nearby Cook Inlet and Kodiak have differed from those in PWS (Amoroso et al., 2017; Hilborn & Eggers, 2000; Malick & Cox, 2016). Decreasing the number of hatchery pink salmon released into nearshore areas that constitute important habitat for wild juveniles may thus be an effective management strategy to reduce potential negative effects of hatchery supplementation on the productivity of wild pink salmon. The observation that the nearby Cook Inlet and Kodiak regions of Alaska have experienced different productivity trends supports the conclusion that the decline in productivity of PWS pink salmon was linked to hatchery releases, and less so to large‐scale changes in ocean conditions linked to the 1988–1989 regime shift. However, because the time series of hatchery releases is confounded with the regime shift, due to very low releases in the 1970s and early 1980s, uncertainty in the hatchery effect is likely larger than indicated by our model.

The question remains whether negative impacts of pink salmon releases on the productivity of wild populations are acceptable considering that hatchery production increases total abundances of pink salmon (Amoroso et al., 2017). Higher hatchery return abundances increase harvest opportunities, but do not appear to stabilize revenue of pink salmon fisheries in PWS (Ward et al., 2018). Furthermore, because Pacific salmon migrate long distances at sea, large‐scale hatchery production may have unintended adverse effects on other species of salmon originating from distant regions (e.g. Connors et al., 2020; Cunningham et al., 2018; Ruggerone et al., 2012). Hatchery salmon interbreeding with wild salmon can also affect the genetic composition and reproductive success of wild salmon (Naish & Hard, 2008; Waples et al., 2020). Hatchery pink salmon that stray into wild streams have lower reproductive success than wild fish (Lescak et al., In Press) and could exacerbate density dependence on the spawning grounds, which may bias estimates of wild recruits/spawner when escapement surveys do not discern between hatchery and wild fish. Quantifying the tradeoffs between industry performance in the fishery supported by the large hatchery program and productivity and abundance of wild salmon populations within and outside PWS are left for future extensions of this work.

In addition to hatchery releases in the year of wild juvenile outmigration, the model included the total run of hatchery and wild fish as a predictor variable, suggesting additional negative effects of returning adults on wild juveniles during their first summer at sea. While the effect was only apparent for the largest returns to PWS (over 50 million returning adults), it is conceivable that hatchery pink salmon compete with wild juveniles both when released into the ocean and as returning adults. Competition among juveniles that have recently entered the ocean and returning adults of the previous brood is consistent with other studies that have suggested that odd–even‐year cycles, as commonly observed in pink salmon returns, are caused by direct interactions between broodlines (Krkosek et al., 2011). Such direct interactions may also be responsible for the switch in odd–even dominance of western Kamchatka pink salmon (Ruggerone & Nielsen, 2009). Our results provide some evidence for such delayed density dependence and indicate that the more abundant broodline may be able to suppress the productivity of the other broodline at extremely high abundances.

While large abundances of pink salmon have previously been linked to declines in the growth and survival of other salmon species in the North Pacific Ocean, in particular sockeye salmon (Cline et al., 2019; Connors et al., 2020; Pyper & Peterman, 1999; Ruggerone & Connors, 2015; Ruggerone & Nielsen, 2004; Ward et al., 2017), competition from other Pacific salmon on the productivity of wild pink salmon has been studied much less, possibly because pink salmon are the most abundant salmonid. Our results suggest a nonlinear effect of sockeye and chum salmon abundance, indicating a potential threshold abundance above which competitors reduce the survival of wild pink salmon. The positive association at low competitor abundances might be a result of shared responses to environmental change at large spatial scales. The abundance of sockeye and chum salmon originating from Kodiak to Washington State was more closely related to pink salmon productivity than basin‐wide abundance including populations from western Alaska and Asia. However, conclusions about the spatial scale of density dependence need to be made with caution due to multicollinearity among candidate time series of competitor abundance, including different spatial aggregates, time lags, and species compositions.

The relationship between ocean temperature and pink salmon productivity has varied over time. Productivity was positively associated with temperature prior to the ocean regime shift, but the relationship was non‐significant after 1988–1989, suggesting that pink salmon productivity under the current regime is less dependent on temperature. This time dependence suggests caution when using observed temperature–productivity relationships to project possible future effects. Large‐scale climate indices such as the PDO and NPGO were considered in our initial modeling efforts but not included in the final models as they were not supported as important predictors of pink salmon productivity. It should also be noted that the Exxon Valdez oil spill in PWS occurred in 1989 during the NE Pacific regime shift. However, previous work found no evidence for a link between wild pink salmon productivity and the Exxon Valdez oil spill in PWS (Ward et al., 2017). Our finding that pink salmon productivity is linked to regional‐scale temperature conditions is also in line with previous studies (Mueter et al., 2002; Springer & van Vliet, 2014). Finally, similar non‐stationary climate–productivity relationships have been reported for other species of Pacific salmon, including sockeye and chum salmon, indicating that salmon productivity was positively related to sea surface temperature in the GoA in the 1970s and 1980s, but unrelated to SST after the 1988–1989 regime shift (Litzow et al., 2019, 2020).

The climate regime shift of 1988–1989 in the Northeast Pacific Ocean further affected the stock–recruitment relationship of the two broodlines in opposite ways. The effect of spawner abundance was much stronger after compared to before the regime shift for the even‐year broodline, suggesting intensified density dependence, whereas it was slightly weaker for the odd‐year broodline. These contrasting responses in the strength of density dependence may be linked to different environmental conditions. Pink salmon egg survival and fry growth are higher for the even‐year compared to the odd‐year broodline at cold incubation temperatures (Beacham & Murray, 1988). It has been suggested that increasing dominance of odd‐year brood lines along the North American coast might be caused by warming freshwater habitats (Irvine et al., 2014). While pink salmon in PWS have historically not shown any clear dominance of either broodline, it appears that an odd‐year dominance has emerged during the past two to three decades. Our results indicate that the climate regime shift altered the spawner–recruit relationships of the broodlines such that population productivity has declined for even‐year pink salmon but increased for odd‐year pink salmon at high spawner abundances. This interaction between density and climate might contribute to recent shifts toward increasing abundances of odd‐year pink salmon in PWS and possibly in other regions, as has been observed throughout the North American range (Irvine et al., 2014).

Factors other than those accounted for in this study can affect wild pink salmon productivity, including environmental conditions during freshwater rearing. For instance, a consequence of a warming climate in Pacific salmon ecosystems is loss of glacier coverage and the resulting changes in freshwater and estuarine habitats (Schoen et al., 2017). Glaciers in western North America are currently losing about 3% of their volume per year (Pitman et al., 2020). The PWS watershed still features significant glacier coverage (~18%), yet continued glacier loss in the region might affect wild pink salmon productivity and interact with other factors, such as those related to competition for resources in freshwater and estuarine habitats. These habitats are modified by glacier retreat via a variety of mechanisms that may increase or decrease the productivity of wild salmon depending on current thermal, biogeochemical, and hydrological conditions (Pitman et al., 2020). While we did not detect any links between variation in pCO2 or pH and PWS pink salmon productivity, this should not be interpreted as conclusive evidence for the lack of an effect, because the reconstructed time series that we used as a proxy for ocean acidification were relatively short and based on a regional ocean biogeochemical model, rather than observations. Furthermore, our analysis assumes that estimated spawner escapements represent true values. While these estimates are likely associated with error due to a number of factors, we believe that the values in this study are the best available as they correct for observer efficiency, stream life, and the proportion of escapement assessed by the aerial survey program (Fried et al., 1998). Observation error can cause bias in estimates of population productivity derived from stock–recruit analyses (Hilborn & Walters, 1992), and observation error in pink salmon escapement estimates can be considerable (Bue et al., 1998; Hilborn et al., 1999). However, our sensitivity analysis suggested that the findings of this study are relatively robust to low and moderate levels of observation error in recruits/spawner estimates.

Our findings highlight that the benefit of higher wild escapements for producing greater recruitment of wild pink salmon diminishes with increasing hatchery releases. Reduced productivity due to hatchery supplementation likely contributes to lower population resilience in wild pink salmon. Further increasing hatchery releases might therefore jeopardize the ability of wild populations to withstand other environmental and/or human stressors.

CONFLICT OF INTEREST

The authors have declared that no conflict of interests exist.

AUTHOR CONTRIBUTIONS

J.O., E.J.W., R.E.B., and M.E.H. designed the study; J.O., E.J.W., R.E.B., M.E.H., S.B.H., M.A.L., G.T.R., and C.H. curated the data; J.O. and E.J.W analyzed the data; J.O. wrote the first draft of the manuscript. All authors contributed to interpretation of results and writing of the manuscript.

Supplementary Material

Click here for additional data file.