(C) PLOS One [1]. This unaltered content originally appeared in journals.plosone.org.
Licensed under Creative Commons Attribution (CC BY) license.
url:https://journals.plos.org/plosone/s/licenses-and-copyright

------------



Tracking and forecasting community responses to climate perturbations in the California Current Ecosystem

['Mary E. Hunsicker', 'Fish Ecology Division', 'Northwest Fisheries Science Center', 'National Marine Fisheries Service', 'National Oceanic', 'Atmospheric Administration', 'Newport', 'Oregon', 'United States Of America', 'Eric J. Ward']

Date: 2022-05

Ocean ecosystems are vulnerable to climate-driven perturbations, which are increasing in frequency and can have profound effects on marine social-ecological systems. Thus, there is an urgency to develop tools that can detect the response of ecosystem components to these perturbations as early as possible. We used Bayesian Dynamic Factor Analysis (DFA) to develop a community state indicator for the California Current Ecosystem (CCE) to track the system’s response to climate perturbations, and to forecast future changes in community state. Our key objectives were to (1) summarize environmental and biological variability in the southern and central regions of the CCE during a recent and unprecedented marine heatwave in the northeast Pacific Ocean (2014–2016) and compare these patterns to past variability, (2) examine whether there is evidence of a shift in the community to a new state in response to the heatwave, (3) identify relationships between community variability and climate variables; and (4) test our ability to create one-year ahead forecasts of individual species responses and the broader community response based on ocean conditions. Our analysis detected a clear community response to the marine heatwave, although it did not exceed normal variability over the past six decades (1951–2017), and we did not find evidence of a shift to a new community state. We found that nitrate flux through the base of the mixed layer exhibited the strongest relationship with species and community-level responses. Furthermore, we demonstrated skill in creating forecasts of species responses and community state based on estimates of nitrate flux. Our indicator and forecasts of community state show promise as tools for informing ecosystem-based and climate-ready fisheries management in the CCE. Our modeling framework is also widely applicable to other ecosystems where scientists and managers are faced with the challenge of managing and protecting living marine resources in a rapidly changing climate.

Data Availability: All time series collected through the CalCOFI and RREAS surveys are available through the ERDDAP data server ( https://coastwatch.pfeg.noaa.gov/erddap/index.html ). The seabird and sea lion time series or points of contact are available through the California Current Integrated Ecosystem Assessment (CCIEA) web dashboard ( https://www.integratedecosystemassessment.noaa.gov/regions/california-current ). The data is third-party data and the authors had no special access or privileges that others would not have.

This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

The goal of our study is to build on this set of novel statistical tools to develop a model of the CCE state that can both track and forecast ecosystem changes in response to climate perturbations. More specifically, we expand the Bayesian implementation of DFA to test the community response to environmental variables within the modeling framework and to develop near-term forecasts of future community states. Using climate and biological data from the central and southern regions of the CCE, our specific objectives were to: (1) summarize environmental and biological variability during 2014–2016 marine heatwave and compare these patterns to past variability; (2) examine whether there is evidence of departures from previous climate patterns and of switches to a new community state during the heatwave; (3) identify relationships, if any, between community variability and climate variables; and (4) test our ability to create one-year ahead forecasts of species responses and the community state based on environmental information. While the focus of our study is the CCE, the approach applied here is widely applicable to the myriad marine ecosystems worldwide that are vulnerable to a rapidly changing climate.

A challenge in summarizing ecosystem responses to perturbations is that time series used to characterize the ecosystem often involve tens to hundreds of variables (species or climate indices); there is often some degree of asynchrony among time series (unevenly or irregularly spaced), and further, each is corrupted by the presence of observation errors. Disentangling these sources of error and separating the signal from the noise is statistically challenging. Traditionally, tools such as Principal Components Analysis (PCA) or nonmetric multidimensional scaling have often been used for identifying leading patterns of variability in multivariate datasets (e.g., [ 36 , 37 ]); however, these approaches are ill-suited to the analysis of time series data that are autocorrelated or non-stationary [ 38 ]. An alternative approach, Dynamic Factor Analysis (DFA), is better suited for identifying shared trends that can be used as a community state indicator. DFA is specifically designed for time series ordination, and avoids many of the problems associated with other multi-variate approaches [ 39 ]. When applied to a collection of multivariate time series, inference in DFA models focuses on estimating a smaller number of temporal patterns (’trends’) that best capture the variation observed. The observed data are then treated as a mixture of these trends [ 40 ]. Ward et al. [ 40 ] recently developed a Bayesian implementation of DFA that offers added flexibility in model aspects over conventional approaches; examples include allowing for extreme “black swan” events (rare and difficult to predict events) [ 41 ] and trend processes that do not follow a random walk. Output from these Bayesian DFA models can also be used to estimate the probability of extreme events occurring or switches among contrasting system states. In the first application of this new method, Litzow et al. [ 42 ] examined shared trends of climate and biology time series in the Gulf of Alaska. Their study did not detect evidence for wholesale community reorganization during the recent northeast Pacific marine heatwave; however, their findings indicated potential for new patterns of ecosystem functioning with continued warming of ocean temperatures.

Indicators of community or ecosystem state are valuable tools for tracking climate-related changes in ecosystem functioning and evaluating those changes within the context of past climate perturbations [ 32 ]. Moreover, combining long-term monitoring surveys and data with modeling frameworks that summarize information across taxa and life stages that respond quickly to climate perturbations could provide early detection of an ecosystem shifting into a novel state. Early detection of such shifts would benefit ecosystem-based and climate-ready fisheries management strategies aimed at mitigating possible deleterious ecological and socio-economic outcomes. There is also a pressing need for forecasts of future ecosystem states to support forward-looking management of LMRs [ 33 – 35 ], including assessments of risk. As climate models and forecasts of ocean conditions continue to improve, there are burgeoning opportunities to develop and test methods that could provide near-term forecasts of community state in relation to ocean conditions.

However, between 2014 and 2016 these ecosystems experienced a marine heatwave that involved the warmest sea surface temperature (SST) and heat content anomalies that had ever been observed over large areas of the North Pacific, with SST anomalies over 6°C [ 13 , 14 ]. It was one of the most extreme heatwaves globally in its combined magnitude, spatial scale, and duration [ 1 , 2 ], and the intense, persistent warming has been attributed to a combination of natural and anthropogenic forcing [ 15 , 16 ]. Several studies have documented myriad biological responses to this event. For example, within the California Current Ecosystem (CCE), there were mass strandings of marine mammals [ 17 ], increased whale entanglements due to shifting prey sources [ 18 ], mass mortality events for marine seabirds [ 17 , 19 , 20 ], a record-breaking domoic acid outbreak [ 21 ], shifts in pelagic macronekton and micronekton communities and species richness [ 22 – 24 ], irruptions of previously rare fishes and invertebrates throughout the California Current [ 25 – 28 ], and extraordinarily high recruitment of rockfishes (genus Sebastes) [ 29 , 30 ] and northern anchovy (Engraulis mordax) [ 31 ]. Yet, to date, there have been few quantitative studies of how the marine heatwave impacted the broader CCE community at multiple trophic levels, and therefore the importance of this extreme event for community-wide patterns of variability, and the persistence of the community response, remains largely unknown.

Climate-driven shifts in community structure tend to involve rapid change across multiple populations that result in switches between contrasting community assemblages that may then persist for decades. A growing number of studies have documented community reorganizations in response to climate drivers (e.g., [ 3 – 7 ]). One of the best-known examples is the widespread northeast Pacific community reorganization that followed the 1976/1977 shift in the Pacific Decadal Oscillation from a cold to warm regime [ 8 , 9 ]. The abrupt change from a cool to warm ocean regime had dramatic implications on ecosystem functioning and living marine resources (LMRs) throughout the region [ 7 , 10 – 12 ]. Since then, northeast Pacific marine ecosystems have experienced several interannual or decadal perturbations that do not appear to have resulted in community-wide shifts of similar magnitude.

Climate perturbations can have strong impacts on ocean ecosystems that in turn affect social and economic components of human communities. These effects may be exacerbated when changes in ocean conditions are more extreme, such as during marine heatwaves (prolonged events of anomalously warm ocean waters). The increasing attention on these extreme events and their impacts (e.g., [ 1 , 2 ]) has invigorated a push for tools that can track and detect as early as possible the response of marine communities to climate-driven perturbations. Early detection, and moreover, near-term forecasts of community shifts could help scientists, managers, and stakeholders better prepare for and respond to the potential consequences of such shifts.

While a wide variety of multivariate or univariate time series methods could be applied to our observed time series to generate forecasts, our objectives were to develop simultaneous estimates of both the community state (i.e., the DFA trend value) and the raw time series (i.e., individual time series summarized by the biology DFA model). We evaluated the ability of our DFA models to generate short-term (one year lead-time) forecasts of community state by first evaluating whether the performance of the biology DFA model was improved when climate time series were included as covariates in the model. If climate time series were found to better explain the variability in the biology time series, these relationships could potentially be used to forecast community trends. For our analysis, we ran the DFA on a subset of the biology data overlapping in time with the climate dataset, i.e., 1981–2017, to make out-of-sample predictions. We used the same LFO-CV procedure described above, with the same forecast period (2009–2017) to compare the biology models with and without a single climate covariate (see S2 Table for all model formulations). In this case, the model used biological and climate data from all preceding years and climate data from the year to be forecast. The six climate covariates from the southern region and the central region of the CCE (12 total) were tested in this analysis. Once the best-supported biology-covariate model was identified, we used that model to make predictions of the community state (i.e., DFA trend value) in 2018 using climate data from that same year and the raw time series of the individual species (i.e., the biology time series summarized by the DFA model). We evaluated forecast skill based on the prediction errors of individual species time series and by comparing the forecasts for 2008–2018 to the 2008–2018 trend values estimated from the biology-covariate model that only included data prior to the forecast year.

After identifying the best-supported DFA model for the climate and biological datasets, we conducted a post-hoc examination of outlier detection and regime shifts. For outlier detection of black swan events, we implemented a method similar to that described in Anderson et al. [ 41 ] and applied it to the climate and biology time series. This approach relies on first differencing the posterior trend mean estimates of the climate and biology trends, x t − x t−1 and then applies a normal density function to identify year-over-year changes that were unlikely to have arisen from a normal distribution (given the process variance). Probabilities can then be assigned to the deviations in each year (e.g., ‘there is a 1:1000 chance of observing a deviation similar to that estimated in year t’). As described in Ward et al. [ 40 ], the presence of regimes can also be estimated by applying hidden Markov models (HMM) to the estimated state indices from a DFA. We evaluated support for regimes and alternate states by using the posterior trend estimates from each model as input. The Bayesian Leave-One-Out Cross Validation information criterion (LOO-CV) [ 61 ] was used to identify the data support for the number of regimes (n = 1–3). The model with the lowest LOO-CV value is deemed the best model.

For each model formulation, we applied the LFO-CV method by first fitting the model to all years of data prior to year T (i.e., training data, years 1, 2, …, (T-1)) and then using the fitted model to predict the trend value in year T (i.e., test data). We repeated this process for 10 years, starting with 2017 as year T and working back to 2008, and then calculated the expected log predictive density (ELPD) across those time steps. The climate and biology models with the highest ELPD were deemed the best supported models. The LFO-CV is a preferred method for evaluating future predictive performance of Bayesian models because it properly accounts for time series structure, and unlike other Bayesian cross-validation methods, does not produce overly optimistic estimates [ 60 ].

In addition, we used LFO-CV to identify the most appropriate error structure for the climate dataset—specifically whether the times series had equal (shared) or unequal (unique) observation errors. For the biology models, we assumed the observation errors were unique by dataset, and our estimates of survey variance supported this assumption.

We ran the DFA on climate datasets (1981–2017) and biological datasets (1951–2017) across the southern and central regions of the California Current combined. Running the analysis at this spatial scale allowed us to capture the broader community response to climate perturbations, compared to running models on each multivariate dataset independently (e.g., time series from a single survey). There are a number of ways to evaluate predictive accuracy of these models. The commonly used Leave-One-Out Cross-Validation (LOO-CV), for example holds each observation out in turn and predictions are made from the remaining data. As our focus was on the temporal nature of the data and forecasting component, we implemented a variant of k-fold cross validation and treated individual years as unique ‘folds’. Because our objectives involved evaluating these models for future predictions, we implemented the Leave-Future-Out Cross Validation Information Criterion (LFO-CV) [ 60 ]. We used this approach to identify data support for (1) the number of latent DFA trends (n = 1–3), (2) first-order autoregressive AR(1) coefficients on the trends ( ϕ estimated with a Normal(0,1) prior), (3) Student-t deviations (i.e., evidence of extreme events, using a prior on the MVT degrees of freedom parameter, v, of v ~ Gamma(2, 0.1)), and (4) a fixed versus estimated trend variance (using a prior on the standard deviation, σ w , of σ w ~ Normal(0, 1)).

Because we implemented the DFA model in a Bayesian setting, we were able to extend this model to include additional features. First, to include extreme events, we relaxed the assumption about process errors w t being drawn from a normal distribution and used a multivariate Student-t distribution (MVT) instead [ 55 ]. We also modified the process equation to consider an optional vector of AR(1) coefficients ϕ on the latent trends. x t = ϕx t − 1 + w t − 1 [ 40 ]. A final modification of the conventional DFA model is that for some models, process variances can be estimated rather than fixed at 1 (maximum likelihood approaches generally use this constraint for identifiability). As implemented in Stan [ 51 , 56 , 57 ], we conducted estimation with three chains, with a warm-up of 2000 samples, followed by 2000 iterations. We used the split-chain potential scale reduction factor [ 58 , 59 ] to assess convergence (Rhat < 1.05). Code to replicate these analyses is deployed as an R package on CRAN (‘bayesdfa’) [ 53 ] and our public GitHub repository, https://github.com/fate-ewi/bayesdfa .

We used a Bayesian version of Dynamic Factor Analysis (DFA) [ 39 , 40 ] using the software Stan and R [ 51 , 52 ] as implemented in the ‘bayesdfa’ package [ 53 ]. DFA is a multivariate statistical tool somewhat analogous to Principal Components Analyses, but for time-series data ( https://cran.r-project.org/web/packages/MARSS/vignettes/UserGuide.pdf ) [ 54 ]. For a collection of time series, the number of estimated ‘trends’ is specified a priori, and DFA estimates these latent trends as independent random walks. In mathematical form, this is expressed as where x t represents the value of latent (unobserved) trends at time t, and the process error deviations w t − 1 are generally assumed to be white noise having arisen from a multivariate normal distribution (with an identity covariance matrix for identifiability). The latent trends are mapped to the observed data through an estimated loadings matrix Z and residual error e t , where y t is the vector of observed states at time t, and the residual error terms e t are assumed to be drawn from a univariate or multivariate normal distribution. Though the covariance matrix of w t is generally fixed [ 39 ], the covariance matrix of e t can be structured; variances may be shared or not across time series, and off diagonal elements may be estimated. The parameter vector b represents optional estimated coefficients relating covariates d t to the observed response. In the context of our DFA modeling, we included climate variables as d t in models where the biological observations were used as the response y t .

We describe the methods in detail below, but in summary our work flow was to (1) apply Bayesian DFA to climate and biology datasets separately and use model selection tools to identify the best supported model and number of shared trends, (2) apply ‘black swan’ and regime detection methods to detect extreme events and alternating community states, respectively, (3) identify whether the CCE community state was strongly correlated with the climate time series (compare performance of the biology models with/without environmental covariates), and (4) evaluate our skill at making predictions of community state and individual species variables. These four steps map on to the four study objectives outlined in the introduction.

The biology time series included in our analysis were selected based on three criteria: first, the measured variables would be expected to show rapid (0- to 1-year lag) responses to climate variability; second, the time series could be updated with no more than one year lag for processing time to increase the speed at which biological responses to perturbation could be detected; and third, the time series were at least 15 years long. A threshold of 15 years allowed us to include long time series that spanned as many climate perturbations as possible and also have enough biological time series to develop an informative indicator of community state. In addition, 15 years is a threshold that has been previously used to define "long oceanographic time series" in the California Current [ 50 ]. The biology time series that met our selection criteria (n = 38) included ichthyoplankton, pelagic young-of-the-year (juvenile fish), squid, and krill abundance; seabird productivity; and California sea lion pup body condition metrics ( S1 Fig , S1 Table ). These 38 time series were collected from four disparate ocean surveys, and span between 22 and 68 years. Datasets collected from surveys that included spatial attributes (e.g., ichthyoplankton and pelagic juvenile fish surveys) were first standardized using Generalized Additive Models to create a univariate time series for each species. While these datasets generally include spatial random sampling, the index standardization accounts for uneven distributions of effort (in space or time). Details on the standardization of individual datasets are included in S1 Appendix . In addition, the biology data were normalized with log transformations where appropriate (all zeros were changed to NAs). For example, if the time series data were assumed to be lognormally distributed (e.g., weight/count data) or the coefficient of variation was > 1, the data were log transformed. All of the time series from an individual dataset (survey) were treated the same, i.e., logged or not. More details on the biology time series used in this study and the associated data sources and log transformations are summarized in S1 Table .

Abundance data for pelagic juvenile groundfishes and invertebrates are collected on the Rockfish Recruitment and Ecosystem Assessment Survey (RREAS). Ichthyoplankton data are collected on the California Cooperative Oceanic Fisheries Investigations (CalCOFI) survey. Seabird reproductive success and California sea lion (Zalophus californianus) pup time series are collected on Southeast Farallon Island and San Miguel Island, respectively. See S1 Table and S1 Fig for detailed information on the individual time series. The base map layer was sourced from NOAA National Geophysical Data Center (2009) ETOPO1 1 Arc-Minute Global Relief Model. NOAA National Centers for Environmental Information (accessed: 19 April 2013, [ 49 ]).

In our analysis, we used oceanographic time series from the southern (n = 6) and central (n = 6) regions of the CCE, derived from a data assimilative configuration of the Regional Ocean Modeling System (ROMS) with 0.1° (~10 km) horizontal resolution and 42 terrain-following vertical levels (oceanmodeling.ucsc.edu) [ 43 ]. From the ROMS output, we generated monthly time series covering 1980–2018 for a suite of variables including sea surface temperature (SST), sea surface height (SSH), isothermal layer depth (ILD), Brunt-Väisälä frequency (BV), a coastal upwelling transport index (CUTI), and a biologically effective upwelling transport index (BEUTI). The ILD is similar to mixed layer depth and defines the depth where temperature deviates by 0.5°C from the surface value. BV is a measure of water column stratification, averaged over the upper 200 m of the water column. CUTI and BEUTI are upwelling indices that quantify vertical transport and nitrate flux through the base of the mixed layer, respectively [ 44 ]. The data were annually averaged (July-June) from the coast to 100 km offshore, with the exception of CUTI and BEUTI, which capture coastal upwelling within 75 km of shore. In the alongshore direction, we calculated averages for two regions with a division at Point Conception, California, separating the southern portion of the CCE (31–34.5°N) from the central region (34.5–40.5°N, Fig 1 ). This is in response to the recognition of Point Conception as a major biogeographic boundary for the California Current System, with differing wind and current patterns north and south of that feature [ 45 , 46 ]. The annual averages were taken from July to June to capture the influence of the El Niño–Southern Oscillation (ENSO), which peaks in winter and is the dominant mode of interannual variability influencing the California Current [ 47 ]. We developed models using ROMS output rather than empirical measurements because they provide full spatial and temporal coverage of surface and subsurface conditions, incorporate available observations, and will enable the use of ROMS forecasts to then forecast biological changes in the CCE. This ocean model is constrained by available satellite and in situ observations to improve its fidelity to nature and has been validated against independent in situ observations [ 43 , 48 ]. Output from this model has been widely used to characterize CCE oceanography, its relation to large scale climate variability, and its influence over the marine ecosystem from phytoplankton to top predators (see Discussion ). More details on the oceanographic time series can be found in S1 Table and S1 Fig .

Overall, the model forecast skill of individual species parameters was moderate to high for half of the species included in the biology-BEUTI model ( S3 Table , S9 and S10 Figs). It is important to emphasize that the source of variability in predictions for each of the original time series is a mixture of the magnitude and uncertainty around the trends and loadings ( x t , Z ), and the magnitude and uncertainty in the estimated covariates ( b ). On the one hand, the time series associated with the highest predictive skill (i.e., lowest prediction errors) included seabird reproductive success (common murre, Cassin’s auklet) and the abundance of juvenile Pacific sanddab Citharichthys sordidus, juvenile halfbanded rockfish Sebastes semicinctus, market squid, and several ichthyoplankton species ( S3 Table , S9 and S10 Figs). On the other hand, forecast skill was lowest (i.e., highest prediction errors) for the abundance of some juvenile rockfishes (chilipepper Sebastes goodei and widow rockfish Sebastes entomelas) and larval fishes (northern anchovy, mesopelagics), which is likely attributed to a lag or mismatch in the timing of the climate-biology relationships. These patterns in forecast skill are similar to those based on the biology-CUTI, -SST, and -ILD models ( S3 Table ). Lastly, the uncertainty around model predictions of species parameters appears to be driven more by the precision of the model coefficients than by the loadings on the community trend (e.g., S11 Fig ).

Given that the biology-BEUTI model was the best supported model over the null model (a model without covariates), we were interested in evaluating the ability of this covariate model to forecast the community state. Comparisons between the community state and the community state forecasts (out-of-sample estimates) indicate that we had skill in forecasting community state one year in advance ( Fig 9 ). Forecasts of the community trend values for ten additional years (2008–2017, Fig 10 ) also indicate that we had some skill for many of the years tested. There are wide confidence intervals around the forecasts; however, given our methodology we can expect that the credible intervals around the trend forecast will be larger than the historical credible intervals ( Fig 9 ). Forecasts have more uncertainty than historical values because the variance of a random walk increases linearly with time [ 66 , 67 ]. Furthermore, our credible intervals are increased because we are additionally (1) propagating full parameter uncertainty across the MCMC draws projecting it, and (2) using a Student-t distribution, which has heavy tails and therefore makes the uncertainty intervals wider than if we used normal distribution.

In comparing models of the biological response with and without climate covariates, we found that several biology models with climate predictors outperformed the biology models that did not include covariates ( S2 Table ). The climate covariate resulting in the best future predictions of community state was BEUTI (central region), followed by CUTI (central region) ( Table 2 , see S2 Table for all models). The coefficients linking BEUTI to observed time series indicate strong, positive relationships between nitrate flux and the reproductive success of seabirds and the abundance of krill in the central California Current ( Fig 8 ). They also indicate strong, negative relationships between nitrate flux and the abundance of juvenile/adult Pacific sardine and larval northern anchovy ( Fig 8 ). The remaining biology-BEUTI relationships were moderate (e.g., ichthyoplankton, market squid Doryteuthis opalescens) to weak (e.g., rockfish spp., Fig 8 ). The biology-CUTI model was similar to the biology-BEUTI model with respect to model structure and estimated species loadings. The estimated coefficients in the CUTI model ( S5 Fig ) also show a similar pattern to those in the biology-BEUTI model. The remaining covariate models only showed weak climate-biology relationships (e.g., S6 and S7 Figs).

While this model provided slight support for heavy-tailed Student-t deviations in the latent trend ( S4 Fig ), we did not detect any black swan events in the community state. We note that the community response to two strong El Niño events (1982–1983 and 1997–1998) and to unusually low productivity conditions (2005) in the central CCE appear similar in magnitude and duration to the response to the 2014–2016 marine heatwave, although the directions of the responses were opposite ( Fig 5 ). Our regime detection analysis also captured the change in the central CCE community in the mid to late 2000s ( Fig 7 ), which may be associated with the large changes in the reproductive success of multiple seabirds (e.g., Cassin’s auklet Ptychoramphus aleuticus, common murre Uria aalge, Brandt’s cormorant Urile penicillatus) and in sea lion pup births around that time ( S1 Fig ). These taxa may have been impacted by changes in the abundance or availability of important prey items resulting from unproductive ocean conditions in the central CCE in 2005 and the below normal SSTs associated with the 2007–2008 La Niña Event [ 63 – 65 ].

The best model invoked two states, and the median probability (and 95% credible intervals) of being in one state versus the other is shown. The figure indicates that ecosystem did not shift into a new state following the marine heatwave.

The estimated trend from this biology DFA model demonstrates a potential shift in community state in the mid-1960s, although there is considerable uncertainty around the trend during this early portion of the time series, likely due to the limited number of observations (ichthyoplankton only) pre-dating the 1970s ( Fig 5 , S1 Fig ). The community state appears to be relatively stable from the late 1970s through the early 2000s, and the trend reached a peak around 2013–2015. Evidence of a community shift early in the time series is supported by our regime detection analysis, which demonstrated that a two-state model best described the latent trend (LOO-CV: one-state = 216.4, two-state = 11.8, three-state = 41.8, Fig 7 ). This shift coincides with a strong increase in the abundance of a few species during that period, including eared blacksmelt (Lipolagus ochotensis), slender blacksmelt (Bathylagus pacificus), northern lampfish (Stenobrachius leucopsarus), which are cool water associated mesopelagic species, as well as a rise in northern anchovy (Engraulis mordax) abundance prior to the shift ( S1 Fig ). Our analysis does not document a shift to a novel community state in response to the recent marine heatwave.

Only time series with ≥ 90% of the loading distributions above or below zero are shown). Loadings with darker shading indicate time series loading most strongly on the biology trend. Cal. = California, Juv. = juvenile fish stage, Juv./adult = juvenile and adult fish stages combined, all other fish are larval fish. See S1 Table and S1 Fig for times series details.

The best model for community variability among our biological time series was also a one-trend model (Model 13 in Table 1 , Fig 5 ). The model formulation was similar to the best climate model, except the observation variances were unique by dataset (survey) and not individual time series. We note that the top two models (Model 13 and 14) showed similar predictive accuracy (Δ ELPD < 1) and only differed with respect to whether the process variance was fixed at 1 or estimated. Here we only show results for the model with a fixed process variance. The biology showed strong coherence in community signal; a majority of the time series (31 of 38) loaded strongly (probability > 0.9) on the single trend and most of them demonstrated loadings in the same direction ( Fig 6 ). The magnitude and direction of the estimated loadings were consistent with the observed high relative abundance of most juvenile groundfishes (rockfish, flatfish), squid, krill, and some ichthyoplankton species during the marine heatwave, and suggest that the reproductive success of some seabird species was higher around the time of the heatwave as well. The few time series loading in the other direction on the trend indicated a reduction in sea lion pup growth rate and lower abundances of juvenile/adult Pacific sardine Sardinops sagax and some ichthyoplankton species (e.g., larval northern anchovy and Pacific hake Merluccius productus) associated with the heatwave.

The best model invoked two states, and the median probability (and 95% credible intervals) of being in one state versus the other is shown. The figure reflects the probability of being in the state associated with warmer conditions versus one with cold conditions.

The climate state during the marine heatwave, as indicated by the DFA trend, was within the bounds of previous observations. While there was support in the best model for heavy-tailed deviations in the climate trend (i.e., Student-t deviations S3 Fig ), our post-hoc examination of outliers detected a single extreme event in the climate state in mid-1998 to mid-1999 (threshold = 0.001), when there was a shift from strong El Niño (1997–98) to strong La Niña (1998–1999) conditions, and not around the time of the heatwave. Application of the Bayesian HMM to the climate trend most supported the presence of two hidden states, reflecting the probability of being in a state associated with warmer conditions versus one with cooler conditions (LOO-CV: one-state = 129.1, two-state = 9.4, three-state = 27.2, Fig 4 ). The LOO-CV did not provide support for a shift to a third novel climate state in the southern and central regions of the CCE during the marine heatwave, however there is a shift back to the previously observed warm state during the marine heatwave.

The model with the highest predictive accuracy (ELPD) of the climate state in the southern and central regions of the CCE was a one-trend DFA model (Model 1 in Table 1 , Fig 2 ). This model included unique observation variances across the six time series, support for heavy-tailed deviations of the latent trend, an AR(1) coefficient on the trend ( S2 Fig ), and an estimated trend variance. Overall, the trend captured a well-documented cooling period in the CCE between 1980 and 2010 (e.g., [ 62 ]) as well as strong El Niño events (e.g., 1982–1983, 1997–1998, 2015–2016) and the 2014–2016 marine heatwave. The trends and loadings indicate that these events were generally associated with weaker upwelling, reduced mixed layer depth, low nutrient flux, and warm, stratified waters (Figs 2 and 3 ). All but one of the climate time series (central ILD) were strongly associated with the single trend, i.e., at least 90% of the loading posterior distributions associated with each time series were above or below zero ( Fig 3 ). The SST, SSH, and BV frequency (water column stratification) time series from the southern and central regions of the CCE loaded positively on this trend ( Fig 3 ). The BEUTI and CUTI time series from both regions of the CCE and the ILD time series from the central region loaded negatively on the trend ( Fig 3 ).

Discussion

We applied a novel set of statistical tools to data from the southern and central regions of the CCE to document the community response to climate perturbations over the past six decades and to create near-term forecasts of community state. Our analysis detected a community response to the 2014–2016 northeast Pacific marine heatwave; however, it did not exceed normal variability within the study timeframe or result in a shift to a novel community state, based on the biological time series investigated here. We identified relationships between community state and multiple climate variables, with nitrate flux through the base of the mixed layer having the strongest correspondence with individual species time series and the shared trend in community variability. Moreover, we demonstrated skill in creating simultaneous one-year lead time forecasts of species responses and community state.

Long-term changes in community state Many studies and anecdotal accounts have documented unexpected biological responses to the 2014–2016 northeast Pacific marine heatwave. Based on the biological time series included in our analysis, the broader CCE community demonstrated a clear response to the marine heatwave (Fig 5). However, our results do not demonstrate a widespread community reorganization beyond the archetypal community structure of this dynamic ecosystem in response to this event. Instead, the mean values for the shared trend in the biology time series, as well as for the shared climate trend, were within the range of previous observations. Many species were present during the marine heatwave that are not typically observed in the CCE. While our analysis could not include these sporadically occurring taxa, due to the large number of zero observations in the historic survey data, the exceptional presences and high abundances of those warm species did not result in a persistent signal among the species included in the DFAs. As additional years of data become available, the DFA models could reveal different outcomes from 2014–2016. However, this is unlikely given that the taxa and life stages used in both studies are known to respond quickly to changes in ocean conditions and given our assumption that the surveys are consistently sampling at the right time and location to fully characterize the short-term response. While our study did not detect a shift in community state in the southern and central CCE during the 2014–2016 heatwave, we did detect a shift in the 1960s. The 1960s shift was likely due to a regime shift previously detected in the southern California ichthyoplankton community [7]. The Peabody et al. [7] study included a much broader suite of ichthyoplankton species than our study which limits our ability to evaluate whether the species driving the shifts are consistent among studies. Previous studies have also documented a shift in response to the 1976/1977 PDO shift (e.g., [7, 68]), but our analyses did not. Our estimated biology trend is consistent with the evidence of this regime shift, however, only ichthyoplankton time series are available prior to the 1970s and there are gaps in the ichthyoplankton data from the late 1960s through the 1970s. The trend estimate therefore has higher uncertainty during this period than elsewhere in the time series. The CCE biology time series included in this study showed strong coherence in community signal in response to regional climate perturbations. Although they span across multiple trophic levels, life-history strategies, and datasets, most of the biological time series loaded in the same direction on the shared trend (Fig 6). In addition, our the CCE shared biology trend and loadings captured an unusual aspect of the 2014–2016 warming events: the abundance of several taxa, including young-of-year rockfish and anchovy, was high during the marine heatwave [22, 29, 31]. By contrast, their abundance was greatly reduced in most previous warm events, including two of the strongest El Niño events on record (1982–1983, 1997–1998) and unusually low productivity conditions (2005–2006) [69]. High abundances of young-of-year rockfish and groundfish, squid, and krill in the CCE are generally associated with more southward transport and subarctic source waters, while abundances are typically far lower in years with more subtropical waters, which are often associated with El Niño and anomalous warm events [29, 70]. The unexpectedly high abundance of these taxa in 2014–2016, despite surface-oriented marine heatwave, may be related to the prevalence of subsurface waters that were more subarctic than subtropical in origin [29] and to some concurrent strong upwelling events, particularly in spring 2015 [71, 72]. Our results were consistent with recent studies of several top predators in the CCE. The DFA trends and loadings indicate a negative response of sea lion pup growth and weight to the 2014–2016 marine heatwave, which also aligns with past work showing that sea lion pup condition covaries with abundance of forage such as larval anchovy and sardine, which provide quality prey to sustain lactation in nursing mothers [73]. Pup condition also improved at the tail end of the marine heatwave when, despite the warm water, anchovy abundance increased dramatically [31]. The trends and loadings suggest that the reproductive success of some seabirds in the central CCE was not diminished by the heatwave, although the heatwave had severe impacts seabird productivity in regions to the north [20].

Regional comparison of the marine heatwave’s effect on community state A compelling outcome of our analysis and a similar analysis applied to Alaskan species by Litzow et al. [42] is that neither detected state changes in North Pacific communities following the massive 2014–2016 marine heatwave, despite the extremely anomalous physical conditions throughout most of the basin and a litany of concurrent biological, ecological, social and economic effects (see Introduction). An important characteristic of both studies is the temporal scale of community analysis (1972–2017 for the Gulf of Alaska (GOA) and 1951–2017 for the CCE). This long temporal scale provides an important context for comparing contemporary change with the magnitude of historical community shifts. In addition, the Bayesian DFA accounts for uncertainty in the shared trends in a way that prevents premature detection of wholesale ecosystem shifts. We note that Suryan et al. [74], fitted a single-trend, non-Bayesian DFA model to a larger set of GOA biological time series (n = 187) over a shorter time span (2010–2018) and found evidence of a well-resolved shift that implied different community states during 2010–2014 and 2015–2018. The different conclusions of Suryan et al. [74] and Litzow et al. [42] studies speak to an inherent tension in retrospective analyses of community change. Limited time series availability means that analyses can be taxonomically and functionally broad (e.g., [74]), or temporally extensive (e.g., [42]), but not both. Each approach has advantages, but direct comparison between the two is difficult. Given the impacts of the 2014–2016 event, and that long-term warming combined with marine heatwaves will push the CCE into novel climate states, we must consider ecological mechanisms that might explain why these communities were apparently resilient to the marine heatwave, along with revisiting methodological details that could further clarify our results.

Environmental covariates DFA models of CCE biology that included a climate covariate performed better than models without one. Nitrate flux into the surface mixed layer (BEUTI) was the best-performing covariate for individual species in addition to the shared trend in the southern and central CCE over the past three decades. Nitrate flux had a strong positive effect on the abundance of krill and some larval fishes and on the reproductive success of seabirds, and a moderate positive effect on some pelagic juvenile fishes, squid, and sea lion pup births. Stronger upwelling magnitude (CUTI), which is correlated with nitrate flux, was the second-best predictor of community variability and had a positive albeit weaker effect on the same suite of species (S5 Fig). These findings are consistent with mechanistic understanding as upwelling increases the supply of nutrients to shallow waters and enhances the productivity at the lower trophic levels, including juvenile rockfishes [75], which affects foraging conditions for higher trophic level species, such as seabirds (e.g., [76]). BEUTI and CUTI had a strong, negative correlations with juvenile/adult Pacific sardine and larval northern anchovy. The relative abundance of Pacific sardine in coastal waters off of Central California has been shown to be lower during periods of strong upwelling [70, 77]. This trend may reflect a change in the production of Pacific sardine or a shift in their relative distribution. In addition, a negative relationship between upwelling and sardine recruitment can generally be explained by the transfer of fish larvae to offshore areas where they have low chance of survival during periods of strong equatorward flow and upwelling [78, 79]. Our understanding of the mechanisms driving anchovy population dynamics is limited [80]. Climate drivers often act in concert to influence community variability, and here we are evaluating the effects of climate variables one at a time. An important next step of this work will be to examine whether including multiple climate covariates further improves the forecast skill of the CCE biology and our community state indicator. However, the individual climate variables are collinear and share information, which affects our ability to makes inference on the covariates. Furthermore, our study is broad synthesis of community indicators and their response to climate perturbations, and should not be interpreted as replacing more detailed investigations into the drivers and mechanistic understanding of the indicators included here.

Community state forecasting skill Our approach for creating simultaneous predictions of species responses and shared ecosystem variability to ocean conditions shows promise for developing near-term forecasts of community state. Our forecasts are based on outputs from the CCE ROMS, which have been used to examine how oceanographic processes affect fish recruitment variability [81, 82] and productivity [83], species habitat suitability [84, 85], and species spatial distributions [86, 87]. The CCE ROMS also supports nowcasts of species distributions based on observed ocean conditions [88–90]. Moreover, multiple efforts are underway in the CCE and other coastal systems to use ROMS outputs to develop short-term forecasts of ocean conditions for uptake by scientists, managers, and other end-users [35, 91–93]. Here, we were able to create forecasts of community state and several individual species parameters one year in advance based on observations of a single climate variable (nitrate flux). Forecast lead times could be extended further by using forecasts of ocean conditions rather than observed conditions, and ocean temperatures in the CCE can be skillfully forecast months to a year in advance, with particularly high skill in the late winter and spring [94]. Future extensions of our work will evaluate whether different combinations of climate variables and time lags might improve our forecasting skill. Using DFA to forecast attributes of community structure in the CCE allows us to create simultaneous forecasts of trends, or ‘ecosystem state’, and raw time series. Our approach could also be applied individually to each dataset in our analysis to generate taxa-specific indicators (e.g., seabird productivity, juvenile fish abundance), though these forecasts would be expected to differ from those with the entire CCE dataset. Similarly, if ecosystem states were not a focus of inference, alternative forecast models could be applied (e.g., ARIMA or non-parametric models) [67]. Forecasts for individual time series from the DFA models used here can be seen as a mixture of the AR forecast on the estimated trends (Fig 9), and linear effects of forecasted climate variables on each time series (Fig 8). Species that have strong associations or loadings on the trend and estimated climate effects that are large in magnitude (e.g., market squid, Pacific sanddabs, shortbelly rockfish Sebastes jordani) are expected to have the most accurate predictions, while those species with weak loadings and weaker effects of climate variables (e.g., California smoothtongue (Leuroglossus stilbius) are expected to have poorer forecast performance. Nonstationary relationships are an important consideration for producing reliable ecological forecasts. While the year-to-year variability in the estimated trend did appear to be stationary in our community models (Figs 5 and 9), the autocorrelation appeared to be nonstationary with the lag-1 autocorrelation between 2000–present being significantly higher (0.82) than over the years 1981–2000 (0.23). In addition to nonstationary variance parameters, future analyses may also consider nonstationary relationships in the covariate relationships, or potential interactions between covariates. A growing number of retrospective analyses have revealed nonstationary relationships among climate and individual species or community-level variables [42, 95–99]. In the northeast Pacific Ocean, these studies have been mostly focused on Alaskan ecosystems with long time series describing climate and biological processes. The best-documented instance of nonstationary relationships among climate and biology time series in the North Pacific centers on a climate shift in the late 1980s [98]. Decades of observational data on either side of that event allow for statistically robust tests for nonstationarity that are not yet available for post-2014–2016 conditions. However, early indications from Alaska suggest the possibility that long-standing relationships between leading climate modes and individual climate and biology time series may have changed following 2014 [99].

[END]

[1] Url: https://journals.plos.org/climate/article?id=10.1371/journal.pclm.0000014

(C) Plos One. "Accelerating the publication of peer-reviewed science."
Licensed under Creative Commons Attribution (CC BY 4.0)
URL: https://creativecommons.org/licenses/by/4.0/


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

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.