Marine heatwaves are not a dominant driver of change in demersal fishes.

Marine heatwaves have been linked to negative ecological effects in recent decades1,2. If marine heatwaves regularly induce community reorganization and biomass collapses in fishes, the consequences could be catastrophic for ecosystems, fisheries and human communities3,4. However, the extent to which marine heatwaves have negative impacts on fish biomass or community composition, or even whether their effects can be distinguished from natural and sampling variability, remains unclear. We investigated the effects of 248 sea-bottom heatwaves from 1993 to 2019 on marine fishes by analysing 82,322 hauls (samples) from long-term scientific surveys of continental shelf ecosystems in North America and Europe spanning the subtropics to the Arctic. Here we show that the effects of marine heatwaves on fish biomass were often minimal and could not be distinguished from natural and sampling variability. Furthermore, marine heatwaves were not consistently associated with tropicalization (gain of warm-affiliated species) or deborealization (loss of cold-affiliated species) in these ecosystems. Although steep declines in biomass occasionally occurred after marine heatwaves, these were the exception, not the rule. Against the highly variable backdrop of ocean ecosystems, marine heatwaves have not driven biomass change or community turnover in fish communities that support many of the world's largest and most productive fisheries.

phylosem: A fast and simple R package for phylogenetic inference and trait imputation using phylogenetic structural equation models.

Phylogenetic comparative methods (PCMs) can be used to study evolutionary relationships and trade-offs among species traits. Analysts using PCM may want to (1) include latent variables, (2) estimate complex trait interdependencies, (3) predict missing trait values, (4) condition predicted traits upon phylogenetic correlations and (5) estimate relationships as slope parameters that can be compared with alternative regression methods. The Comprehensive R Archive Network (CRAN) includes well-documented software for phylogenetic linear models (phylolm), phylogenetic path analysis (phylopath), phylogenetic trait imputation (Rphylopars) and structural equation models (sem), but none of these can simultaneously accomplish all five analytical goals. We therefore introduce a new package phylosem for phylogenetic structural equation models (PSEM) and summarize features and interface. We also describe new analytical options, where users can specify any combination of Ornstein-Uhlenbeck, Pagel's-δ and Pagel's-λ transformations for species covariance. For the first time, we show that PSEM exactly reproduces estimates (and standard errors) for simplified cases that are feasible in sem, phylopath, phylolm and Rphylopars and demonstrate the approach by replicating a well-known case study involving trade-offs in plant energy budgets.

Range edges of North American marine species are tracking temperature over decades.

Understanding the dynamics of species range edges in the modern era is key to addressing fundamental biogeographic questions about abiotic and biotic drivers of species distributions. Range edges are where colonization and extirpation processes unfold, and so these dynamics are also important to understand for effective natural resource management and conservation. However, few studies to date have analyzed time series of range edge positions in the context of climate change, in part because range edges are difficult to detect. We first quantified positions for 165 range edges of marine fishes and invertebrates from three U.S. continental shelf regions using up to five decades of survey data and a spatiotemporal model to account for sampling and measurement variability. We then analyzed whether those range edges maintained their edge thermal niche-the temperatures found at the range edge position-over time. A large majority of range edges (88%) maintained either summer or winter temperature extremes at the range edge over the study period, and most maintained both (76%), although not all of those range edges shifted in space. However, we also found numerous range edges-particularly poleward edges and edges in the region that experienced the most warming-that did not shift at all, shifted further than predicted by temperature alone, or shifted opposite the direction expected, underscoring the multiplicity of factors that drive changes in range edge positions. This study suggests that range edges of temperate marine species have largely maintained the same edge thermal niche during periods of rapid change and provides a blueprint for testing whether and to what degree species range edges track temperature in general.

Empirical orthogonal function regression: Linking population biology to spatial varying environmental conditions using climate projections.

Ecologists and oceanographers inform population and ecosystem management by identifying the physical drivers of ecological dynamics. However, different research communities use different analytical tools where, for example, physical oceanographers often apply rank-reduction techniques (a.k.a. empirical orthogonal functions [EOF]) to identify indicators that represent dominant modes of physical variability, whereas population ecologists use dynamical models that incorporate physical indicators as covariates. Simultaneously modeling physical and biological processes would have several benefits, including improved communication across sub-fields; more efficient use of limited data; and the ability to compare importance of physical and biological drivers for population dynamics. Here, we develop a new statistical technique, EOF regression, which jointly models population-scale dynamics and spatially distributed physical dynamics. EOF regression is fitted using maximum-likelihood techniques and applies a generalized EOF analysis to environmental measurements, estimates one or more time series representing modes of environmental variability, and simultaneously estimates the association of this time series with biological measurements. By doing so, it identifies a spatial map of environmental conditions that are best correlated with annual variability in the biological process. We demonstrate this method using a linear (Ricker) model for early-life survival ("recruitment") of three groundfish species in the eastern Bering Sea from 1982 to 2016, combined with measurements and end-of-century projections for bottom and sea surface temperature. Results suggest that (a) we can forecast biological dynamics while applying delta-correction and statistical downscaling to calibrate measurements and projected physical variables, (b) physical drivers are statistically significant for Pacific cod and walleye pollock recruitment, (c) separately analyzing physical and biological variables fails to identify the significant association for walleye pollock, and (d) cod and pollock will likely have reduced recruitment given forecasted temperatures over future decades.

Spatial heterogeneity contributes more to portfolio effects than species variability in bottom-associated marine fishes.

Variance of community abundance will be reduced relative to its theoretical maximum whenever population densities fluctuate asynchronously. Fishing communities and mobile predators can switch among fish species and/or fishing locations with asynchronous dynamics, thereby buffering against variable resource densities (termed 'portfolio effects', PEs). However, whether variation among species or locations represent the dominant contributor to PE remains relatively unexplored. Here, we apply a spatio-temporal model to multidecadal time series (1982-2015) for 20 bottom-associated fishes in seven marine ecosystems. For each ecosystem, we compute the reduction in variance over time in total biomass relative to its theoretical maximum if species and locations were perfectly correlated (total PE). We also compute the reduction in variance due to asynchrony among species at each location (species PE) or the reduction due to asynchrony among locations for each species (spatial PE). We specifically compute total, species and spatial PE in 10-year moving windows to detect changes over time. Our analyses revealed that spatial PE are stronger than species PE in six of seven ecosystems, and that ecosystems where species PE is constant over time can exhibit shifts in locations that strongly contribute to PE. We therefore recommend that spatial and total PE be monitored as ecosystem indicators representing risk exposure for human and natural consumers.

A geostatistical state-space model of animal densities for stream networks.

Population dynamics are often correlated in space and time due to correlations in environmental drivers as well as synchrony induced by individual dispersal. Many statistical analyses of populations ignore potential autocorrelations and assume that survey methods (distance and time between samples) eliminate these correlations, allowing samples to be treated independently. If these assumptions are incorrect, results and therefore inference may be biased and uncertainty underestimated. We developed a novel statistical method to account for spatiotemporal correlations within dendritic stream networks, while accounting for imperfect detection in the surveys. Through simulations, we found this model decreased predictive error relative to standard statistical methods when data were spatially correlated based on stream distance and performed similarly when data were not correlated. We found that increasing the number of years surveyed substantially improved the model accuracy when estimating spatial and temporal correlation coefficients, especially from 10 to 15 yr. Increasing the number of survey sites within the network improved the performance of the nonspatial model but only marginally improved the density estimates in the spatiotemporal model. We applied this model to brook trout data from the West Susquehanna Watershed in Pennsylvania collected over 34 yr from 1981 to 2014. We found the model including temporal and spatiotemporal autocorrelation best described young of the year (YOY) and adult density patterns. YOY densities were positively related to forest cover and negatively related to spring temperatures with low temporal autocorrelation and moderately high spatiotemporal correlation. Adult densities were less strongly affected by climatic conditions and less temporally variable than YOY but with similar spatiotemporal correlation and higher temporal autocorrelation.

Estimating partial regulation in spatiotemporal models of community dynamics.

Niche-based approaches to community analysis often involve estimating a matrix of pairwise interactions among species (the "community matrix"), but this task becomes infeasible using observational data as the number of modeled species increases. As an alternative, neutral theories achieve parsimony by assuming that species within a trophic level are exchangeable, but generally cannot incorporate stabilizing interactions even when they are evident in field data. Finally, both regulated (niche) and unregulated (neutral) approaches have rarely been fitted directly to survey data using spatiotemporal statistical methods. We therefore propose a spatiotemporal and model-based approach to estimate community dynamics that are partially regulated. Specifically, we start with a neutral spatiotemporal model where all species follow ecological drift, which precludes estimating pairwise interactions. We then add regulatory relations until model selection favors stopping, where the "rank" of the interaction matrix may range from zero to the number of species. A simulation experiment shows that model selection can accurately identify the rank of the interaction matrix, and that the identified spatiotemporal model can estimate the magnitude of species interactions. A 40-yr case study for the Gulf of St. Lawrence marine community shows that recovering grey seals have an unregulated and negative relationship with demersal fishes. We therefore conclude that partial regulation is a plausible approximation to community dynamics using field data and hypothesize that estimating partial regulation will be expedient in future analyses of spatiotemporal community dynamics given limited field data. We conclude by recommending ongoing research to add explicit models for movement, so that meta-community theory can be confronted with data in a spatiotemporal statistical framework.

Within-lake habitat heterogeneity mediates community response to warming trends.

Climate change is rapidly altering many aquatic systems, and life history traits and physiological diversity create differences in organism responses. In addition, habitat diversity may be expressed on small spatial scales, and it is therefore necessary to account for variation among both species and locations when evaluating climate impacts on biological communities. Here, we investigated the effects of temperature and spatial heterogeneity on long-term community composition in a large boreal lake. We used a five-decade time series of water temperature and relative abundance of fish species captured in the littoral zone throughout the summer at 10 discrete locations around the lake. We applied a spatial dynamic factor analysis (SDFA) model to this time series, which estimates the sensitivity of each species to changing water temperature while accounting for spatiotemporal variation. This analysis described the trend in community composition at each sampling location in the lake, given their different trends in temperature over time. The SDFA indicated different magnitude and direction of species responses to temperature; some species increased while others decreased in abundance. The model also identified five unique trends in species abundance across sites and time, indicating residual dynamics in abundance after accounting for temperature effects. Thus, different regions in the lake have experienced different trajectories in community change associated with different rates of temperature change. These results highlight the importance of considering habitat heterogeneity in explaining and predicting future species abundances, and our model provides a means of visualizing spatially-explicit temporal variation in species' dynamics.

Integrating count and detection-nondetection data to model population dynamics.

There is increasing need for methods that integrate multiple data types into a single analytical framework as the spatial and temporal scale of ecological research expands. Current work on this topic primarily focuses on combining capture-recapture data from marked individuals with other data types into integrated population models. Yet, studies of species distributions and trends often rely on data from unmarked individuals across broad scales where local abundance and environmental variables may vary. We present a modeling framework for integrating detection-nondetection and count data into a single analysis to estimate population dynamics, abundance, and individual detection probabilities during sampling. Our dynamic population model assumes that site-specific abundance can change over time according to survival of individuals and gains through reproduction and immigration. The observation process for each data type is modeled by assuming that every individual present at a site has an equal probability of being detected during sampling processes. We examine our modeling approach through a series of simulations illustrating the relative value of count vs. detection-nondetection data under a variety of parameter values and survey configurations. We also provide an empirical example of the model by combining long-term detection-nondetection data (1995-2014) with newly collected count data (2015-2016) from a growing population of Barred Owl (Strix varia) in the Pacific Northwest to examine the factors influencing population abundance over time. Our model provides a foundation for incorporating unmarked data within a single framework, even in cases where sampling processes yield different detection probabilities. This approach will be useful for survey design and to researchers interested in incorporating historical or citizen science data into analyses focused on understanding how demographic rates drive population abundance.

Predicting life history parameters for all fishes worldwide.

Scientists and resource managers need to know life history parameters (e.g., average mortality rate, individual growth rate, maximum length or mass, and timing of maturity) to understand and respond to risks to natural populations and ecosystems. For over 100 years, scientists have identified "life history invariants" (LHI) representing pairs of parameters whose ratio is theorized to be constant across species. LHI then promise to allow prediction of many parameters from field measurements of a few important traits. Using LHI in this way, however, neglects any residual patterns in parameters when making predictions. We therefore apply a multivariate model for eight variables (seven parameters and temperature) in over 32,000 fishes, and include taxonomic structure for residuals (with levels for class, order, family, genus, and species). We illustrate that this approach predicts variables probabilistically for taxa with many or few data. We then use this model to resolve three questions regarding life history parameters in fishes. Specifically we show that (1) on average there is a 1.24% decrease in the Brody growth coefficient for every 1% increase in maximum size; (2) the ratio of natural mortality rate and growth coefficient is not an LHI but instead varies systematically based on the timing of maturation, where movement along this life history axis is predictably correlated with species taxonomy; and (3) three variables must be known per species to precisely predict remaining life history variables. We distribute our predictive model as an R package, FishLife, to allow future life history predictions for fishes to be conditioned on taxonomy and life history data for fishes worldwide. This package also contains predictions (and predictive intervals) for mortality, maturity, size, and growth parameters for all described fishes.

Density-dependent changes in effective area occupied for sea-bottom-associated marine fishes.

The spatial distribution of marine fishes can change for many reasons, including density-dependent distributional shifts. Previous studies show mixed support for either the proportional-density model (PDM; no relationship between abundance and area occupied, supported by ideal-free distribution theory) or the basin model (BM; positive abundance-area relationship, supported by density-dependent habitat selection theory). The BM implies that fishes move towards preferred habitat as the population declines. We estimate the average relationship using bottom trawl data for 92 fish species from six marine regions, to determine whether the BM or PDM provides a better description for sea-bottom-associated fishes. We fit a spatio-temporal model and estimate changes in effective area occupied and abundance, and combine results to estimate the average abundance-area relationship as well as variability among taxa and regions. The average relationship is weak but significant (0.6% increase in area for a 10% increase in abundance), whereas only a small proportion of species-region combinations show a negative relationship (i.e. shrinking area when abundance increases). Approximately one-third of combinations (34.6%) are predicted to increase in area more than 1% for every 10% increase in abundance. We therefore infer that population density generally changes faster than effective area occupied during abundance changes. Gadiformes have the strongest estimated relationship (average 1.0% area increase for every 10% abundance increase) followed by Pleuronectiformes and Scorpaeniformes, and the Eastern Bering Sea shows a strong relationship between abundance and area occupied relative to other regions. We conclude that the BM explains a small but important portion of spatial dynamics for sea-bottom-associated fishes, and that many individual populations merit cautious management during population declines, because a compressed range may increase the efficiency of harvest.

Hierarchical analysis of taxonomic variation in intraspecific competition across fish species.

The nature and intensity of intraspecific competition can vary greatly among taxa, yet similarities in these interactions can lead to similar population dynamics among related organisms. Variation along the spectrum of intraspecific competition, with contest and scramble competition as endpoints, leads to vastly different responses to population density. Here we investigated the diversity of intraspecific competition among fish species, predicting that functional forms of density-dependent reproduction would be conserved in related taxa. Using a hierarchical model that links stock-recruitment parameters among populations, species, and orders, we found that the strength of overcompensation, and therefore the type of intraspecific competition, is tightly clustered within taxonomic groupings, as species within an order share similar degrees of compensation. Specifically, species within the orders Salmoniformes and Pleuronectiformes exhibited density dependence indicative of scramble competition (overcompensation) while the orders Clupeiformes, Gadiformes, Perciformes, and Scorpaeniformes exhibited dynamics consistent with contest competition (compensation). Maximum potential recruitment also varied among orders, but with less clustering across species. We also tested whether stock-recruitment parameters correlated with maximum body length among species, but found no strong relationship. Our results suggest that much of the variation in the form of density-dependent reproduction among fish species may be predicted taxonomically due to evolved life history traits and reproductive behaviors.

Space-time investigation of the effects of fishing on fish populations.

Species distribution models (SDMs) are important statistical tools for obtaining ecological insight into species-habitat relationships and providing advice for natural resource management. Many SDMs have been developed over the past decades, with a focus on space- and more recently, time-dependence. However, most of these studies have been on terrestrial species and applications to marine species have been limited. In this study, we used three large spatio-temporal data sources (habitat maps, survey-based fish density estimates, and fishery catch data) and a novel space-time model to study how the distribution of fishing may affect the seasonal dynamics of a commercially important fish species (Pacific Dover sole, Microstomus pacificus) off the west coast of the USA. Dover sole showed a large scale change in seasonal and annual distribution of biomass, and its distribution shifted from mid-depth zones to inshore or deeper waters during late summer/early fall. In many cases, the scale of fishery removal was small compared to these broader changes in biomass, suggesting that seasonal dynamics were primarily driven by movement and not by fishing. The increasing availability of appropriate data and space-time modeling software should facilitate extending this work to many other species, particularly those in marine ecosystems, and help tease apart the role of growth, natural mortality, recruitment, movement, and fishing on spatial patterns of species distribution in marine systems.

The importance of spatial models for estimating the strength of density dependence.

Identifying the existence and magnitude of density dependence is one of the oldest concerns in ecology. Ecologists have aimed to estimate density dependence in population and community data by fitting a simple autoregressive (Gompertz) model for density dependence to time series of abundance for an entire population. However, it is increasingly recognized that spatial heterogeneity in population densities has implications for population and community dynamics. We therefore adapt the Gompertz model to approximate, local densities over continuous space instead of population-wide abundance, and allow productivity to vary spatially using Gaussian random fields. We then show that the conventional (nonspatial) Gompertz model can result in biased estimates of density dependence (e.g., identifying oscillatory dynamics when not present) if densities vary spatially. By contrast, the spatial Gompertz model provides accurate and precise estimates of density dependence for a variety of simulation scenarios and data availabilities. These results are corroborated when comparing spatial and nonspatial models for data from 10 years and -100 sampling stations for three long-lived rockfishes (Sebastes spp.) off the California, USA coast. In this case, the nonspatial model estimates implausible oscillatory dynamics on an annual time scale, while the spatial model estimates strong autocorrelation and is supported by model selection tools. We conclude by discussing the importance of improved data archiving techniques, so that spatial models can be used to reexamine classic questions regarding the existence and magnitude of density. dependence in wild populations.

Using spatiotemporal species distribution models to identify temporally evolving hotspots of species co-occurrence.

Identifying spatiotemporal hotspots is important for understanding basic ecological processes, but is particularly important for species at risk. A number of terrestrial and aquatic species are indirectly affected by anthropogenic impacts, simply because they tend to be associated with species that are targeted for removals. Using newly developed statistical models that allow for the inclusion of time-varying spatial effects, we examine how the co-occurrence of a targeted and nontargeted species can be modeled as a function of environmental covariates (temperature, depth) and interannual variability. The nontarget species in our case study (eulachon) is listed under the U.S. Endangered Species Act, and is encountered by fisheries off the U.S. West Coast that target pink shrimp. Results from our spatiotemporal model indicated that eulachon bycatch risk decreases with depth and has a convex relationship with sea surface temperature. Additionally, we found that over the 2007-2012 period, there was support for an increase in eulachon density from both a fishery data set (+40%) and a fishery-independent data set (+55%). Eulachon bycatch has increased in recent years, but the agreement between these two data sets implies that increases in bycatch are not due to an increase in incidental targeting of eulachon by fishing vessels, but because of an increasing population size of eulachon. Based on our results, the application of spatiotemporal models to species that are of conservation concern appears promising in identifying the spatial distribution of environmental and anthropogenic risks to the population.

A Bayesian approach to identifying and compensating for model misspecification in population models.

State-space estimation methods are increasingly used in ecology to estimate productivity and abundance of natural populations while accounting for variability in both population dynamics and measurement processes. However, functional forms for population dynamics and density dependence often will not match the true biological process, and this may degrade the performance of state-space methods. We therefore developed a Bayesian semiparametric state-space model, which uses a Gaussian process (GP) to approximate the population growth function. This offers two benefits for population modeling. First, it allows data to update a specified "prior" on the population growth function, while reverting to this prior when data are uninformative. Second, it allows variability in population dynamics to be decomposed into random errors around the population growth function ("process error") and errors due to the mismatch between the specified prior and estimated growth function ("model error"). We used simulation modeling to illustrate the utility of GP methods in state-space population dynamics models. Results confirmed that the GP model performs similarly to a conventional state-space model when either (1) the prior matches the true process or (2) data are relatively uninformative. However, GP methods improve estimates of the population growth function when the function is misspecified. Results also demonstrated that the estimated magnitude of "model error" can be used to distinguish cases of model misspecification. We conclude with a discussion of the prospects for GP methods in other state-space models, including age and length-structured, meta-analytic, and individual-movement models.

Modeling structured population dynamics using data from unmarked individuals.

The study of population dynamics requires unbiased, precise estimates of abundance and vital rates that account for the demographic structure inherent in all wildlife and plant populations. Traditionally, these estimates have only been available through approaches that rely on intensive mark-recapture data. We extended recently developed N-mixture models to demonstrate how demographic parameters and abundance can be estimated for structured populations using only stage-structured count data. Our modeling framework can be used to make reliable inferences on abundance as well as recruitment, immigration, stage-specific survival, and detection rates during sampling. We present a range of simulations to illustrate the data requirements, including the number of years and locations necessary for accurate and precise parameter estimates. We apply our modeling framework to a population of northern dusky salamanders (Desmognathus fuscus) in the mid-Atlantic region (USA) and find that the population is unexpectedly declining. Our approach represents a valuable advance in the estimation of population dynamics using multistate data from unmarked individuals and should additionally be useful in the development of integrated models that combine data from intensive (e.g., mark-recapture) and extensive (e.g., counts) data sources.

Assessing the quality of life history information in publicly available databases.

Single-species life history parameters are central to ecological research and management, including the fields of macro-ecology, fisheries science, and ecosystem modeling. However, there has been little independent evaluation of the precision and accuracy of the life history values in global and publicly available databases. We therefore develop a novel method based on a Bayesian errors-in-variables model that compares database entries with estimates from local experts, and we illustrate this process by assessing the accuracy and precision of entries in FishBase, one of the largest and oldest life history databases. This model distinguishes biases among seven life history parameters, two types of information available in FishBase (i.e., published values and those estimated from other parameters), and two taxa (i.e., bony and cartilaginous fishes) relative to values from regional experts in the United States, while accounting for additional variance caused by sex- and region-specific life history traits. For published values in FishBase, the model identifies a small positive bias in natural mortality and negative bias in maximum age, perhaps caused by unacknowledged mortality caused by fishing. For life history values calculated by FishBase, the model identified large and inconsistent biases. The model also demonstrates greatest precision for body size parameters, decreased precision for values derived from geographically distant populations, and greatest between-sex differences in age at maturity. We recommend that our bias and precision estimates be used in future errors-in-variables models as a prior on measurement errors. This approach is broadly applicable to global databases of life history traits and, if used, will encourage further development and improvements in these databases.

Rigorous meta-analysis of life history correlations by simultaneously analyzing multiple population dynamics models.

Correlations among life history parameters have been discussed in the ecological literature for over 50 years, but are often estimated while treating model estimates of demographic rates such as natural mortality (M) or individual growth (k) as "data." This approach fails to propagate uncertainty appropriately because it ignores correlations in estimation errors between parameters within a species and differences in estimation error among species. An improved alternative is multi-species mixed-effects modeling, which we approximate using multivariate likelihood profiles in an approach that synthesizes information from several population dynamics models. Simulation modeling demonstrates that this approach has minimal bias, and that precision improves with increased number of species. As a case study, we demonstrate this approach by estimating M/k for 11 groundfish species off the U.S. West Coast using the data and functional forms on which pre-existing, peer-reviewed, population dynamics models are based. M/k is estimated to be 1.26 for Pacific rockfishes (Sebastes spp.), with a coefficient of variation of 76% for M given k. This represents the first-ever estimate of correlations among life history parameters for marine fishes using several age-structured population dynamics models, and it serves as a standard for future life history correlation studies. This approach can be modified to provide robust estimates of other life history parameters and correlations, and requires few changes to existing population dynamics models and software input files for both marine and terrestrial species. Specific results for Pacific rockfishes can be used as a Bayesian prior for estimating natural mortality in future fisheries management efforts. We therefore recommend that fish population dynamics models be compiled in a global database that can be used to simultaneously analyze observation-level data for many species in life history meta-analyses.

Spatial variation buffers temporal fluctuations in early juvenile survival for an endangered Pacific salmon.

Spatial, phenotypic and genetic diversity at relatively small scales can buffer species against large-scale processes such as climate change that tend to synchronize populations and increase temporal variability in overall abundance or production. This portfolio effect generally results in improved biological and economic outcomes for managed species. Previous evidence for the portfolio effect in salmonids has arisen from examinations of time series of adult abundance, but we lack evidence of spatial buffering of temporal variability in demographic rates such as survival of juveniles during their first year of life. We therefore use density-dependent population models with multiple random effects to represent synchronous (similar among populations) and asynchronous (different among populations) temporal variability as well as spatial variability in survival. These are fitted to 25 years of survey data for breeding adults and surviving juveniles from 15 demographically distinct populations of Chinook salmon (Oncorhynchus tshawytscha) within a single metapopulation in the Snake River in Idaho, USA. Model selection identifies the most support for the model that included both synchronous and asynchronous temporal variability, in addition to spatial variability. Asynchronous variability (log-SD = 0·55) is approximately equal in magnitude to synchronous temporal variability (log-SD = 0·67), but much lower than spatial variability (log-SD = 1·11). We also show that the pairwise correlation coefficient, a common measure of population synchrony, is approximated by the estimated ratio of shared and total variance, where both approaches yield a synchrony estimate of 0·59. We therefore find evidence for spatial buffering of temporal variability in early juvenile survival, although between-population variability that persists over time is also large. We conclude that spatial variation decreases interannual changes in overall juvenile production, which suggests that conservation and restoration of spatial diversity will improve population persistence for this metapopulation. However, the exact magnitude of spatial buffering depends upon demographic parameters such as adult survival that may vary among populations and is proposed as an area of future research using hierarchical life cycle models. We recommend that future sampling of this metapopulation employ a repeated-measure sampling design to improve estimation of early juvenile carrying capacity.