Improving short-term recruitment forecasts for coho salmon using a spatiotemporal integrated population model

Fishery managers often rely on forecasts of future population abundance to set allowable harvest quotas or exploitation rates. While there has been substantial research devoted to identifying environmental factors that can predict recruitment for individual populations, such correlations often degrade over time, thereby limiting their utility for management. Conversely, examining multiple populations at once to detect shared, spatially structured patterns can offer insights into their recruitment dynamics that are advantageous for forecasting. Here, we develop a population dynamics model for natural origin coho salmon (Oncorhynchus kisutch) stocks in Washington State that leverages spatial and temporal autocorrelation in marine survival to improve one-yearahead forecasts of adult returns. Executed in a Bayesian hierarchical integrated modelling framework, our spatiotemporal approach incorporates multiple data types and shares information among stocks to estimate key biological parameters that are informative for forecasting. Retrospective evaluation of one-year-ahead forecast skill indicated that the spatiotemporal integrated population model (ST-IPM) outperformed existing forecasts of Washington State coho salmon returns by 25–38 % on average. Moreover, the ST-IPM estimates parameters that were previously non-identifiable for many stocks, and propagates uncertainty from multiple contributing data sources into model forecasts. Our results add to a growing body of work demonstrating the utility of spatiotemporal and integrated approaches for modelling population dynamics, and the framework developed here has broad applications to the assessment and management of coho salmon in Washington State and elsewhere throughout their range.


Introduction
A central challenge to forecasting fish population dynamics lies in anticipating environmentally-driven variation in recruitment (Cushing, 1982;Walters and Martell, 2004). While retrospective analyses can often detect relationships between environmental conditions and fish production for individual populations, such correlations are often weak and diminish over time, thereby limiting their utility for tactical management (Drinkwater and Myers, 1987;Myers, 1998;Walters and Collie, 1988, but see Hare et al., 2010;Scheuerell and Williams, 2005;Tommasi et al., 2017). Among the main reasons why such relationships may be unreliable are spurious correlations among autocorrelated time-series, non-stationarity, and the confounding effect of multiple latent processes acting on fish stocks simultaneously (Kilduff et al., 2014;Litzow et al., 2019;Mueter et al., 2002b;Wells et al., 2017). Alternatively, recruitment may be better understood by examining multiple stocks at once to detect shared, spatially structured patterns (Myers and Mertz, 1998;Peterman et al., 1998;Pyper et al., 2001). Not only does a multi-population approach reduce the risk of spurious correlations, but spatial coherence in stock dynamics integrates across the many ecosystem processes that may be jointly influencing recruitment (Walters and Martell, 2004).
While recruitment patterns can be detected from spawner-recruit residuals or survival estimates, such data are often noisy and of limited availability and scope (Myers et al., 1995). Inferences from short time-series with large measurement errors can be misleading (Clark and Bjørnstad, 2004;De Valpine and Hastings, 2002;Walters and Ludwig, 1981), and reliance on a single data type can exclude complementary information available from other sources. This is particularly disadvantageous in a multi-population context where biological quantities of interest may be informed by different types of data across stocks. Conversely, 'integrated' approaches to population modelling can use all available information to estimate parameters through a joint likelihood that captures and propagates the uncertainties in each contributing data source (Deriso et al., 1985;Fournier and Archibald, 1982;Maunder and Punt, 2013). Similar to state-space models (Valpine and Hilborn, 2005), integrated population models (IPMs) describe the data as noisy realizations of underlying biological processes which are represented as latent, unobserved states (Scheuerell et al., 2020). When structured hierarchically, IPMs can facilitate sharing information among data rich and data poor populations (Jiao et al., 2011;Punt et al., 2011) and are increasingly being used to estimate spatially structured biological dynamics (Cao et al., 2019;Grüss et al., 2019;Kristensen et al., 2014;Punt, 2019) and forecast population trajectories (Buhle et al., 2018).
In anadromous Pacific salmon (Oncorhynchus spp.), recruitment is typically defined as the number of mature adults that return from a given year-class (Ricker, 1954), which depends heavily on cohort survival during marine residency. Considered a critical period in the life cycle, conditions experienced during early marine residency are particularly influential in determining year-class strength (Beamish and Mahnken, 2001;Beamish et al., 2004; but see Ruggerone and Connors, 2015). Consequently, identifying drivers of salmonid early marine survival is a subject of substantial research interest (reviewed in Beamish, 2018;Chittenden et al., 2009;Pearcy, 1992). Such investigations have identified relationships between marine survival and environmental conditions expressed at basin (e.g. Pacific Decadal Oscillation (PDO), North Pacific Gyre Oscillation (NPGO); Di Lorenzo et al., 2008;Kilduff et al., 2015;Mantua et al., 1997), regional (e.g. coastal upwelling, sea surface temperature; Kilduff et al., 2014;Koslow et al., 2002), and local (e.g. estuaries; Mahnken et al., 1998;Teo et al., 2009) scales. As a result, recruitment can be correlated among populations at spatial scales that match those of the dominant oceanic features affecting survival (Mueter et al., 2002b(Mueter et al., , 2002a. Here, we develop an integrated population model that leverages spatial correlations in marine survival to improve state-wide forecasts for Washington natural origin coho salmon (O. kisutch) returns. Preseason forecasts for these stocks are used to determine allowable harvest rates each year such that under-forecasting can lead to foregone harvest opportunities, while over-forecasting may risk overfishing (Pacific Fishery Management Council, 2016). Most Washington coho salmon spend their first eighteen months in freshwater, after which they migrate to the ocean where the majority will spend another eighteen months before returning to spawn at age three. Historically, preseason forecasts for many of these stocks were based on sibling regressions (Peterman, 1982) which used observed returns of jacks -precocious males that mature after only six months at seato predict returns of 'adult' (age three) males and females originating from the same cohort. The rationale behind this approach is that if early marine survival is a key determinant of year-class strength, then returns of jacks that matured early but still experienced this critical period can inform the survivorship of the entire cohort. However, lack of reliable jack abundance data and weakening of the jack-to-adult relationship now limits the performance of sibling regressions for natural origin coho salmon in Washington. Consequently, many forecast methods currently in use rely instead on environmental indicators to predict marine survival (e.g. Rupp et al., 2012b;Zimmerman, 2018). Unfortunately, environmentally based forecasts have failed to predict large fluctuations in abundance in some years (e.g. Wainwright, 2021). In this study, we explore an alternative forecasting approach that relies on spatial and temporal autocorrelation in marine survival rather than sibling or environmental relationships. Executed in a Bayesian hierarchical integrated modelling framework, our approach incorporates multiple data types as well as prior information, and facilitates sharing information among populations. Retrospective evaluation of one-year-ahead forecast accuracy from 2002 to 2017 indicated that the spatiotemporal IPM (ST-IPM) outperformed existing adult return forecasts by 25-38 % on average. Our results emphasize the utility of integrated and spatiotemporal approaches for modelling population dynamics, and the framework developed here has broad applications to the assessment and management of coho salmon in Washington State and elsewhere throughout their range.

Coho salmon life history
Mature coho salmon in Washington typically migrate upriver in late summer and fall to spawn between October and December in small streams and mainstem channels of larger rivers (Ohlberger et al., 2019). The year that a cohort is spawned is referred to as its 'brood year'. Embryos produced in brood year y overwinter in the gravel and emerge as fry in the spring of year y + 1. Most juveniles then spend ~ 1 year rearing in freshwater before outmigrating to the ocean as smolts in spring of the following year (y + 2) from early April to early June. Most coho salmon spend roughly eighteen months at sea before returning to spawn in their natal habitats in the fall of year y + 3 (Bradford et al., 2000). However, some males mature and return to spawn as 'jacks' after only six months at sea. Because jack abundance data are generally of lower quality for many natural origin coho salmon stocks in Washington State and jacks comprise only a small portion of total cohort recruitment (Quinn, 2005), only age three individuals are considered in the present study.

Populations and data
Our analysis includes data from thirty-six coho salmon management units (henceforth 'stocks') throughout the Salish Sea (Puget Sound, southern Strait of Georgia, and Strait of Juan de Fuca) and Washington coast ( Fig. 1, Table S1) which are currently forecasted for management purposes. These management units are defined in the Coho Fisheries Regulation Assessment Model (FRAM;Pacific Fishery Management Council Model Evaluation Workgroup, 2008), and may represent single spawning populations, or aggregations of multiple spawning populations. While all of these stocks are of natural origin (henceforth 'wild'), unknown numbers of hatchery origin fish may also be present on the spawning grounds and counted towards the escapement.
The data types included in our integrated model are (1) quantity of habitat occupied by each population (stream length), (2) adult escapement counts, (3) harvest abundance, (4) smolt outmigration counts, and (5) coded wire tag (CWT) marine survival estimates (Table S1). The quantity of occupied habitat was obtained from SalmonScape, a Washington Department of Fish and Wildlife interactive web map of species distributions (https://apps.wdfw.wa.gov/salmonscape/) and included all habitat known or presumed to be used by coho salmon for spawning and/or rearing. Adult escapement and harvest numbers were generated by FRAM (Pacific Fishery Management Council Model Evaluation Workgroup, 2008), which obtains its escapement data from the Washington Department of Fish and Wildlife (WDFW) Salmon Conservation and Reporting Engine (SCoRE; (https://fortress.wa.gov/dfw/score/sco re/species/coho.jsp?species=Coho) before aggregating it at the management unit spatial scale at which FRAM operates. These escapement data were generally derived from redd counts or area-under-the-curve estimates of live spawners expanded to account for un-surveyed areas and times. Smolt outmigration estimates were based on smolt trapping from WDFW and tribal comanagers (e.g., Anderson et al., 2019). Finally, CWT marine survival estimates were collated as part of the Salish Sea Marine Survival Project (https://marinesurvivalproject.com/) as described in Zimmerman et al. (2015). The abundances of the stocks included in this study were highly variable, with arithmetic mean (± standard deviation) annual returns ranging from 264 ± 227 fish for Port Gamble Bay, up to 166, 630 ± 87, 520 fish for the Snohomish River, and a median average return of 9705 fish across all stocks (Table S2).

Model design
The ST-IPM specifies returns of adult (age three) coho salmon as the product of smolt production and marine survival. For stocks and years where smolt outmigration data are collected, these observations are available to inform predictions of adult returns. When these data are absent, smolt outmigration is predicted as a density-dependent function of spawner abundance (escapement) (Barrowman et al., 2003). While spawner-recruit relationships have limited utility for short-term forecasting (Walters, 1989) and watersheds become fully saturated at relatively low coho salmon spawner abundance (Bradford et al., 2000), escapement-based estimates of smolt production may provide a useful baseline for predicting variation in adult returns that is not explained by ocean mortality. Moreover, propagating uncertainty from the adult-to-smolt production relationships and data into model forecasts could be advantageous for managers wishing to consider forecast uncertainty in harvest control decisions.
Marine survival can be estimated using CWT data when available or inferred from the difference between smolt outmigration counts and adult returns for a given cohort. However, such information does not become available until after the fishing season, and thus cannot be directly used in forecasting. While CWT and smolt outmigration data can be used to estimate properties of marine survival time-series such as historical averages, autocorrelation, and recent trends that may be useful for forecasting (Winship et al., 2015), these data are only available for a handful of stocks. For stocks that lack both CWT and smolt abundance data, the influences of juvenile production versus marine survival on adult returns cannot be distinguished.
The rationale behind the modelling approach presented here is to leverage the spatially correlated marine survival patterns of coho salmon (Coronado and Hilborn, 1998;Zimmerman et al., 2015) to share information among stocks. By assuming that nearby populations experience similar patterns of ocean mortality, retrospective estimates of marine survival time-series can be reconstructed even for stocks that lack informative data. Furthermore, by introducing prior information from previous studies (Barrowman et al., 2003;Thorson et al., 2015a) and adopting a hierarchical model design (Jiao et al., 2011;Punt et al., 2011), we can further refine estimates of marine survival and smolt production, allowing improved inference on statistical properties of these time-series that may be useful for forecasting.

Smolt production
For each coho salmon stock (i), the number of smolts (R i,y ) migrating to the ocean in a given calendar year (y) was assumed to depend on the number of spawners (escapement) two calendar years earlier (S i,y− 2 ) according to a Beverton-Holt function (Beverton and Holt, 1957): where S i,y− 2 and R i,y represent state estimates (i.e. model predictions) of the true spawner and smolt abundance respectively, a i is a stock-specific productivity term describing the slope (in smolts produced per spawner) of the function at the origin, R max i is the maximum asymptotic smolt production, expressed as the number of smolts produced per kilometer of habitat, and d i is the kilometers of stream habitat for each stock (Barrowman et al., 2003). While the Ricker function is also commonly used for some Pacific salmon species, overcompensation is unlikely to occur in coho salmon (Barrowman et al., 2003;Bradford et al., 2000). Deviations from the Beverton-Holt adult-to-smolt relationship (ε i,y ) were assumed to follow a lognormal distribution (Peterman, 1981) with a mean of zero and a variance of σ 2 Ri . To facilitate sharing information and allow for correlation between α and R max , these parameters were modelled hierarchically as random variables arising from a common bivariate lognormal distribution that was shared among stocks (Buhle et al., 2018): where log(μ θ ) is the vector of lognormal hypermeans for the Beverton-Holt parameters (log(μ θ ) = log(μ a ), log(μ Rmax )) and Σ θ is the variancecovariance matrix. To improve posterior sampling efficiency and mitigate the bias that can result from estimating hierarchical models using Monte Carlo methods (Betancourt, 2016;Monnahan et al., 2017), we implemented a multivariate extension of the non-centered parameterization, with the variance-covariance matrix Σ θ decomposed into the Cholesky factor (L Ωθ ) of the correlation matrix Ω θ and a vector of error where: and: In a non-centered parameterization, z θ i is a vector of standard scaling factors for each parameter and stock (z a i , z Rmax i ) that follows a normal distribution with a mean of zero and a standard deviation of one  (Betancourt, 2016;Monnahan et al., 2017). Of the thirty-six management units included in this study, smolt outmigration data of any kind were available for only nine, and among these, the data were typically available for only a subset of years (Table S1). Consequently, we chose to use informative priors based on the posterior distributions reported by a previous hierarchical analysis of coho salmon adult-to-smolt production (Barrowman et al., 2003). Note that that Barrowman et al. (2003) assume a 1:1 sex ratio and report the posterior mean of R max in units of female smolts/km, which is doubled here to represent R max in total smolts/km: To explore prior sensitivity, we compared the resulting posteriors to those from model fits that used vague prior distributions (Fig. S1). The Cholesky factor of the correlation matrix was drawn from an LKJ prior distribution (Lewandowski et al., 2009): where η is a shape parameter that specifies the expected degree of correlation between α and R max , which we fixed at 2, representing a weakly informative prior expectation of weaker correlation between parameters (Stan Development Team, 2020). Stock-specific standard deviations of Beverton-Holt adult-to-smolt errors (σ R i ) were also modelled hierarchically among management units using a non-centered approach:

Marine survival
Marine survival was assumed to be density-independent (but see Emlen et al., 1990) and unrelated to variation in freshwater production (but see Chasco et al., 2021;Haeseker et al., 2012;McCormick et al., 2009). While there is evidence of density-dependence in the marine phase (e.g. Ruggerone and Connors, 2015), this is generally observed in more abundant species such as sockeye (O. nerka) and pink (O. gorbuscha) salmon, and operates based on the aggregate density of salmon in ocean foraging areas rather than the abundance of any one stock (Ohlberger et al., 2019;Pyper and Peterman, 1999;Ruggerone and Nielsen, 2004).
Information on a population's marine survival can come from CWT data or be inferred from the difference between smolt outmigration and adult returns. Of the thirty-six stocks included in our analysis, CWT data were available for only fifteen, and smolt outmigration data were available for an additional two (excluding the seven stocks for which both data exist) (Table S1). Because marine survival of southern coho salmon populations can be spatially correlated at relatively fine spatial scales (Zimmerman et al., 2015), we specified marine survival anomalies as a spatial Gaussian field (a Gaussian process in two or more dimensions, e.g. Ward et al., 2015;Webster et al., 2020) to facilitate sharing information among stocks. For each stock i, marine survival over time was expressed as a mean-reverting lag-1 autoregressive (AR-1) process of the form: where λ i,y is the marine survival in year y for stock i and μ λi is the mean of the logit marine survival time series for stock i, which followed a hierarchical normal distribution among stocks with hyperparameters μ μ λ and σ μ λ . In model year-one, there is no previous state estimate to inform that year's marine survival so the time-series for each stock was initialized at ψ i , which was hierarchically normally distributed among stocks with a mean of μ ψ and a variance of σ ψ 2 . The autocorrelation terms for each stock (ϕ i ) were drawn from uniform prior distributions bounded between -1 and 1. The marine survival deviations for each stock (ξ i,y ) were multivariate normally distributed with a mean of zero and variancecovariance matrix of Σ ξ . The covariance between stocks i and j was expressed as a function of Euclidean distance between their marine entry locations according to a squared exponential kernel of the form: where, x i and x j are the coordinates of the marine entry point for stocks i and j in eastings and northings, γ 2 is the marginal variance of the function, ρ is the length scale, and σ d is the error standard deviation, which is applied only when i = j according to the Kronecker delta function δ i,j .

Harvest and escapement
Of the individuals that survive natural marine mortality (R i,y λ i,y ), a portion are harvested in fisheries: where C i,y is the number of individuals from stock i harvested in year y according to the exploitation rate u i,y , which was specified as a multivariate random walk: The logit harvest rates for model year-one (ϑ i ) were initialized hierarchically as in eq. 12 with hyper-parameters μ ϑ and σ ϑ . Harvest process errors (ε i,y ) were multivariate normally distributed with a mean of 0 and variance-covariance matrix Σ ε , which was parameterized with a single variance term (σ 2 ε ) on the diagonal, and covariance (ρ ε σ 2 ε ) on the off-diagonal elements (e.g. Holmes et al., 2012). The number of individuals returning to spawn in a given year (S i,y ) was then calculated as the difference between the total number of smolts that survived natural ocean mortality in the previous year (R i,y− 1 λ i,y− 1 ) minus those that were harvested (C i,y ).
State estimates of spawning abundance were then used recursively in the subsequent estimation of smolt production (eq. 1).

Likelihoods
The observed smolt (J i,y ), escapement (E i,y ), and harvest (H i,y ) abundance data were assumed to follow lognormal likelihoods: where σ J , σ E and σ H represent the observation error terms for the smolt, escapement and harvest data respectively, and R i,y , S i,y , and C i,y are the model-generated state estimates of smolt, escapement, and harvest abundance. The escapement observation error term (σ E ) could not be reliably estimated, and so was fixed at 0.2 and subject to testing with alternative values (e.g. Fleischman et al., 2013). For the CWT data, the estimated number of tagged fish that were recovered (n i,y ) was assumed to follow a beta-binomial likelihood with a number of trials (N i,y ) equal to the total number of tagged fish that were released, and a probability of recovery (p i ,y ).
To allow for extra-binomial variance in tag recoveries (for instance, due to incomplete sampling of the harvest and escapement), the probability of recovery was assumed to follow a beta distribution, implemented as a conjugate prior to the binomial: The shape parameters of the beta distribution (α i,y ,β i,y ) were specified in terms of mode (λ') and concentration (κ): where λ ' i,y corresponds to the adjusted state estimate of marine survival. For ten of the stocks in our study, available CWT data were collected from an adjacent hatchery rather than the wild stock itself. Hatchery and wild coho salmon populations have been shown to exhibit similar trends and interannual patterns in marine survival, but average mortality is typically lower for wild stocks (Coronado and Hilborn, 1998). As such, a hierarchically distributed offset term (e.g. Ohlberger et al., 2019) was applied to the likelihood to account for hatchery-specific deviations in survival from associated wild populations: Here, τ i is the marine survival offset term for stock i which is applied if the marine survival data for that stock is based on a hatchery ( h i = 1). For stocks with marine survival data based on CWT recoveries of wild fish ( h i = 0) no adjustment is required. The hatchery offset terms themselves were hierarchically distributed: where μ τ is the average offset term among stocks and σ τ is the standard deviation.

Model estimation and validation
Posterior sampling was performed via Hamiltonian Monte Carlo (HMC) No-U-turn sampling (NUTS) through the Stan model building software (Stan Development Team, 2020), implemented in R (R Core team, 2015) via the Rstan package (Gelman, 2014). Sampling occurred using five HMC chains with lengths of 2000 iterations (for simulated forecast trials) to 20,000 iterations (for fits to the complete data set). The first half of samples was discarded as a "warmup" and each subsequent sample was saved to build the posterior distribution. Convergence was assessed using the Gelman-Rubin diagnostic (Gelman & Rubin, 1992) and effective number of samples, as well as trace-plots and autocorrelation plots of HMC chains. Posterior sampling was monitored for divergent transitions and low Bayesian Fraction of Missing Information (BFMI), neither of which were present in fits to the complete data set. Model goodness of fit was assessed by examining model fits to the observed data (Fig. S2) and comparing the model's predictive distributions to observed data (posterior predictive check, Fig. S3). Model performance was assessed by simulating data with known parameter values and evaluating the model's ability to recover them (Fig. S4-S6). Prior influence was examined by comparing prior and posterior distributions for model parameters (Fig. S7). A complete glossary of all model terms and prior distributions can be found in Table S3.

Forecast evaluation
The ST-IPM generates one-year-ahead adult return forecasts for a given calendar year y by multiplying the preceding year's marine survival (λ i,y-1 ) and smolt outmigration estimates (R i,y-1 ) for each stock. In simulated forecast trials, we assume that model fits are informed only by data that would realistically be available to biologists and managers at the time of forecast development. According to the current management cycle, preseason forecasts of adult returns for the fall of year y are prepared in January of year y, when the most recent smolt outmigration data available are from year y-1. However, while harvest data from fall of year y-1 will have been collected by January of year y, these data are seldom formalized and entered into widely accessible databases at this point in time. As such, we assume that only harvest data up to year y-2 are available for conditioning model forecasts of year y. Similarly, only CWT recoveries from outmigrating smolts tagged in year y-3 are assumed to be available for forecasts of year y due to the constraints of collecting and processing these data. The time at which the escapement data become available is variable among management units, but we assume that at least preliminary escapement estimates for year y-1 would be available in January of year y to inform forecast development.
To provide both sufficient training data to condition the model and enough simulated forecasts to reliably calculate performance metrics, we produced one-year-ahead forecasts of adult returns for 2002− 2017. Forecast accuracy was calculated using the arithmetic mean absolute scaled error (MASE): Where P i,y is the forecasted adult return (sum of model-estimated harvest (C i,y ) and escapement (E i,y )) for stock i in year y, A i,y is the observed return (the sum of the observed harvest and escapement), and 1 n− 1 ∑ n 1 ⃒ ⃒ A i,t − A i,t− 1 | measures the degree of interannual variability between training years (t) in the observed adult returns during the training period of length n years. MASE has a number of advantages over alternative accuracy metrics (e.g. Mean absolute predictive error (MAPE), Root mean square error (RMSE)), including scale independence, symmetry, insensitivity to outliers, and interpretability (Hyndman and Koehler, 2006;Ward et al., 2014). We also compared forecast accuracy using RMSE, and Median Symmetric Accuracy (MSA) of the Log Accuracy Ratio (LAR) (Morley et al., 2018), which did not qualitatively alter our findings (Table S4). For all metrics, forecast skill was evaluated with respect to both the observed adult abundance data, as well as state estimates produced by the ST-IPM conditioned on all years' data (Table S4, Fig. S9-S10). We compared the forecast skill of the ST-IPM to published records of past forecasts for these stocks that were agreed upon by the state and tribal co-managers. The methods used to generate these forecasts vary among management units and over time, and the specific details of each approach are not necessarily publicly documented. In addition to the published forecasts, we also compared the performance of the ST-IPM to state-space implementations of several common univariate time-series models fitted to the adult return data, including a random walk and lag -1 autoregressive (AR-1) and moving average (MA-1) models.

Smolt production
There was substantial variability among coho salmon stocks in both smolt productivity (α) and capacity (R max ) (Figs. 2 and 3). Median smolt capacity among stocks (μ Rmax ) was estimated to have a median value of 687 smolts/km (95% credible interval = 487 to 960 smolts/km), with a lognormal standard deviation (σ Rmax ) of 0.76 (95% credible interval = 0.57 to 1.0). Note that we do not apply a bias correction factor to lognormal hypermean parameters (i.e. log(μ Rmax ) log(μ a ), log(μ σR )) such that their exponentiated values (μ Rmax , μ a and μ σR ) are interpretable as medians of the among-stock distributions. Our μ Rmax values differ from those reported by Barrowman et al. (2003) (μ Rmax = 1437 smolts/km, σ Rmax = 0.64), indicating lower median asymptotic smolt capacity per km, and greater variability among the stocks in our study. Area 12− 12B, the Puyallup and Samish Rivers, Area 12C-12D, and the Snohomish River exhibited the largest R max estimates, while Port Gamble Bay, Area 7− 7A, Grays Harbor, and the Dungeness and Green Rivers had the lowest (Figs. 2,3). Median smolt productivity among stocks (μ a ) was estimated to be 342 smolts/spawner (95% credible interval = 246 to 764 smolts/spawner) with an among-stock standard deviation (σ a ) of 0.84 (95% credible interval = 0.6 to 1.15). Despite the use of informative priors, these values differ substantially from those reported by Barrowman et al. (2003) (μ a = 71.52 smolts/spawner , σ a = 0.43), indicating greater median smolt productivity, and greater variability among the stocks in our study. Area 12A, the Nooksack, Humptulips, and Dungeness Rivers, and Port Gamble Bay exhibited some of the largest estimates of smolt productivity, while the Baker and Samish Rivers, East Juan de Fuca, Puyallup and Deschutes stocks had the lowest (Figs. 2,3). The among-population median of the Beverton-Holt adult-to-smolt error terms (μ σR ) was 0.36 (95% credible interval = 0.25 to 0.47), with a lognormal standard deviation (σ σR ) of 0.78 (95 % credible interval = 0.54 to 1.15). Areas 12A and 13A, Lake Washington, the Skokomish River, and Area 13B were estimated to have the greatest variation about the adult-to-smolt relationship, while the Queets, Skagit, Quillayute, Snohomish, and Hoh Rivers showed the least (Figs. 2,3)  (Fig. 4).

Marine survival
Average marine survival across all stocks and years (Logit − 1 (μ μ λ )) was estimated to be 0.051 (95% credible interval = 0.033 to 0.076), with an among-stock standard deviation of 0.55 (95% credible interval = 0.3 to 0.84). Most stocks exhibited strong lag-1 temporal autocorrelation in marine survival (ϕ i ) with only 3 stocks exhibiting a median estimate of ϕ i below 0.25, and 27 stocks with estimates greater than 0.5. Estimation of the Gaussian field indicated spatial autocorrelation in marine survival anomalies with a length scale (ρ) of 78.5 km (95 % credible interval = 63.4 to 98.1 km) (Fig. 4). This value of ρ implies a median correlation of 0.5 at a Euclidean distance of ~ 93 km, and a correlation of less than 0.1 between stocks whose marine entry points were ~168 km or more apart (Fig. 4). Clusters of stocks with highly correlated marine survival corresponded geographically to the Strait of Juan de Fuca, northern Puget Sound and the Strait of Georgia, central and southern Puget Sound, and the Washington coast (Figs. 1,5 ).

Forecast performance
Evaluation of one-year-ahead forecast skill from 2002 to 2017 averaged across the thirty-six stocks considered in this study indicated that the MASE for the existing published forecasts was 0.93 while that of the ST-IPM was 0.70, a difference of 25 % (Figs. 6,7). The random walk, AR-1, and MA-1 models exhibited MASE values of 0.74, 0.72 0.76, and respectively (Table S4, Fig. S8). Importantly, the ST-IPM did not always outperform the published forecasts, exhibiting slightly greater MASE values for the Stillaguamish, Skokomish, Hoh, Elwha River, East Juan de Fuca, Area 12C-12D, Area 11, Area 10E, and Area 10 stocks. The ST-IPM and published forecasts also differed somewhat in their tendency to produce large (>|50|%) biases. The existing forecast methods overforecasted by 50% or more in 187 instances (stock-year combinations) and under-forecasted by 50% or more in 139. Conversely, the ST-IPM over-forecasted by 50% or more in 205 instances and under-forecasted by 50 % in 86 instances. The published forecasts have been developed using different approaches over time, and it is possible that aggregating over sixteen years could mask recent improvements in methodology. However, restricting our comparison of forecast skill to only the last five years in our study (2013-2017) did not qualitatively alter the relative performance of each method, resulting in a MASE of 0.83 for the published forecasts versus 0.61 for ST-IPM. Using alternative metrics of forecast skill had little qualitative effect on the overall assessment of forecast performance (Table S4). Similarly, evaluation of forecast accuracy with respect to the state estimates of adult returns instead of observed data did not qualitatively alter our findings (Table S4, Fig. S9-S10). MASE for the published forecasts with respect to state estimates of adult returns was 1.1, compared to 0.79 for the ST-IPM, and 0.83, 0.81, and 0.86 for the random walk, AR-1, and MA-1 models respectively (Table S4, Fig. S9-S10).

Discussion
We developed a spatiotemporal integrated population model (ST-IPM) to forecast adult returns of wild Washington State coho salmon. In retrospective evaluations of one-year-ahead forecast skill with respect to both state estimates and observed data, the ST-IPM outperformed the existing published forecasts for these stocks by ~25-38% on average, depending on the specific metrics used. There are several features of the ST-IPM that likely contribute to its forecast skill. A hierarchical and integrated design allows the ST-IPM to incorporate multiple data types and share information among stocks (Buhle et al., 2018;Jiao et al., 2011) while stage-specific modelling of the life cycle disentangles the effects of juvenile production from ocean mortality (Rose, 2000;  Scheuerell et al., 2020). Additionally, leveraging spatial correlations among populations can improve estimation of shared, environmentally-driven processes, particularly for those that lack informative data (Thorson et al., 2013). Collectively, these characteristics of the ST-IPM facilitate estimation of key biological quantities such as marine survival and temporal autocorrelation therein, which we found to be substantial (>0.5) for most stocks in our analysis.
Temporal autocorrelation can be a particularly useful property for predicting future states (Johnson et al., 2016;Punt, 2011;Winship et al., 2015). By assuming that conditions in the near future will be similar to recent observations, forecasting via autocorrelation represents an implicit treatment of environmental effects on population dynamics (Haltuch et al., 2018). Implicit approaches to environment-recruitment modelling avoid the challenges of identifying explicit functional relationships between environmental variables and stock dynamics, and are generally less prone to prediction error (Johnson et al., 2016;Punt, 2011;Winship et al., 2015). For Pacific salmon in particular, ecosystem impacts on ocean survival can be the result of many interacting factors that have indirect, nonlinear, or cascading effects on cohorts during early marine residency (Emmett et al., 2006;Schroeder et al., 2014;Tucker et al., 2016;Wells et al., 2017Wells et al., , 2016. Such complex environmental dynamics will be difficult to predict using explicit mechanistic models, but may manifest as spatiotemporal autocorrelation in affected biological processes such as growth and survival (Mueter et al., 2002a;Mueter et al., 2002b;Peterman and Dorner, 2012;Pyper et al., 2005).
While our analysis indicates that the ST-IPM outperforms the existing published forecasts on average, there are some stocks for which other forecast approaches appear better suited. The degree to which any alternative method outperformed the ST-IPM for a given stock was generally minor, but nonetheless it may be beneficial for managers to compare among models for each stock, or consider an ensemble approach (Jardim et al., 2020;Stewart and Hicks, 2018). Although we found that the ST-IPM generally outperformed the random walk, moving average (MA-1), and autoregressive (AR-1) models (Fig. S8, S10), the improvements were often minor and simpler models offer advantages in ease of implementation and transparency to stakeholders that may offset Fig. 3. Estimated adult-to-smolt production relationships for coho salmon stocks. In each panel, the median deterministic portion of the Berverton-Holt adult-tosmolt function is shown as a solid blue line, with the 50% and 95% credible intervals shown as dark and light shaded blue boundaries respectively. Where available, observed adult-to-smolt data are plotted as dark grey circles. Individual adult-to-smolt state estimates are shown in all panels as blue circles.
a small loss of forecast skill. Previous research has shown that simple autoregressive forecasts outperform more complex models (Ward et al., 2014), so it is not necessarily surprising that the AR-1 exhibited comparable forecast skill to the ST-IPM despite the greater complexity of the latter. Given the large uncertainty in the adult-to-smolt relationships for many stocks, this component of the ST-IPM may contribute little to, or possibly detract from forecast accuracy. However, for managers interested in considering uncertainty in harvest control decisions (e.g. Privitera- Johnson and Punt, 2020), propagating uncertainty from the adult-to-smolt production relationships into estimates of forecast uncertainty may nonetheless be desirable. Furthermore, the integrated life cycle design of the ST-IPM offers additional functionality (described below) that may be valuable to managers beyond its use as a forecast model.
It is important to note that even in cases where the ST-IPM does improve forecast skill, it may not necessarily lead to better management outcomes. Closed loop simulation analysis of sockeye salmon (O. nerka) stocks has demonstrated that implementation error can negate the potential benefits of improved preseason forecasts (Dorner et al., 2009). However, Walters (1989) found that the value of preseason forecasts for Pacific salmon populations can be substantial when opportunities for in-season fisheries management are limited, as is the case for Washington coho stocks. Understanding the relative importance of preseason forecast accuracy versus other factors (e.g. harvest control rules, environmental conditions) to achieving management objectives should be an important consideration in identifying research and management priorities (Rupp et al., 2012a). It is worth noting that forecast error for several stocks remained substantial under every modelling approach considered here, and management outcomes for these stocks may be more tractably improved by developing in-season management capacity, or adopting harvest control rules that explicitly account for uncertainty (Privitera- Johnson and Punt, 2020). We recommend that future studies evaluate management outcomes across a suite of harvest control rules, forecast methods, and environmental scenarios using closed-loop simulation analysis.
There are several next steps that could be taken to continue development of the ST-IPM. The design of the spatial Gaussian field could be improved by specifying covariance as a function of marine 'over-water' rather than Euclidean distance (e.g. Hocking et al., 2018) and allowing for anisotropic covariance, such that decorrelation distance varies depending on direction (Thorson et al., 2015b). The ST-IPM's current treatment of observation error is also incomplete, as there are errors in the stock assignments of harvested fish that we did not explicitly consider. Greater transparency in the data inputs and assumptions used to generate FRAM harvest estimates will be necessary to appropriately propagate these uncertainties. Furthermore, the stocks included in our analysis likely differ substantially in the precision and bias of their smolt outmigration and adult escapement counts, which the ST-IPM does not account for. The ST-IPM's estimates of adult-to-smolt production may also be biased by unknown levels of hatchery-origin spawners present in the wild escapement (Falcy and Suring, 2018). Future model developments could be made to estimate the prevalence of hatchery-origin fish within the spawning population (e.g. Buhle et al., 2018), although the data available to do so may be limited for many of the stocks included in this study.
The presence of hatchery-origin spawners, errors in FRAM stock assignments, and variable data quality may have contributed to our estimates of smolt productivity and capacity differing substantially from those of Barrowman et al. (2003). Given that Barrowman et al. (2003) focused on populations with little or no hatchery-origin spawners and high quality smolt and adult enumeration data, we recommend that our estimates (i.e. μ Rmax , μ α ) be interpreted more cautiously against theirs from a biological standpoint. The presence of hatchery-origin spawners misattributed as recruits would tend to positively bias estimates of productivity by inflating recruitment resulting from low parental spawner abundances. Systematic misattribution of harvest among populations due to errors in FRAM fishery stock assignments could have similar effects. Furthermore, our estimates of Beverton-Holt parameters and hyperparameters were often highly uncertain (Fig. 2) and sensitive to the priors that were assumed (compare Fig. 2 to Fig. S1). Future efforts to address data quality issues and differentiate hatchery and wild  spawners could potentially resolve some of the differences between our estimates of productivity and capacity and those reported by Barrowman et al. (2003). However, there are also plausible biological reasons for the discrepancies, such as the fact that Barrowman et al. (2003) generally included smaller streams with high quality rearing habitat, potentially explaining the generally lower values of smolt capacity estimated in our study.
While the ST-IPM does not currently incorporate environmental information (other than freshwater habitat size) the model structure can readily accommodate covariates in both the adult-to-smolt and marine survival components (Maunder and Thorson, 2019;Maunder and Watters, 2003;Miller et al., 2016;Schirripa et al., 2009;Subbey et al., 2014). Freshwater covariates such as river discharge (Lawson et al., 2004;Mathews and Olson, 1980;Ohlberger et al., 2018) and habitat quality (Sharma and Hilborn, 2001) may explain some variation in adult-to-smolt production and offer predictive power for stocks that lack smolt outmigration observations. Similarly, ocean indicators could be evaluated for their ability to 'soak up' autocorrelation in marine mortality or explain independent residual variation. Failing to account for temporal and spatial autocorrelation can impede detection of robust correlations between environmental conditions and biological processes (Dormann, 2007;Walters and Martell, 2004) such that the ST-IPM may serve as a useful framework for identifying and evaluating such relationships. While beyond the scope of the present study, future work could use the ST-IPM to evaluate plausible environmental variables and functional forms (e.g. linear, nonlinear, non-stationary etc.). Inclusion of marine covariates could improve forecast performance (Logerwell et al., 2003), or may simply be informative towards process-level understanding of salmonid marine survival (Beamish et al., 2000;Quinn et al., 2005;Sharma et al., 2013;Zimmerman et al., 2015). While the ST-IPM was used for short-term forecasting here, there are several other applications in which it may also be well-suited. As a life cycle model, the ST-IPM generates estimates of stage-specific productivity and capacity, as well as the spawner abundance that produces maximum sustained yield (S MSY ) (Moussalli and Hilborn, 1986;Ohlberger et al., 2019). As such, the ST-IPMparticularly if coupled with FRAMcould lead to an improved stock assessment framework for coho salmon that uses all available data, shares information among stocks, and propagates uncertainties into estimates of stock status and management reference points. The ST-IPM also produces state estimates of marine survival (adjusted for hatchery and sampling biases) for many wild coho salmon stocks throughout Washington State that lack coded wire tag data. Such estimates are not only useful for understanding and projecting the dynamics of these stocks (Buhle et al., 2018), but may also be instructive in continuing investigations of Salish Sea marine survival. Finally, by providing a cohesive structural model of coho salmon population dynamics with parameter estimates and associated uncertainties, the ST-IPM can readily serve as an operating model for management strategy evaluation (MSE, Punt et al., 2016). Such analyses could investigate the performance of alternative harvest control policies under a range of environmental scenarios, consider the impacts of data quality and availability on management performance, or explore effects of habitat alterations and interventions at various stages of the life cycle. It would be useful for future studies to pursue these applications of the ST-IPM while continuing its development and operationalization as a forecasting tool.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.