Investigating the effects of climatic variability on biological diversity, productivity, and stability is key to understanding possible futures for ecosystems under accelerating climate change. A critical question is, how does climatic variability influence juvenile recruitment of different fish species and life histories that use estuaries as nursery habitats? Here we examined spatiotemporal abundance trends and environmental responses of 18 fish species that frequently spend the juvenile stage rearing in the San Francisco Estuary, CA, USA. First, we constructed multivariate autoregressive state-space models using age-0 fish abundance, freshwater flow (Flow), and sea surface temperature (SST) data collected over four decades. Next, we calculated coefficients of variation (CV) to assess portfolio effects (1) within and among species, life histories (anadromous, marine opportunist, estuarine dependent species), and the whole community, and (2) within and among regions of the estuary. We found that species abundances varied over space, across time (increasing, decreasing, dynamically stable), and in 83% of cases, in response to environmental conditions (wet/dry, cool/warm). Anadromous species responded strongly to flow in the upper estuary; marine opportunist species responded to flow and/or SST in the lower estuary; and estuarine dependent species had diverse responses across the estuary. Overall, the whole community when considered across the entire estuary had the lowest CV, and aggregate life histories and species provided strong biological insurance to the portfolio (2.4 to 3.5-fold increases in stability, respectively). Spatial insurance also increased stability, although to a lesser extent (up to 1.6-fold increases). Our study advances the notion that fish recruitment stability in estuaries is controlled by biocomplexity – life history diversity and spatiotemporal variation in the environment. However, under climate change, intensified drought and marine heatwaves may increase the risk of multiple consecutive recruitment failures by synchronizing species dynamics and trajectories via Moran effects and ultimately erode estuarine nursery function.

These data were derived from the following resources available in the public domain: (1) Department of California Fish and Game San Francisco Bay Study (URL: https://filelib.wildlife.ca.gov/Public/BayStudy/); (2) NOAA/UCSD Shore Stations Farallon Islands station (URL: https://shorestations.ucsd.edu/data-farallon/); and (3) California Department of Water Resources dayflow (URL: https://data.cnra.ca.gov/dataset/dayflow). All datasets were filtered to the same time periods (April to October 1980 to 2018) and summarized as mean annual values: (1) San Francisco Bay Study: log-transformed mean annual catch per unit effort of 18 age-0 fishes at 35 core stations ('CPUE'); (2) zscore-transformed Farallon Islands: mean annual sea surface temperature ('SST') that incorporates interpolated missing values from seasonal autoregressive integrated moving average (ARIMA) models-- see associated script; (3) Dayflow: zscore-transformed mean annual Net Delta Outflow ('Flow').

Metadata is available in the public domain: (1) Department of California Fish and Game San Francisco Bay Study (URL: https://filelib.wildlife.ca.gov/Public/BayStudy/); (2) NOAA Shore Stations Farallon Islands station (URL: https://shorestations.ucsd.edu/data-farallon/); and (3) California Department of Water Resources dayflow (URL: https://data.cnra.ca.gov/dataset/dayflow). README files are available for each data type.

R Markdown scripts are provided to reproduce model results based on the derived datasets: (1) Multivariate autoregressive state space (MARSS) model script specifies species-level model structures, fits, diagnostics, AICc comparisons, and coefficients; (2) Autoregressive integrated moving average (ARIMA) script specifies seasonal model structure and results, to impute missing daily SST values before summarizing into annual means from April-October 1980-2018.