(C) ProPublica.
This unaltered story was originally published at ProPublica.org. [1]
Licensed under creative commons by-nc-nd/3.0 [2]


Climate change threatens Chinook salmon throughout their life cycle

Author Name, ProPublica

2022-05

Data informing survival estimates

Spawner age and abundance estimates were compiled through large-scale collaborative efforts between states, tribes, and coordinating bodies94,95,96. Stage-specific survival estimates were obtained from multiple sources97,98,99,100, originating from tagging and detection records downloaded from PTAGIS.org, a regional collaborative database (Supplementary Table 2). We used detection records only of fish identified as wild from known population sources. No hatchery fish were released in any of the streams supporting populations included for this analysis. Environmental covariates used in the model included air temperature, streamflow and temperature, SST, and coastal upwelling (Supplementary Table 3).

For the smolt-to-adult analysis, we considered data for all out-migrating yearling spring/summer Chinook salmon tagged in the Snake River Basin detected from 2000 to 2015 at Bonneville Dam—the furthest downstream dam on the Columbia River. We marked a fish as having survived the marine stage if it was detected at Bonneville Dam (or farther upstream) as an adult. The data included (i) the last detection date at Bonneville Dam as juveniles, (ii) transportation classification (transported or in-river), and (iii) whether the fish was detected in the Columbia River as an adult. Fish were identified by Columbia Basin Research (CBR) as having migrated volitionally or having been transported by barge. All PIT-tag files are available on the CBR website (http://www.cbr.washington.edu/dart/cs/data/nmfs_sar/). We excluded all fish with an unknown rearing type (i.e., hatchery versus wild), geographic regions with fewer than 200 individuals (over the 16 years), those fish released or tagged below the confluence of the Snake and Columbia Rivers, and fish that returned to spawn in less than 1 year. In addition, we excluded fish that passed Bonneville Dam prior to April 9th or later than July 8th (<0.14% of the total observations), because there was insufficient data to inform the temporal autocorrelation at the margins of the migration period. In total, there were 122,415 wild individuals for the SAR analysis.

Life-cycle model structure

We employed a stochastic, age-structured model modified from previous life‑cycle models37,38,101. The model as depicted in Fig. 8 has five annual time‑steps, based on the typical maximum 5-year generation time of Snake River spring/summer Chinook salmon. These time‑steps correspond approximately to life-stage transitions.

Fig. 8: Diagram of life-cycle model. Environmental covariate survival models that were fit directly to PIT-tag data are shown in blue (S tributary , S mainstem , S SAR , S upstream ). Life stages and fitted transition parameters are shown in black, with associated equations for reference in the Supplementary Methods. Full size image

The first time‑step (“spawner to parr” stage) spans spawning (August/September), incubation, and early parr rearing. Survival through the second timestep (S 2 ) includes both tributary rearing (S tributary ) from summer (July or August) to the following spring when migrating juveniles pass Lower Granite Dam (April/May), and continue (S mainstem ) through the Snake and Columbia River hydrosystem.

The third timestep includes ocean entry (May/June) and the first winter and spring in the ocean. Some Chinook salmon return to spawn in their third year (jacks), but most females and mature adult males stay in the ocean for one or two more years, during which ocean survival is represented as S o . The number of fish in the ocean each year is a latent variable, fit by detections of survivors when they re-enter freshwater (“smolt to adult return,” or S sar ). Upstream migration survival through the hydrosystem (S upstream ) and from Lower Granite Dam to spawning areas (S prespawn ) are captured in the fourth or fifth timestep. Older females tend to lay more eggs, which is reflected in the fecundity parameter, F 5 .

To calculate “effective spawners,” we combined different age classes that returned to spawn in the same year as the weighted sum of three (weight = 0), 4- (weight = 1), and 5‑year‑old fish (weight = F 5 ). This sum was estimated by an expansion of the number of redds (nests) counted during spawning surveys45.

Statistics and reproducibility

We fit the life-cycle model in two steps. First, we fit individual life‑stage relationships with covariates using a variety of methods for different life stages to generate a posterior distribution for stage-specific parameter estimates. In addition to these parameters, the life-cycle model introduced some additional parameters that could not be directly fit to the data. We, therefore, conducted a second step to calibrate the life cycle model using the Approximate Bayesian Computing approach102,103. All analyses were conducted in RStudio (ver. 1.3.1073) using R (ver. 3.6.2)104.

In the calibration step, we simulated population time series under recent climatic conditions using the life-cycle model. We compared the resulting spawner and smolt time series to those observed and ranked parameter sets by deviance from a Kolmogorov–Smirnov test. We selected the top 0.2% of 500,000 sets of parameter combinations (Supplementary Fig. 1). All resulting parameter sets met the criterion of producing a spawner distribution that was not statistically different from the observation dataset (P value > 0.05 from Kolmogorov–Smirnov test). This process maintained the appropriate correlation structure within each parameter set.

Spawner to parr (S 1 ) and parr to smolt (S tributary ) stages

We fit adult recruits per spawner for eight populations in a hierarchical Bayesian framework using multiple likelihood equations. These equations reflected stages that could be compared directly with data. Briefly, we fit a two-stage Gompertz function105 (see Supplementary Methods) to solve the two stages simultaneously, combined with independently estimated survival rates for later stages. Individual population coefficients (productivity and capacity parameters for both stages as well as coefficients for temperature and flow) were assumed to be random samples from an underlying normal distribution106 (priors shown in Supplementary Table 4).

To determine the best environmental covariates, we compared the estimated predictive error of alternative models using a leave-one-out, cross-validation method for Bayesian models using the LOO package107. Because of correlations among the climate variables tested, multiple models had similar support from the data (Supplementary Table 5). Our primary objective was to identify divergence in potential responses to climate change. Therefore, we selected two models with covariates that showed differing trajectories with climate change (Fig. 2). Model 1 included summer air temperature and fall streamflow, while model 2 included summer streamflow only.

Juvenile survival through the Columbia River hydrosytem (S mainstem )

We estimated juvenile survival from Lower Granite to Bonneville and arrival day at Bonneville Dam using the COMPASS model108. In the COMPASS model, equations for each dam and riverine reach used measured hourly flow, temperature, and spill to predict fish survival and migration rate. Dams had additional equations to determine the proportion of fish that passed using the spillway, turbine, or bypass system as a passage route. The COMPASS model also tracked the proportion of fish loaded into barges for transport and release downstream from the hydropower system.

Smolt-to-adult return (S sar )

We used a mixed-effects logistic regression model to determine the effects of ocean entry date and environmental covariates on the probability that an individual fish would return as an adult to Bonneville Dam. The model109 included random effects of day and day‑by‑year interaction, which followed an auto‑regressive process. Using Akaike Information Criterion, we selected variables with high importance based on model weights, which included a large-scale measure of SST (SSTarc) and a more local, coastal measure of SST (SSTwa), as well as spring upwelling. We applied separate models for fish that had migrated through the mainstem river vs. those that had been transported (Supplementary Table 6).

Adult upstream survival (S upstream )

For the adult upstream survival model, we used generalized additive mixed models (GAMMs) to evaluate the effects of both anthropogenic and environmental covariates on spring/summer Chinook salmon survival47. To run the model in simulation mode, all the non-environmental covariates (fisheries catch, proportion transported as juveniles) had distributions similar to those during the baseline period of 2004–2016. Survival from the hydrosystem to spawning streams (S prespawn ), was held constant due to lack of appropriate data with which to fit a relationship.

Calibration step

The model included a set of maturation parameters that could not be estimated directly from the data. Therefore, we treated them as tuning parameters for the life‑cycle model as a whole. These parameters partitioned total smolt-to-adult survival (SAR) into age-specific survival rates (S) and propensity to return at a given age (b) as follows:

$$S_{SAR} = \frac{{b_3\times N_3 + b_4\times N_4 + N_5}}{{N_2}}$$ (1)

where b 3 denotes the proportion of all 3-year-old fish (N 3 ) that returned as jacks, b 4 is the proportion of all 4‑year‑olds (N 4 ) that returned to spawn, and N 5 is all 5-year-old fish. We calculated the number of fish in each age class as follows:

$$N_3 = S_3 \times N_2$$ (2)

$$N_4 = (1 - b_3) \times N_3 \times S_0$$ (3)

$$N_5 = (1 - b_4) \times N_4 \times S_0$$ (4)

where S 3 denotes the survival rate of all smolts (N 2 ) over the first winter in the ocean and S 0 is the survival rate over subsequent years at sea. Fish that stay in the ocean longer have additional mortality, but there is a fecundity advantage (F) for older spawners because larger fish have higher fecundity. In the model, the effective number of spawners reflects the age distribution of female spawners, which return as either 4‑ or 5‑year olds, with the latter having a fecundity advantage. A percentage of fish returned as 6‑year olds, but it was so small that we pooled them with the 5‑year olds.

We fit the tuning parameters (S 0 , b 3 , b 4 , and F 5 ) using a modified Approximate Bayesian Computing approach102,103. We applied this method by first generating a prior distribution for each parameter. The priors for parameters that range between 0 and 1 (S 0 , b 3 , b 4 ) had a beta distribution centered on the mean and included the full range used by Zabel et al.38 (Table 1). For F 5 , we used the normal distribution with mean from Kareiva et al.98 with sd = 0.1. These values were randomly combined with the freshwater parameter sets to create 500,000 combinations of all ten parameters.

We ran the life-cycle model to tune each population using each of these parameter sets in a different simulation. The life-cycle model was forced by historical meteorological conditions from 2000 to 2015 (Supplementary Table 3). For juvenile mainstem survival, we used COMPASS reconstructions of juvenile migration that incorporated actual river management practices. Note that historical river management was extremely variable over the juvenile migration stages and not comparable to the conditions we projected in climate simulations. In particular, transportation rates were higher, the spill was lower, and dam‑passage survival was lower than in the corresponding actions proposed by the U.S. Army Corps of Engineers et al.110 in their environmental impact statement. We ran the SAR model in retrospective mode, which used the fitted estimates of historical random effects and observation error. We simulated other life stages in the calibration as we would in future projections.

We ran 500,000 iterations with different parameter combinations through the historical reconstruction. We selected this number to produce at least 1000 sets of parameter values for each population that produced spawner distributions similar to those observed (Fig. 9). Parameter values from the top 0.2% of these parameter sets are shown compared with their prior distributions in Supplementary Fig. 1. The resulting models produced Kolmogorov–Smirnov P values from 0.16 to 0.76, each of which rejected the null hypothesis that modeled and observed time series came from different distributions (P < 0.05).

Fig. 9: Historical time series of spawner abundance. Heavy lines (black) show empirical estimates modeled from redd‑count data for population status reviews, using methods described in ref. 39. Broken lines (red) show median abundance from life‑cycle model (LCM) simulations. The empirical estimates account for observation error and include 95% confidence limits, shown in green, while the 95% bounds from the life-cycle model are shown in yellow. Full size image

Climate scenarios

Detrended climate simulations

In creating future climate scenarios, our aim was to maintain the statistical properties of the environmental data driving survival in various life stages. Large-scale oceanic and atmospheric drivers affect marine and freshwater environments simultaneously. To account for this, we estimated an unstructured covariance matrix for all the freshwater and marine environmental covariates using the Template Model Builder (TBM) libraries for R111. We used a multivariate state-space model

$${\mathbf{x}}_t \sim {\mathrm{MVN}}\left( {\rho ^x{\mathbf{x}}_{t - 1},{\mathbf{Q}}} \right)$$ (5)

$${\mathbf{y}}_t\sim {\mathrm{N}}\left( {{\mathbf{x}}_t,0.001} \right)$$ (6)

where \({\mathrm{x}}_{\mathrm{t}}\) was a vector of environmental processes at time t, \({\rho}^{\mathrm{x}}\) was the correlation between vectors from successive time-steps, and Q was an unstructured covariance matrix. For n covariates, the number of estimated correlation coefficients in the unstructured covariance matrix was equal to \(0.5 \times n \times \left( {n - 1} \right)\). Observation error for the observation model was fixed at 0.001, essentially implying that covariate data were measured without error.

Climate change scenarios

Our objective was to specifically explore how the relationships among various climate drivers interacted with population dynamics to shape population trajectories. To do this, we first simulated 1000 time series of 75 years of detrended environmental conditions that followed the model-fits of the covariance relationships in the historical record (Supplementary Fig. 2), although different matrices were generated for different models.

Second, we added offsets to these detrended simulations according to GCM projections. Each offset (trend) consisted of a single time series added to all 1,000 detrended simulations. Then, for each climate scenario, we input the resulting combinations of detrended-plus-offset environmental conditions into the life-cycle model as forcing factors.

For each model/population/climate scenario, we ran 1000 iterations. Climate trends were created from the median (GCM 50 ) and interquartile range across 10–80 GCM projections, depending on the covariate (Supplementary Table 1). We tracked geometric mean population abundance in 10-year intervals. We also assessed the first year (if any) in which a population in a given simulation fell below a quasi-extinction threshold of adult abundance (QET50). The QET50 was passed when the running mean of spawners, measured at the spawning stream, dropped below 50 individuals in any 4-year period45.

We explored projected climate trends from two carbon emissions scenarios, representation concentration pathway (RCP) 4.5 and 8.5112,113. For marine variables, we used time-series output from GCMs directly. For freshwater variables, we used output that had been statistically downscaled and then processed through hydrological and stream temperature models70,114,115. For each variable, 10–80 time series per emissions scenario were available (Supplementary Table 1).

We used four steps to extract the relevant trends from GCMs. First, we calculated the variable mean (e.g., the average for a given month) over a baseline period centered on 2015 (2005–2025) within each time series individually. Second, we calculated anomalies from the 2005–2025 baseline period within each time series. Third, we created a 20-year running mean of anomalies for each time series. We smoothed each of these individual time series using a 20-year running mean to reduce interannual variation that was already accounted for by the TMB model. Fourth, we calculated the 25th, 50th, and 75th quantiles across all time series for each covariate in each year from the smoothed anomalies.

We repeated this process for each climate scenario, producing the trends shown in Fig. 2. These trends represented the spread across climate models of the low, medium, and high rates of change in each covariate. The smoothed trends were nearly linear because decadal and sub-decadal variability averaged out across models for most variables. However, decadal patterns were still apparent in the smoothed trends for flow, due to its very high natural variability as well as the few GCMs driving flow projections.

Although temperature and flow are physically linked (air temperatures are input into the hydrological models, and both air temperature and flow are input into stream temperature models) we found the trends for temperature were much more linear than those for flow. Thus, there was no consistent relationship between them (i.e., the rank order of temperature projections was not strongly correlated with the rank order of flow projections). We, therefore, treated the temperature and flow offsets as independent, and retained decadal variation in flow, which was low relative to annual variation. Exploration of alternative methods indicated this choice did not influence population responses.

Thus, in summary, we added each climate trend to the raw simulation of a detrended climate produced by the covariance model. All new time series retained the autocorrelation, variance, and covariance of the historical environment with forcing from greenhouse gas emissions added.

Life-cycle model simulations

Running the model in simulation mode was relatively straightforward for stages that had been fitted on an annual timestep (S 1 and S tributary ) because they only required an annual input of environmental conditions. However, submodels for upstream and downstream survival required daily environmental inputs. We, therefore, ran a separate step to produce survival estimates from numerous representations of daily time series, which were then grouped by annual spring mean temperature and flow. For each annual timestep, we sampled randomly from the appropriate disaggregated submodel output (Supplementary Methods).

We ran the life-cycle model from 2015 to 2089, applying climate trends to the detrended climate for each environmental covariate on an annual timestep. We repeated 1000 iterations for each population, per model, and climate scenario. Our first results stemmed from a single combination of freshwater and marine covariates (Model 1: summer temperature + fall flow for freshwater, SSTarc in winter for transported fish, and SSTarc in winter + SSTwa in summer for in-river fish to show the full-time series of population response to climate scenarios (Figs. 4 and 5).

Sensitivity analyses

In our first sensitivity analysis, we compared population outcomes from Model 1 with a different freshwater covariate model (Model 2: summer flow + Model 1 marine covariates). Finally, we exchanged the marine covariates in Model 3: summer temperature + fall flow for freshwater, SSTwa in summer for transported fish, and SSTarc in spring + spring upwelling for in-river fish (Fig. 6). Thus, we explored all combinations of the top two models for freshwater and marine stages, respectively.

In our second sensitivity analysis, we applied climate trends for the ensemble mean of RCP 8.5 to one life stage, while the other life stages experienced a detrended climate (Fig. 7). We cycled through survival for parr‑to‑smolt stage (S tributary ), downstream migration (S mainstem ), smolt‑to‑adult return ratio (S sar ), and upstream spawning migration (S upstream ). In each case, we reported the extent of population decline as the ratio of geometric mean population size in 2060–2069 divided by mean abundance in 2020–2029. We ran 1000 simulations per life stage and calculated the mean change in abundance for each population.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.
[END]

[1] URL: https://www.nature.com/articles/s42003-021-01734-w
[2] URL: https://creativecommons.org/licenses/by-nc-nd/3.0/us/
   URL: https://www.propublica.org/steal-our-stories

ProPublica via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/propublica/

You are viewing proxied material from magical.fish. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.