Species fluctuations sustained by a cyclic succession at the edge of chaos

Abstract

Although mathematical models and laboratory experiments have shown that species interactions can generate chaos, field evidence of chaos in natural ecosystems is rare. We report on a pristine rocky intertidal community located in one of the world's oldest marine reserves that has displayed a complex cyclic succession for more than 20 y. Bare rock was colonized by barnacles and crustose algae, they were overgrown by mussels, and the subsequent detachment of the mussels returned bare rock again. These processes generated irregular species fluctuations, such that the species coexisted over many generations without ever approaching a stable equilibrium state. Analysis of the species fluctuations revealed a dominant periodicity of about 2 y, a global Lyapunov exponent statistically indistinguishable from zero, and local Lyapunov exponents that alternated systematically between negative and positive values. This pattern indicates that the community moved back and forth between stabilizing and chaotic dynamics during the cyclic succession. The results are supported by a patch-occupancy model predicting similar patterns when the species interactions were exposed to seasonal variation. Our findings show that natural ecosystems can sustain continued changes in species abundances and that seasonal forcing may push these nonequilibrium dynamics to the edge of chaos.

Keywords: chaos; coexistence; cyclic succession; rocky intertidal community; rock–paper–scissors dynamics.

Fig. 1.

The rocky intertidal community. (A) Aerial view of the study site at Goat Island Bay. (B) Cyclic succession at the rocky intertidal site. First, barnacles settle on bare rock, and second, crustose algae invade. Third, mussels settle on top of the barnacles and crustose algae, forming a dense carpet that smothers the barnacles and algae underneath. Fourth, the mussels detach, bare rock becomes available again, and the cycle restarts. Drawn by Jan van Arkel (University of Amsterdam, Amsterdam, The Netherlands). (C) Time series were obtained from a permanent grid consisting of 20 plots and five nodes (A–E). The percentages of cover of barnacles, crustose algae, and bare rock were monitored in the plots, whereas mussel cover was estimated from photographs of the nodes. Ten plots and five nodes within the red line were used for the time series analysis.

Fig. 2.

Time series observations of (A) sea surface temperature, (B) barnacles, (C) crustose algae, (D) mussels, and (E) bare rock.

Fig. 3.

Cross-wavelet spectra of all species pairs. Cross-wavelet spectra of (A–F) the observed time series (Fig. 2 B–E) and (G–L) the model predictions (Fig. 5 H–K). The spectra show how common periodicities in the fluctuations of two species (y axis) change over time (x axis). Color indicates cross-wavelet power (from low power in blue to high power in red), which measures to what extent the fluctuations of the two species are related. Black contour lines enclose significant regions, with >95% confidence that cross-wavelet power exceeds red noise. Arrows indicate phase angles between the fluctuations of the two species. Arrows pointing right represent in-phase oscillations (0°), and arrows pointing upward indicate that the first species lags the second species by a quarter period (90°). Shaded areas on both sides represent the cone of influence, where edge effects may distort the results.

Fig. 4.

Time series of the local Lyapunov exponent (LLE). The LLE is calculated as the average rate of trajectory divergence (or convergence) over a time span of 180 d.

Fig. 5.

Model simulations with the patch-occupancy model. (A–E) Without seasonal forcing, the species show damped oscillations that eventually settle down to equilibrium, and (F) the attractor is a stable spiral point. (G–K) With seasonal forcing, the species show sustained irregular oscillations characterized by (L) a strange attractor typical of chaotic systems. For the seasonally forced simulations, the first 30 y were discarded to remove initial transients. The chaotic attractor in L was obtained by running the model for 100 y. Parameter values are provided in SI Appendix, Table S1.