Resilience and stability of a pelagic marine ecosystem

(a) Data collection

In order to characterize the ocean–atmospheric, hydrographic, physical and biotic conditions in the SCCS, we performed a data inventory of representative variables (electronic supplementary material, table S1). These variables were selected based on their ecological importance, the length of the time series and the completeness of the dataset. To reflect the ocean–atmospheric conditions affecting the regional climate in the area, winter averages of the Pacific decadal oscillation (PDO) [18], the tropical multivariate El Niño southern oscillation (ENSO) index (MEI), as well as the North Pacific gyre oscillation index, i.e. reflecting basin-scale variations in wind-driven upwelling and horizontal advection [19], were used. To account for hydrodynamic conditions, coastal upwelling, open-ocean (wind stress curl-driven) upwelling [20], sea-level height (detrended), alongshore transport, as well as dynamic height were used. Physical conditions were represented by spring averages of 0–100 m temperature, salinity and oxygen concentrations, as well as water column density through the average sigma theta across the regular area monitored by the California cooperative oceanic fisheries investigations (CalCOFI) programme (electronic supplementary material, figure S1). Nutrient conditions were represented by average concentrations of nitrate in the mixed layer. Since nitrate was consistently sampled only from 1984 onwards, we extended the time series backwards until 1951 based on modelled estimates. These were derived from a generalized additive model (GAM) with upwelling, temperature and sea level as predictors, demonstrating highly significant relationships overall (electronic supplementary material, table S2) and a high degree of explained deviance (79.4%).

Finally, a set of biotic variables representing the dominant species or functional groups across multiple trophic levels were compiled. In order to ensure spatio-temporal overlap, the variables were primarily based on CalCOFI data obtained in the same season (spring) and area (electronic supplementary material, figure S1). Although some species perform pronounced feeding migrations, notably albacore (Thunnus alalunga), Pacific sardine (Sardinops sagax) and Pacific hake (Merluccius productus), and prefer slightly different spawning habitats, the SCCS remains a core distribution area in which zooplankton and fish species of both warm-water (southern) and cold-water (northern) affinities overlap and interact [21–23]. The mean spring chlorophyll a in the mixed layer, spring biomass of copepods, divided into small (less than 0.5 mm) and large (greater than 0.5 mm) taxa, salps, chaetognaths and euphausiids, the latter divided into warm-water (southern) and cold-water (northern) species [21], were chosen to characterize primary producers and zooplankton, respectively. Note that owing to lack of chlorophyll a prior to 1984 we extended the time series backwards until 1951 with modelled estimates (electronic supplementary material, table S2), i.e. based on a GAM explaining a high degree of deviance from 1984 to 2010 (87.8%). The prey fish community was represented by ichthyoplankton samples of mesopelagic fish, divided into a warm-water (southern) and cold-water (northern) assemblage [22], as well as pelagic fish, the latter by jack mackerel (Trachurus symmetricus), Pacific sardine, northern anchovy (Engraulis mordax) and Pacific mackerel (Scomber japonicus) derived from analysis of ichthyoplankton samples and available stock assessments (electronic supplementary material, table S1). Note that in the absence of sardine stock assessments during the period of low stock size from 1963 to 1980, hindcast model estimates were used [24].

The dominant fish predators in the area were represented by ichthyoplankton samples of Pacific hake and spawning stock biomass (SSB) estimates of albacore. The albacore SSB estimates were extended from 1951 to 1965 based on their significant linear relationship (p < 0.001; n = 42; R² = 0.70) with the 3 year running mean of total landings (electronic supplementary material, figure S2). Furthermore, SSB estimates of Pacific hake were used to represent an aggregated predatory fish biomass. Finally, marine birds and mammal predators were represented by total seabird densities across the CalCOFI area during spring [25] and annual counts of sea lion pups in California [26]. Owing to the limited length of these time series, the variables were not included in the following principal component analysis (PCA). Instead, correlations of these variables with the main principal components (PCs) were computed and compared with the analogous loadings of all other biotic variables included in the PCA (electronic supplementary material, tables S3–S4).

(b) Long-term ecosystem dynamics

PCA was used to extract the dominant modes of variability in the time series [15]. Beforehand, missing values were replaced by 2-year averages (or 4-year averages in the case of 2 consecutive missing years). In order to assess the sensitivity of the PCA to this replacement, we performed a bootstrap routine where missing values were replaced by a random number with mean and variance corresponding to a 5-year sliding window. The mean and 95% confidence intervals of the resulting PCs after 1000 iterations were then compared with the actual PCs (electronic supplementary material, figure S3). To improve linearity all biotic variables, except stock assessment estimates, were log₁₀(x + 1) transformed. In order to study the temporal response of the biota relative to environmental drivers, we performed a PCA on biotic variables only (electronic supplementary material, table S3). As a complement, a PCA using all the time series was performed (figure 2). Furthermore, we identified potential abrupt transitions in the dominant modes of variability, characterized by the first three PCs (electronic supplementary material, figure S4), using a sequential algorithm designed to detect statistically significant shifts in the mean level and the magnitude of fluctuations in time series by using modified two-sided Student's t-tests (STARS) [27]. We used a significance level of p = 0.05 and a cut-off length of 15 years. More information on the STARS method and Excel add-in software is available online at www.beringclimate.noaa.gov.

• Download figure
• Open in new tab
• Download PowerPoint

(c) Drivers of ecosystem dynamics

We used GAMs [28] to examine the relationship between the first three PCs of the biotic PCA and a set of abiotic predictors representing potential drivers of ecosystem dynamics in the SCCS (electronic supplementary material, table S5). The following formulation with the PCs as response variables was used

where a is the intercept, s is the spline smoother function, V_i (i = 1 to n) is a number of abiotic predictors and ε is the error term. We applied a step-wise fitting routine by first fitting the models to each predictor, then sequentially refitting the models after adding the most significant predictor in each run (electronic supplementary material, table S5). Second, we minimized the generalized cross-validation criterion and used partial F-tests to find the best possible set of predictors for each model (electronic supplementary material, table S6). The spline smoother function (s) was constrained to three degrees of freedom (k = 3), to allow for nonlinearities but to restrict flexibility during model fitting. To account for potential threshold-dependent relationships, involving hysteresis [12], we compared the final GAMs with threshold GAMs (TGAMs), allowing for non-additive effects of the explanatory variables below and above a certain threshold value (Φ) estimated from the data [29]. To compare the performance of TGAMs relative to GAMs, it is necessary to account for the additional parameter used for the threshold search [29]. Thus, we used the genuine CV to compare models (electronic supplementary material, table S7), which equals the average squared leave-one-out prediction errors; the leave-one-out prediction is obtained by removing 1 year at a time from the model fitting and predicting its value from the resulting model. After fitting and selection, the final model residuals were checked for normality and serial dependency (electronic supplementary material, figure S5).

(d) Trophic cascades and modes of trophic control

Trophic cascades are generally characterized by negative relationships between adjacent trophic levels [30–32]. To examine the occurrence of trophic cascades and the degree of bottom-up versus top-down forcing, we performed linear regression and correlation analysis between adjacent trophic levels, consisting of chlorophyll a and aggregated biomass of zooplankton, small pelagic fish and predatory fish (table 1). In addition, we performed an analysis between observed nitrate availability and chlorophyll a (1984–2010), as well as a separate analysis between chlorophyll a and the sum of the two dominant zooplankton species, Calanus pacificus and Metridia pacifica [33]. Note that we excluded years with missing values in order to avoid potential bias owing to averaging.

(e) Functional complementarity and the degree of compensation

In order to investigate the degree of compensation between functional groups in each trophic level, we used correlation and the variance ratio (VR) test [34]. A VR less than 1 indicates that the sum of covariances among species is negative, hence indicating asynchrony, whereas a VR greater than 1 occurs when the sum of covariances among species is positive, i.e. indicating synchronous (correlated) dynamics among species or functional groups within each trophic level. We assessed the statistical significance of the VR through bootstrapping [35], by randomly shuffling the times-series values for each species or functional group independently in order to break any dependency between species fluctuations, as well as serial dependency (autocorrelation) in the time series. The VR statistics were compared with the resulting 95% confidence intervals (table 2). All statistical analyses were conducted using the R software version 2.15.1 (www.r-project.org).

(a) Long-term ecosystem dynamics

The temporal dynamics of the SCCS show pronounced multi-decadal variability (figure 2), as well as higher frequency fluctuations at a range of timescales (electronic supplementary material, figure S6). Based on a PCA of only biotic variables (electronic supplementary material, table S3), the dominant mode of variability (PC1; explaining 19.2% of the total variance) is characterized by moderate fluctuations until the mid-1970s, followed by abrupt decline and sudden return to positive values in the late 1990s (figure 3a). Likewise, PC2 (r² = 0.179) shows minor fluctuation during the first two decades, but more pronounced decadal variability throughout the remainder of the period (figure 3c). Finally, PC3 (r² = 0.138) illustrates higher interannual variability and peak values during the mid-1990s (figure 3e). The long-term dynamics of the PCs were robust to the replacement of missing values as the PCs were well contained within the narrow 95% confidence intervals obtained through random bootstrapping (figure 3 and electronic supplementary material, figure S3). The significant break points in PCs (e.g. 1965, 1976, 1987 and 1998) coincide with abrupt shifts documented in the North Pacific and correspond to large-scale changes in ocean–atmospheric forcing, e.g. manifest as the PDO and ENSO [36–38]. The occurrence of breakpoints indicates abrupt transitions in the SCCS but provides no guidance regarding the underlying drivers or the type of ecosystem response [12,13,16]. Nonlinear threshold models show that nitrate concentrations, wind stress curl-driven upwelling and temperature are the best single predictors of ecosystem state, represented by PC1, PC2 and PC3, respectively (electronic supplementary material, tables S5–S6). Interestingly, the response of the ecosystem to these key drivers is optimally modelled by linear relationships (figure 3b,d,f and electronic supplementary material, table S7).

• Download figure
• Open in new tab
• Download PowerPoint

The positive loading of chlorophyll a on PC1 (electronic supplementary material, table S3) is directly explained by nitrate availability, whereas a suite of interrelated hydrographic processes, e.g. oxygen, temperature and sea level (electronic supplementary material, table S8), likely underlie the negative loading of mesopelagic fish and alternating dominance of warm and cold-water euphausiids, pelagic fish and predator species with opposite loadings on PC1 [20–22,24,39,40]. Likewise, the significant positive and negative correlations of sea lions (r = 0.821, p < 0.001) and seabirds (r = −0.423, p = 0.04), respectively, with PC1 suggests opposite responses to changing hydrographic conditions (e.g. upwelling) and/or to indirect changes in the prey community [25]. The strength of curl-driven upwelling underlies the shift and opposite loadings of small and large copepods on PC2, as well as the positive effect on total biomass of zooplankton and pelagic fish [20]. The positive loadings of sardine and mesopelagic fish and the negative loadings of zooplankton on PC3 may be explained by a combination of temperature, oxygen conditions and alongshore transport (electronic supplementary material, table S6) [24,39]. Sea lions and seabirds were not correlated with PC2 and PC3 (electronic supplementary material, table S4). Furthermore, we find no evidence of trophic cascades, as adjacent trophic levels are not negatively correlated (table 1). Instead, we find positive correlations across all trophic levels, indicating bottom-up control. Trophic cascades have been shown to be associated with abrupt transitions and hysteresis in pelagic ecosystems [41,42] and are considered hallmarks of top-down control and ecological instability [43,44]. Hence, the pronounced multi-decadal variability and occurrence of significant breakpoints in the SCCS does not indicate functional changes, involving trophic cascades and hysteresis [17], but rather a linear response to positively autocorrelated climate variables [45,46]. These drivers regulate the ecosystem dynamics via bottom-up and physical control [20,47] without profound changes in functioning.

(b) Functional complementarity and the underlying mechanisms of stability

To provide a mechanistic understanding of the processes promoting resilience and stability in the SCCS, we investigated the importance of statistical averaging versus functional complementarity across multiple trophic levels. In the SCCS, and other upwelling areas worldwide, small pelagic fish, e.g. Pacific sardine and northern anchovy, typically show pronounced multi-decadal fluctuations [20,24,36,48], a feature shared by predators, e.g. albacore and Pacific hake, as well as mesopelagic fish and euphausiids [21,39,40]. The high degree of population variability is reflected by elevated mean coefficients of variation (CV), compared with lower CVs of aggregated community biomass or abundance (table 2). These findings are consistent with experimental studies on grasslands demonstrating high population, but low community variability at high diversity [49]. Furthermore, a variance-ratio test indicates significant compensatory dynamics, except for mesopelagic fish where values are just within the confidence interval (table 2).

Our results show strong temporal asynchrony between key species or functional groups within each trophic level (figure 4). This is perhaps best illustrated by the sequential succession of small pelagic fishes demonstrating stable aggregate dynamics despite decadal-scale declines and near absences of individual species (figure 4c). We suggest that functional complementarity within trophic levels is the primary mechanism by which diversity maintains function and promotes resilience and stability in the SCCS. The functional complementarity is manifested by opposite responses of functionally similar species to positively autocorrelated drivers [20,24,36]. These drivers involve local or regional forcing, such as upwelling and temperature, impacting key population processes (e.g. growth, mortality and recruitment) or patterns of transport and advection, causing dispersal and range shifts of complementary species. The opposite response to these drivers leads to asynchronous population dynamics that serves to reduce variability and maintain functioning at the community level. Likewise, opposite responses to climate drivers have been shown for lower trophic levels, e.g. euphausiids [21,40] and copepods [50]. Functional complementarity may therefore underlie the long-term stability of total mesozooplankton biomass, despite ENSO and other perturbations (figure 4f). Although northern and southern mesopelagic fish respond differently to temperature and alongshore transport, the lower degree of compensatory dynamics is likely due to a similar response to deep-water oxygen [39,51]. Furthermore, the degree to which functional complementarity operates at the level of mammal and bird predators is unclear. The negative correlation between sea lions and total sea bird densities may indicate opposite responses to drivers, including changes in prey availability [25,26], yet further research applying a comparable degree of taxonomic resolution is needed to support this notion.

• Download figure
• Open in new tab
• Download PowerPoint

Over evolutionary timescales, competition has surely shaped niche partitioning and functional complementarity, e.g. manifested by different gill morphology, prey preferences and habitat use [20], but unlikely serves as the present cause of asynchrony. Hence, our results conform to theory showing the importance of complementarity, rather than direct competition in promoting stability [52]. Our results highlight a potential difference with terrestrial ecosystems. There, based primarily on experiments, compensatory dynamics are considered rare [53], whereas statistical averaging instead is thought to promote stability [49]. Although controlled experiments limit natural processes related to species composition and diversity (i.e. dispersal, invasion and colonization), the contrast is likely attributed to different variance spectra among terrestrial and marine environments [54]. In marine ecosystems, adaptations to more red-shifted drivers (i.e. higher variance at lower temporal frequencies) should favour compensatory dynamics among species. Furthermore, shorter generation times relative to the characteristic timescale of environmental variability affect the response to climate drivers [24,40,45]. The relative importance of statistical averaging and functional complementarity in promoting community stability and resilience across marine and terrestrial ecosystems merits further attention.

(c) Conclusions and management implications

Anthropogenic impacts, notably overexploitation and climate change, threaten important ecosystem functions and services [1], including provisioning, regulating, supporting and cultural services (electronic supplementary material, table S9). Therefore, maximizing resilience plays a central role in ecosystem management [1,55]. Although theory provides a conceptual basis for understanding resilience [12,17] and the positive effect of BEF [5–7], its applicability to management requires empirical insights regarding the underlying mechanisms. Using modelling and a unique long-term dataset, we have shown that the pronounced multi-decadal variability and abrupt transitions in the SCCS does not represent fundamental changes in ecosystem functioning, but a linear response to key environmental drivers channelled through bottom-up control. We show strong temporal asynchrony within multiple trophic levels and conclude that functional complementarity is the primary mechanism promoting resilience and long-term stability. Functional complementarity may therefore buffer against low functional redundancy (i.e. a limited number of species that exhibit similar ecological functions), which has been shown for many coastal and marine fish communities [56]. In addition to functional complementarity adaptive prey preferences and feeding modes are common in the SCCS [20,57–59]. By efficiently using a wide range of prey in a dynamic physical and trophic environment, adaptive foraging and omnivory may contribute to reducing short-term population variability and increase food-web stability [60,61]. Assessing species responses to key environmental drivers and mapping the degree of niche overlap within ecological communities will enhance our understanding of functional complementarity and the stabilizing role of diversity. More importantly, it may serve to identify species or functional groups that merit special protection, as well as guide appropriate management actions, in order to enhance ecosystem resilience and ensure the productivity and sustainability of ecosystems in the face of a changing climate.