Historical ecological datasets from a coastal marine community of crustose coralline algae (CCA) enabled the documentation of ecological changes in this community over 30 years in the Northeast Pacific. Data on competitive interactions obtained from field surveys showed concordance between the 1980s and 2013, yet also revealed a reduction in how strongly species interact. Here, we extend these empirical findings with a cellular automaton model to forecast ecological dynamics. Our model suggests the emergence of a new dominant competitor in a global change scenario, with a reduced role of herbivory pressure, or trophic control, in regulating competition among CCA. Ocean acidification, due to its energetic demands, may now instead play this role in mediating competitive interactions and thereby promote species diversity within this guild.
Changes in climate can alter species interactions by affecting individual responses to resource availability, including seasonal patterns [1–3]. Much work in this area has focused on species-level effects, such as latitudinal  or microhabitat variation, including a focus on intertidal marine systems . However, the scaling of species or population response to potential change in inter-species interactions has most often been studied in terrestrial systems in the context of phenology, or seasonal timing of life-history events. For example, trophic mismatch emerging from phenological shifts has disrupted plant–pollinator interactions as plant and insect phenologies respond disjointedly to an earlier and more variable start of the growing season [6–10].
Major changes in climate are also ongoing in the marine environment, many of them related to the process of ocean acidification. The relative abundances of dissolved carbon dioxide (CO2), bicarbonate (HCO3−) and carbonate (CO32–) in seawater are altered as high concentrations of atmospheric CO2 cause the oceans to take up increasing CO2 [11,12]. Because the relative concentrations of these three forms of inorganic carbon change with pH, increasing acidification (decreasing pH) favours bicarbonate over carbonate. Marine organisms rely on different forms of dissolved inorganic carbon for different physiological processes. Increased carbon dioxide and bicarbonate availability associated with ocean acidification may boost algal productivity [13,14] (but see [14,15]), however; lower carbonate saturation increases the seawater solubility of carbonate structures, including the shells and skeletons of many marine organisms [16–21] (but see [15,22]).
The coralline red algae, a group of marine macroalgae (Corallinales, Rhodophyta), make skeletons composed of high-Mg calcite, which is thought to be the mineral form most susceptible to ocean acidification [17,18]. Here, we focus on a guild of intertidal crustose coralline algae (CCA) that form a network of overgrowth interactions, the outcome of which is affected primarily by species morphology and growth rates [23–26]. Competition for space is important in structuring marine hard-bottom communities [27–29], typically occurring as interference competition between sessile plants and invertebrates [24,29].
CCA are an ideal system for the use of two-dimensional cellular automaton models. CCA encrust rocks and other hard substrates on the seafloor and intertidal regions. They reach only a few millimetres in height and grow radially outward in the absence of substrate irregularity and disturbances . CCA compete for space by means of planar overgrowth interactions, the outcome of which depends on the thallus thicknesses and lateral growth rates of the competing species . Further, the area of the thallus scales directly with photosynthetic productivity and reproductive investment, making two-dimensional size a predictor of fitness in some crusts .
The presence of mobile invertebrate grazers has historically been important in this system. Invertebrate grazers include chitons, limpets and urchins, who are able to rasp through the calcified coralline algal skeleton . Overall thallus thickness, which can protect CCA individuals from grazing wounds, exists in a trade-off with lateral growth rate due to the restraints of resource allocation [25,26]. When grazers are present, coralline crusts spend more energy on wound-healing mechanisms, and these conditions favour crusts that grow slowly both vertically and laterally over thin, fast-growing strategies . During the 1980s, experimentally increased grazer presence from 0 to 10 and 100% of its natural abundance was found to change the coralline algal competitive network from a pure hierarchy to an increasingly intransitive network . Here, we define a hierarchy when all species of higher rank have a high probability (arbitrary cut-off of p > 0.8) of overgrowing all species of a lower rank. By contrast, we define an intransitive network when at least one species of lower rank overgrows at least one species of a higher rank, causing the absence of a strict competitive hierarchy . For our purposes, intransitivity can therefore come in the form of a rank reversal or variation in competitive outcome between species such that a strict ranking or hierarchy is absent.
Ocean acidification experienced by natural CCA communities has also been shown to promote competitive intransitivity, both in the presence and in the absence of grazers . Observed changes in the coralline algal competitive dynamics at Tatoosh Island, WA, USA, have been linked to thallus thickness and growth strategy over a period of acidification [32–34]. Based on altered observations of ecological dynamics under ocean acidification, we parametrized a cellular automaton model with empirical data from historical and recent experiments [24,33]. Using historical (1980s) and modern (2010s) overgrowth frequencies as approximations for pre- and post-acidification states, we predicted patterns of species coexistence under the continued stress of ocean acidification. Specifically, we asked whether the current community is representative of a new equilibrium state, whether all species would remain in the community at equilibrium and what their relative abundances would be. We also asked whether trophic control remained strong when coupled with acidification stress, or if, alternatively, acidification has altered the control of producer community structure by consumers.
2. Material and methods
(a) Model overview
Simulations were run using a cellular automaton model (see electronic supplementary material, figure S4). The model was run using local interaction rules , parametrized using species abundances and species interaction strengths observed from field surveys and field experiments compared over a 30-year time span during which seawater pH has declined [24,33]. Whereas species abundances and competitive interactions have been well described, colonization and disturbance rates remain relatively unknown. Thus, the relative contributions of external and internal sources of colonization were estimated from historical species abundance patterns. Further details of the model are described below.
(b) Parametrization of overgrowth probabilities
Our model compared the outcome from experimental and survey data on overgrowth interactions among a guild of CCA that were collected in the 1980s  versus those quantified in the 2010s  at Hedophyllum Cove on the northeastern face of Tatoosh Island in Washington, USA (48.4° N, 128.7° W). The species included in this guild are the CCA Lithophyllum impressum, Lithothamnion phymatodeum, Pseudolithophyllum muricatum and Pseudolithophyllum whidbeyense. Experimental and survey data on competitive overgrowth consist of overgrowth win–lose frequency counts from previous studies (see electronic supplementary material, table S1). Data from the 1980s were collected as either surveys of overgrowth competition in unmodified communities with grazers present (historical scenario H1) or experimental transplants with grazers removed (H2) . Data to match these scenarios were collected in the 2010s (modern scenarios M1, M2) and additionally from experimental transplants with grazers present at natural abundance (M3) . Changes in pH isolated from other changes to the seawater environment have been documented at this field site over the 30-year interval [36–38], allowing us to compare the consequences of a change in competitive outcome through time, both with and without grazers.
(c) Model simulation procedure
Hard substrate within the intertidal zone was represented by a grid of 500 × 500 cells, with each containing only one species at a time. The 500 × 500 cell grid was determined to represent an area of 1.25 m2 using average growth rates from observations across all species , and a single-model time step is therefore equal to 1 year. To avoid edge effects in model dynamics, the grid was wrapped to form a torus by connecting opposite edges. At the start of a simulation, the empty grid, representing bare substrate, was populated at random with 50 individuals of each species placed at random on this two-dimensional system.
Next within a time step, a focal cell was chosen at random. If the focal cell was empty, colonization was possible. Colonization could occur either from the local species pool, determined by the abundance of each species on the simulation grid, or from an external propagule rain. The relative proportion of local to external colonization could be varied between simulations by assigning a value 0–1 to the parameter ProbLocalColonization, where a value of 0 represented 100% external colonization.
A random number 0–1, LocalOrMeta, was drawn from a uniform distribution and used to determine whether colonization is local or external to the focal cell within a given time step. If LocalOrMeta was less than ProbLocalColonization, then colonization proceeds from the local species pool. Under local colonization, the probability that any species will colonize the focal cell, ColonizeProb, was determined at random from an underlying distribution of the current abundance. Otherwise, if LocalOrMeta was greater than ProbLocalColonization, colonization was externally sourced according to the estimated colonization rate (see electronic supplementary material, figure S4).
If the focal cell was occupied at the start of the time step, then there was a possibility of disturbance acting on that cell. To simulate disturbance, a random number between 0 and 1 was drawn from a uniform distribution. This number was then compared with a set probability of disturbance, called the Disturbance Rate. If the random number was smaller than the Disturbance Rate, then the focal cell became empty. On the other hand, if the random number was larger than the Disturbance Rate, then the focal cell retained its original occupant.
If the focal cell was occupied at the start of the time step and disturbance had not occurred, then any of its eight neighbours (including diagonal) were identified. One of the neighbouring cells was chosen at random and its occupant species identified. To account for the fact that an ecological event did not necessarily occur at every time step everywhere on the substrate in the natural system, we specified a simulation procedure to determine whether or not overgrowth competition would take place between two given cells within a particular time step. The overgrowth outcome for a chosen neighbouring species and the species occupying the focal cell was determined from a random draw from a uniform distribution (0, 1), allowing the focal cell to defend itself from invasion. If the overgrowth probability was larger than the randomly generated probability, then the neighbour invaded the focal cell. Otherwise, the occupants of both cells remained unchanged.
Only one of these processes—colonization, disturbance or overgrowth competition—occurred at a focal cell within a single time step. Within a time step, each species within the grid was chosen as the focal cell at random, without replacement (i.e. each cell is the focal cell only once). Following each event (colonization, disturbance or competition), new species abundances were tallied. The simulation was run over 1000 time steps for each scenario, which was enough to generate equilibrial species abundances over many time steps.
(d) Parametrization of colonization and disturbance scenarios
Historical abundances in nature are a result of unknown external colonization and disturbance rates. To estimate these, we used a hill-climbing optimization procedure  to find best-fit rates using historical data from natural communities (1980s with grazers) on overgrowth probabilities and relative species abundances determined by surveys (table 1).
To optimize external colonization rates, we set the probability of internal colonization equal to zero (i.e. a fully externally sourced system). To begin, external colonization rates, expressed as probabilities between 0 and 1, were drawn at random from a uniform distribution. The cellular automaton was run using overgrowth parameters from scenario ‘H1’ for 100 time steps, which in each case was long enough to reach a stable species abundance distribution. The model-derived relative abundances of each species were compared with historical species abundances using the sum of differences between simulated and historical abundances of all five species. This number, referred to as model ‘fitness’, was saved. Subsequently, the hill-climbing algorithm chose a random increment between −0.05 and 0.05 by which to modify the colonization rate of a given species. This was repeated for all five species' colonization rates and for the disturbance rate. Disturbance was bounded between 0.01 and 0.5. Colonization rates of all species were unbounded.
The cellular automaton was re-run using the new colonization and disturbance rates. If the simulation fitness number was better than the saved fitness, it was saved along with the colonization and disturbance parameters used to generate the best-fit abundance. We repeated this hill-climbing algorithm for 1000 steps and determined the best overall fit and parameters. This entire procedure of hill-climbing for 1000 steps was repeated 1000 times to arrive at a best-fit solution for both external colonization and disturbance rates (see electronic supplementary material, table S2).
We parametrized disturbance and colonization rates by matching the output of historical model ‘H1’ (competition rules defined by historical surveys in the presence of grazers) to historical species abundance data, and applied these disturbance and colonization rates to all model scenarios. By doing so, we generated model outcomes that led to stable coexistence of all species, reflecting field conditions.
(e) Community metrics
To determine the effect of trophic interactions on community diversity, the model was run under a given set of conditions until a mature community developed, then re-parametrized using different overgrowth scenarios. For example, to simulate the addition of grazers to a mature grazer-free community, the output grid obtained from a grazer-free scenario (H2, historical grazer-free, or M2, modern grazer-free) was used as an input parameter to a new simulation, this time using overgrowth scenarios in the presence of grazers (H1 and M3, historical and modern, respectively). Similarly, to simulate the onset of ocean acidification scenarios, the output grid obtained from a historical simulation (H1) was used to seed an ocean acidification simulation (M3) (table 2).
For all model scenarios, we quantified the abundances of all five species in the system in terms of the number of cells each occupied. We assessed the diversity of the assemblage with Simpson's diversity index, calculated using the vegan package in the R statistical language (R v. 3.1.0 GUI v. 1.64 Mavericks build).
Comparison of model outputs allows us to detect linkages between changes in the environment and the response of the algal community. An effective way to do this is to simulate the community under a given scenario until it reaches a mature state, when simulated species abundances have become stable. When we subsequently apply altered competitive rules, we can interpret the dynamic response of community composition as a time-integrated response of the effects of changing those rules on the community (figure 1; for baseline simulated abundances, see electronic supplementary material, figure S2).
Historical (1980s) competitive networks showed a strong dependence on grazer presence for generation of intransitive competitive networks that led to species coexistence . Model scenarios H1 and H2 tested the historical effect of grazers, who have been shown to play a significant role in structuring the community (figure 1a,b) [24,25]. This comparison showed only a small effect of grazer addition on species abundances, with an increase in the abundance of L. phymatodeum and decrease in that of L. impressum. Similarly, comparing scenarios M1 and M2 reveals that grazer presence in the modern community also increases abundance of L. phymatodeum and decreases that of L. impressum (figure 1c). Unlike in the historical community, competitive networks respond with surprisingly little change in competitive networks and none in the degree of network intransitivity  (figure 1d).
Simulated changes to competitive networks documented over the past 30 years, comparing H1 to M3, led to large changes in competitive dominance and rapid corresponding shifts in long-term species abundances (figure 1e,f). The magnitude of this effect overshadowed that of grazer effects on species abundances in the modern community (figure 1c,d), largely driven by the rapid decline in abundance of the historical competitive dominant, P. muricatum (figure 1e). Notably, this rapid change occurred in the only comparison scenario that included large changes in species rank (figure 1f). Furthermore, empirical competitive networks obtained by competition experiments (M3) exposed a lag time in the natural community (M1) caused by memory effects in slow-growing, long-lived populations, such as these CCA . Model results were then used to provide insight on the possible transitional states of current CCA communities and to determine whether new species dominance and abundance patterns are likely to emerge as a result of altered competitive networks.
Output abundance patterns of the different model scenarios revealed the emergence of a new competitive dominant (L. phymatodeum) in the modern experimental scenarios, both with and without grazers (M3, M2) but not in the modern survey scenario (M1) where P. muricatum continued to dominate the community (see electronic supplementary material, figure S2). These model results matched the competitive dominant indicated in each corresponding competitive network. All scenarios led to stable coexistence of all five species.
The collapse of P. muricatum observed in modern scenarios was robust to the presence or absence of grazers. Adding modern interaction rules to historical baselines with (figure 2a–c) and without (figure 2d–f) grazers revealed the emergence of a new competitive dominant in each case. In both comparisons, abundances of L. impressum and L. phymatodeum increased. With grazers present, L. phymatodeum became the new dominant at both higher (50%) and low (0%) degrees of external colonization (figure 2a,b). With grazers absent, external colonization played a role in determining dominance in the modern scenario. When colonization was purely local, L. phymatodeum was the primary space occupier, though its abundance was only slightly higher than that of L. impressum. Lithophyllum impressum became the dominant space occupier when 50% external colonization was permitted (figure 2e), despite L. phymatodeum being predicted as the competitive dominant by interaction diagrams (figure 2f). Neither L. impressum nor L. phymatodeum dominated the community to the extent that P. muricatum did historically.
The permitted amount of external colonization to the system had the largest effect on historical community abundance, and smaller probabilities of external colonization extended the time it took for the community to respond to the onset of acidification (figure 2, comparing a with b and d with e). Decreasing the proportion of external colonization created monocultures, both in the presence and absence of grazers across all scenarios. Note also that exclusively local colonization did not permit establishment of P. muricatum as a dominant space occupier in historical simulations (figure 2a,c; electronic supplementary material, figure S2). Generally, the relative contribution of local colonization needed to be high to obtain a community that was dominated (more than 50% abundance) by one species in the modern, low-pH community. Within the modelled modern scenarios, the minimum ‘threshold’ amount of local colonization required for a community to become dominated by a single species was higher in the presence of ocean acidification stress (see electronic supplementary material, figure S3). Thus, stress promoted coexistence in this community because dominance required increased local colonization (or a more closed system, where colonization, is increasingly dependent on local abundance) as pH increased.
All scenarios had similar community diversity, with values of Simpson's diversity index near 0.8 when colonization was 100% externally sourced, and decreasing to 0.4 when colonization was strictly local (figure 3). The functional response of each scenario to reductions in external colonization differed, with modern experimental scenarios showing the highest robustness to reduced external colonization, and historical experiments and surveys showing the greatest sensitivity in species richness to reducing external propagule supply.
The strength of grazer mediation of CCA competitive interactions and promotion of intransitivity has decreased over time. In past studies, grazer presence and abundance have promoted coexistence within CCA by increasing the intransitivity of competitive networks, and maintained biodiversity by preventing or delaying the establishment of a competitive dominant [24,25]. Comparing the outputs of models parametrized with 2010s experimental data with (M3) and without (M2) grazers shows the effect of grazers in modern communities. The presence or absence of grazers did not change patterns of species abundance in modelled communities, even though the interaction networks have changed (see electronic supplementary material, figure S2). This may be due to the fact that while species' competitive rankings have changed slightly, the degree of intransitivity has not changed by the same magnitude as in historical scenarios, in which the removal of grazers led to an entirely hierarchical network . Network intransitivity is thus likely to be the primary driver of the differences in species abundance patterns as opposed to species rank.
We can also observe the effects of grazers on competitive dynamics by examining changes in species abundance patterns when a simulated mature community in the absence of grazers is subjected to the addition of grazers via changed competition rules (figure 1a–d). The transition over 30 years from historical to modern competitive networks generates the largest change in species abundance patterns (figure 1e,f). In this case, species placement in the relative hierarchy rather than overall network intransitivity may be driving these changes. This contrast with drivers of species abundance in response to grazer presence (M2 versus M3, figure 1c,d) stems from the relatively large displacement of species' competitive rankings between historical and modern scenarios, in which P. muricatum has been completely displaced (H1 versus M3, figure 1e,f).
Similarly, disturbance can act to promote coexistence by affecting all competitors equally [40,41], delaying the establishment of stable ecological dynamics and the eventual dominance of a competitively superior species . Disturbance may also act via compensatory mortality, where disturbance has a larger effect on the eventual competitive dominant [43,44]. Though our parametrized disturbance rate was relatively high (42.5%) to match observed 1980s abundance patterns, it was held constant over all model scenarios. Thus, disturbance probably contributed to the stability of our scenarios, but could not solely account for differences in diversity and abundance metrics between models because disturbance rates remained identical across scenarios.
Modern, low-pH conditions suggest the emergence of a new dominant species, L. phymatodeum, in experimental surveys and confirmed by our model (M3). Comparison of experimental and survey data as well as models parametrized by each dataset enable us to assess whether the modern natural community is in a state of transition. Neither surveyed species abundances nor modelled abundances parametrized by competition parameters obtained from census data (M1) reveal a reduced competitive rank of P. muricatum or the emergence of L. phymatodeum as the new dominant (see electronic supplementary material, figure S2), although both indicate increased species abundance of L. phymatodeum (see electronic supplementary material, figures S1 and S2).
Taken together, our models indicate that the current community will continue to respond to stress imposed by changes in seawater chemistry, although it has so far been slow to change. The slow dynamics observed in the non-experimental scenario (M1) corroborates lagged responses observed in other studies of response of long-lived, slow-growing plants to environmental changes [45,46]. Such lagged responses, or disequilibrium dynamics , are generated as an artefact of species demographic rates and not due to the magnitude of response to a particular stressor [46–48]. This effect is often exacerbated due to the prevalence of clones in long-lived plant communities, which increases genetic persistence . Indeed, such responses point to a need for both experimental and model approaches to understand change, even when long-term data are available [45,46,49].
One way to resolve the question of dynamics with respect to competition driven by slow growth rates is through experimental scenarios . In situ experimental transplants eliminated carry-over effects from previous or ongoing competition among species and were able to mediate changes to interaction strengths by effectively resetting competition to reflect current competitive abilities. Thus, M1 versus M3 comparisons indicate, first, that the modern natural community is not at equilibrium, taking the M3 model abundances to reflect stable conditions under current competitive rules (see electronic supplementary material, figure S2). Second, the current community has reacted more strongly to the reduced competitive dominance of P. muricatum than predicted by any of the model scenarios using empirical interaction strengths (table 1; see electronic supplementary material, figure S1). Modern survey abundances were averaged over two nearby sites to control for any potentially anomalous population dynamics at a single site. This should have provided a conservative estimate on changes in species abundances, yet points to a decrease in the population size of P. muricatum.
These conclusions highlight the importance of interaction strengths in determining the effects of climate change on communities. Effects of climate change on community function and biodiversity are expected to scale up from species-level responses. However, a recent meta-analysis reveals that species richness and biodiversity metrics can remain unchanged even when species-level patterns indicate fundamental community reorganization . In our model CCA system, although L. phymatodeum becomes the dominant in this guild, it does not reach the same abundance in the community that P. muricatum did historically, allowing competitively inferior species to coexist at slightly higher abundances than in historical scenarios. These findings contrast with those from other studies that predict ecosystem simplification under ocean acidification stress at natural CO2 vent systems [51,52]. Whether this is due to the extreme carbonate changes at vents compared with the relatively small changes in this system remains an area for study. Our results, however, suggest that species interactions will be the underlying driver for biodiversity changes as ocean acidification continues [53,54].
In both historical and modern scenarios, the persistence of the competitive dominant was enhanced when colonization was predominantly local, a common structuring feature of benthic marine communities . In our community, this happens as the dominant species takes over more space and increases its reproductive output, which scales with two-dimensional size. This suggests that environmental drivers that influence reproductive potential (and thus recruitment) have been important in structuring CCA communities. In modern scenarios, the minimum ‘threshold’ amount of local colonization required for a single species to become dominant in terms of space occupancy becomes higher, meaning that species coexistence is facilitated except in the case of a nearly closed system (see electronic supplementary material, figure S3). Going forward, our models thus indicate that environmental drivers of algal reproductive potential may become less important in structuring future CCA communities.
Despite empirical and model evidence of altered community dynamics within this guild, the prognosis for maintenance of local coralline algal diversity is uncertain. When dynamics are followed to equilibrium, current levels of ocean acidification may act as a form of intermediate disturbance in this system. By this mechanism, acidification stress may be playing the role traditionally held by grazers in this system, that of reducing the competitive dominance of a single species and promoting local biodiversity through compensatory effects and competitive release. However, if acidification intensifies in the future ocean, it may have stronger effects on all competitors and drive down population sizes and diversity. Further, it is not known whether communities will actually reach the equilibrium abundances predicted by our model scenarios. If the differences currently observed between natural community surveys and model outputs persist due to additional physiological stress responses of coralline algae, acidification may have already become too stressful to allow persistence of all species.
Other recent work has identified a negative indirect effect of acidification on calcified grazers via reduced nutritional quality of algal food sources rather than a direct detrimental effect of acidification on grazers . It seems therefore that ocean acidification could increase bottom-up forcing, with effects on algae working up the food web, as well as decrease trophic forcing as we have found in our study. On the other hand, trophic interactions may play an increasingly important role in buffering effects of environmental disturbances that may propagate up the food web . The structuring role of trophic interactions in the future ocean may then depend on several factors in the face of climate change , including the continued absence of phenological shifts in macroalgae and grazers, minimal direct negative effects of acidification on grazers, and the relative importance of bottom-up versus top-down trophic control in different coastal systems.
The datasets supporting this article are all previously published and available in the supplementary material of McCoy & Pfister . Model code and input parameters have been uploaded as part of the electronic supplementary material.
S.J.M. designed the study, conducted modelling work and data analysis, and drafted the manuscript; S.A. conducted modelling work and advised data analysis, and helped draft the manuscript; C.A.P. advised study design and data analysis, and helped draft the manuscript. All authors gave final approval for publication.
We have no competing interests.
Support came from the United States National Science Foundation Doctoral Dissertation Improvement Grant NSF DDIG DEB-1110412 (C.A.P. and S.J.M.), the ARCS Foundation (S.J.M.), NSF-OCE 09-28232 (C.A.P.) and NSF 11-48867 (S.A.). This research was conducted with United States Government support under and awarded by DoD, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a and by the National Science Foundation Graduate Research Fellowship (GRFP-1144082). S.J.M. is currently supported by a Marie Curie International Incoming Fellowship within the 7th European Community Framework Programme (grant agreement FP7-PEOPLE-2012-IIF No. 330271).
We thank the Makah Tribal Nation for access to Tatoosh Island and R. T. Paine for thought-provoking discussions and the availability of historical data. Many thanks to J. T. Wootton, J. E. Cohen and E. Sander for comments and suggestions on the model structure, and to I. Miller for computing support.
- Received October 23, 2015.
- Accepted February 8, 2016.
- © 2016 The Author(s)
Published by the Royal Society. All rights reserved.