Impact of the extinct megaherbivore Steller's sea cow (Hydrodamalis gigas) on kelp forest resilience

1. Introduction

The time frame of concern for most conservation projects is limited to a temporal window spanning less than a century (Dietl and Flessa, 2011). When contemplating interventions, conservation efforts tend to focus on the recent past and near future. This strategy can be problematic in a world where ecosystems have been changing for centuries to millennia, climate is warming rapidly, and ecological baselines are shifting (Dayton et al., 1998), making the recent past an unreliable guide to conditions preceding anthropogenic impact. Understanding what diversity and resilience once looked like for a given ecosystem is critical if we aim to regenerate ecosystems capable of persistence in an uncertain future. Rigorous exploration of future impacts from climate warming and other factors, including sequential reintroduction of species, or rewilding, is essential to increase the probability of regenerating past states that may have been more resilient. Here, we propose a Past-Present-Future (PPF) approach to conservation rooted in mathematical modeling, illustrating this method with a study of the giant kelp ecosystem in the northeastern Pacific Ocean.

Giant kelp forests are one of the most productive marine ecosystems in the world (Foster et al., 2013), distributed in cool temperate coastal regions of both hemispheres (Steneck et al., 2002; Graham et al., 2007). Northern Pacific kelp forests are particularly iconic, and studies on these ecosystems have become fundamental to the concepts of keystone predators, foundation species, and alternative ecological states (Estes and Palmisano, 1974; Filbee-Dexter and Scheibling, 2014; Castorani et al., 2018). Many kelp-dominated communities include sea urchins as major kelp grazers. The kelp-urchin system typically exhibits two states; forest and barrens (Mann, 1977; Filbee-Dexter and Scheibling, 2014). Various circumstances can lead to uncontrolled urchin grazing, resulting in a landscape denuded of kelp. The productivity and biodiversity of such communities are so drastically reduced that the community state is referred to as an “urchin barrens” (Graham, 2004). The transition from a forest to barrens is well-understood in North Pacific forests: predatory sunflower sea stars (Pycnopodia helianthoides) and sea otters (Enhydra lutris) control sea urchin populations, facilitating the development of dense kelp forests dominated by the species Macrocystis pyrifera and Nereocystis luetkeana. Reduction of either predatory species may drive a state transition from forest to barrens. While hydrodynamic forces are an important environmental control on kelp forest biomass in certain contexts (Reed et al., 2011; Bell et al., 2015) and are included in the model (Supplementary Section 2.1.3), here we focus on evaluating how changing ecological interactions dictate kelp forest system state.

Widespread anthropogenic extirpation of the otter resulted in the conversion of numerous forests to barrens in the nineteenth and twentieth centuries. More recently, marine heat waves, and the emergence of sea star wasting disease (SSWD) in the North Pacific devastated populations of P. helianthoides, facilitating another widespread transformation of forests to barrens (Hamilton et al., 2021; McPherson et al., 2021).

1.1. The North Pacific kelp-urchin system

Kelp forests consist of dense concentrations of giant kelp that grow from holdfasts anchored to hard substrate at depths up to 15 m. Robust stalks extend from the substrate to the surface, bearing blade-like fronds that capture light. In M. pyrifera-dominated forests, blades form a dense canopy, capturing up to 90% of incident sunlight (Wing et al., 1993; Detmer et al., 2021). Forests have diverse assemblages of other organisms, high overall primary productivity, baffle wave energy and storm surges, and are conventionally considered to be the “healthy” state (Figure 1A). Adult kelp, however, do inhibit other photosynthesizers—phytoplankton, benthic algae, and juvenile kelp (Graham et al., 1997)—because of their high light absorption and competition for substrate space. Barrens, in contrast, consist of hard substrate dominated numerically by urchins, and largely devoid of kelp or understory algae (Graham, 2004) (Figure 1B).

FIGURE 1

Figure 1. Observed and hypothetical Macrocystis pyrifera forest states. (A) Modern kelp forest state showing a dense stand of adult kelp, sparse understory algae, purple sea urchins, sunflower sea stars and a sea otter. (B) The urchin barrens state, with no kelp, a greater abundance of sea urchins, sparser understory algae, and no sea otters. (C) A hypothetical historical state, similar to (A), but with Steller's sea cow present, a more open kelp canopy, and more abundant understory algae.

Transition of a forest to barrens is often abrupt relative to the duration of the forest, and preceded by an eruption of actively foraging urchins. Typically, a large fraction of urchins live crowded into substrate crevices when predators are present, capturing and consuming the abundant drift detritus produced by the fragmentation of kelp and other macroalgal species. Two factors can increase the fraction of foraging urchins: a reduction of predation pressure or a decline in the availability of drift detritus. Lower predation pressure results in growth of the urchin population and a reduction of available crevice space and food availability, while a reduction of drift detritus leads directly to starvation, which may determine if urchins remain sheltered, or instead become mobile foragers (Kriegisch et al., 2019). In both cases, urchins will begin to leave crevices to graze directly on kelp, resulting in the death of individual kelp. This may create a positive feedback whereby the loss of adult kelp further reduces the availability of drift detritus, prompting more urchins to become mobile foragers (Karatayev et al., 2021). Furthermore, as kelp density decreases, urchin starvation increases, leading to a decline in the nutritive condition of the urchins. Otters generally avoid predation on less nutritive urchins (Smith et al., 2021), exacerbating and accelerating the transition from forest to barrens. Barrens can transition back to forest if urchin populations decline sufficiently because of starvation. Feehan et al. (2018) have shown that offshore export of particulate kelp detritus is an important food source for urchin larvae. Hence the devastation of kelp productivity by adult urchins could generate a negative feedback on the urchin population. Nevertheless, urchins appear capable of subsisting at low tissue mass for long periods of time under low nutrient conditions, a condition termed “zombie urchins” (Spindel et al., 2021). Transition to forest can also be facilitated by an increase of sunflower sea star populations, or by high recruitment of juvenile kelp (Williams et al., 2021).

1.2. Current state

The dominant state of giant kelp communities in the northern Pacific over the past several decades has been maintained by healthy populations of sunflower sea stars, and the return of stable sea otter populations in some areas of their former range. The situation changed dramatically in 2013–2014 with the outbreak of SSWD on both coasts of North America, and a catastrophic shift to a reduced kelp canopy and dominance of urchin barrens communities (Rogers-Bennett and Catton, 2019). The disease quickly devastated sunflower sea star populations, with the species disappearing entirely in large parts of its range. An additional factor in the decline of the forests may have been the onshore expansion of a persistent body of warm water, colloquially known as “the blob”, around the same time. This persistently cohesive and anomalously warm body of water, first detected in the sub-polar northern Pacific in late 2013, expanded in 2014 to encompass almost the entire west coast of North America, with sea surface temperatures reaching 2.5–3.0^oC above average. The marine heat wave (MHW) persisted from 2013 to 2016 in some coastal areas (Mogen et al., 2022). It has been speculated that the warm temperatures had a detrimental effect on kelp density and productivity, with greater severity in northern California compared to central and southern areas (Rogers-Bennett and Catton, 2019). Certainly the coincidence of the MHW, SSWD, and the possibility of synergistic effects cannot be overlooked (Burt et al., 2018; Rogers-Bennett and Catton, 2019).

The heat wave and rapid decline of sea stars (Hamilton et al., 2021) were followed by widespread “outbreaks” of urchins, with increases in the numbers of actively foraging urchins and apparently local population sizes. These outbreaks resulted in the transition of numerous kelp forests to urchin barrens; for example, bull kelp (Nereocystis luetkana) was reduced by more than 90% (Rogers-Bennett and Catton, 2019). In Monterey Bay, the hypothetical setting for this modeling exercise, the outbreak of S. purpuratus that began in 2014 resulted in significant losses of kelp, and while there is evidence of recovery, that recovery has been dominated by the bull kelp, and not the formerly dominant giant kelp, M. pyrifera (Smith and Tinker, 2022). The loss of kelp productivity and alteration of the habitat has had cascading effects throughout forest communities; abalone (Haliotis) populations were reduced to 20% by 2017 (Rogers-Bennett and Catton, 2019). That state has persisted up to the present, although some recovery of kelp has been noted since summer 2021 along the coast of California (Bland, 2021). To date, the sunflower sea star remains functionally extinct from Baja California to Cape Flattery, Washington, and has disappeared entirely from large parts of its range (Hamilton et al., 2021).

1.3. A paleobiological alternative

The heavily studied giant kelp ecosystem may be historically young in parts of the North Pacific, at least in its present state. A megaherbivorous mammal, the Steller's sea cow (Hydrodamalis gigas), inhabited kelp ecosystems in the Comander Islands in historical times. This giant sirenian was first described by Georg Wilhelm Steller, a scientist who accompanied Russian commercial voyages to the islands in 1741. The species was quickly exploited as a source of fresh meat and is believed to have become extinct by 1768. Individuals attained very large sizes, up to 9 m in length, and are estimated to have weighed up to 10 tons. H. gigas was an obligate kelp feeder, apparently incapable of submerging as observed by Steller (Steller, 1751) because of its high buoyancy, and therefore browsed the giant kelp canopy available at the surface (Estes et al., 2016; Bullen et al., 2021). There are no recorded observations of the species beyond the Aleutian Islands, but fossils of H. gigas, and other members of the genus that might represent different species, are known from the Pliocene of California and Baja California, and the post-Pliocene of Japan (Domning, 1976). The genus may have evolved opportunistically with the late Neogene regional expansion of kelp forests. The rarity (34 partial skeletons have been recorded) of fossil remains makes it impossible to discern whether the genus once encompassed this entire range, including the Arctic during interglacial intervals, or if instead the range was fragmentary, responding perhaps to variable oceanographic conditions after the onset of northern Hemisphere cooling (~2.4 mya). Similarly, previous efforts to split the fossils into a phyletic lineage of multiple species is questionable (Domning, 1976), particularly given the sparsity of specimens and reliance on overall body size as a primary differentiating trait among the putative species. Nevertheless, all occurrences of Hydrodamalis, fossil and historical, imply a dependency on high productivity coastal habitats, and given body size and the angle of the snout, a feeding habit of surface browsing (Domning, 1976; Estes et al., 2016). Several hypotheses have been proposed to explain rarity, geographic restriction, and extinction of H. gigas, but none are conclusive (reviewed in Supplementary Section 1).

The impacts of megaherbivores on community structure and ecosystem functioning have been well-established (Owen-Smith, 1988; Gill, 2014; Hyvarinen et al., 2021), and their conservation is currently an issue of major concern. Given our understanding that these animals often contribute significantly to the dynamics of their communities, including stability, resilience, and state transitions, it is probable that H. gigas (and congeneric species) played a significant ecological role in kelp forest communities. Recently, Bullen et al. (2021) presented six hypotheses of how H. gigas may have structured giant kelp communities. Several of the hypotheses predict that grazing of kelp fronds by sea cows opened up the canopy, increasing light intensity in both the water column and on the seafloor (Figure 1C). One result would be increased production of non-kelp primary producers, including a higher biomass of understory algae, because the kelp canopy is known to suppress the abundances of understory algal species (Reed and Foster, 1984; Detmer et al., 2021). The authors also suggest that sea otter-urchin dynamics could have been altered if those conditions led to an increase in invertebrate prey diversity and biomass.

1.4. Ecosystem function of H. gigas—A new hypothesis

Here, we present and test an additional hypothesis of the functional role of H. gigas—that sea cow browsing would have increased the biomass of understory algae, and the concentration of drift detritus, thereby increasing the resilience of kelp forests. We define resilience as the capacity of a system, upon being perturbed, to either remain unchanged, or to recover to its state prior to disturbance. The particular case in which a system remains unchanged when perturbed is referred to as resistant. To test the hypothesis, we develop a mathematical model that describes the dynamics of a hypothetical kelp forest community off the central California coast, comprising interactions among the giant kelp (M. pyrifera), purple sea urchins (S. purpuratus), an understory alga favored by urchins (Chondrocanthus corymbiferus), sea otters (E. lutris), sunflower sea stars (P. helianthoides), and Steller's sea cow (H. gigas; Figure 2). We put forth two classes of models, referred to as “Historical” and “Modern”, with the former including the sea cow and the latter lacking the megaherbivore.

FIGURE 2

Figure 2. The complete model system, illustrating system components and interactions. Solid, red arrows indicate a positive or reinforcing impact in the direction of the link, e.g., kelp fronds have a positive impact on sea cow populations. Blue, broken arrows conversely indicate a negative, or dampening impact, e.g., kelp fronds have a negative impact on understory algae. The overall impact of a chain or loop is indicated as the product of the signs (positive or negative) of the links along the path. E.g., Steller's sea cow had a positive impact on understory algae, calculated as the product of the two negative links on the path Steller's sea cow-kelp fronds-understory algae (Organismal silhouette's courtesy of Phylopic, with the exception of the sunflower sea star. The giant kelp is attributed to H. N. Eyster.).

Although kelp density is controlled by numerous factors (Foster et al., 2013), for the purposes of the model we distill these to dominant factors, including temperature, light, hydrodynamic forces, spatial competition with understory algae, browsing of fronds by Steller's sea cows, and grazing by purple sea urchins. Understory algal density is controlled by spatial competition with kelp on the substrate, light attenuation by the kelp canopy, and sea urchin grazing. Sea urchin density in the model is controlled by sea otter and sunflower sea star predation, and the availability of drift detritus from kelp and understory algae. Browsing by sea cows is added to this description of the natural system to build the Historical model.

2. Methods

2.1. Model system

The Historical and Modern models of central Californian kelp forest are systems of biotic interactions (Figure 2) and physical parameters. Parameters are values that are fixed, or whose variation is independent of other model components, for example incident sunlight or storm events. Variables are partially or fully determined by interactions or dependence on parameters or other variables, for example species populations sizes or densities. The goals of the models developed here are: (1) to describe ecological interactions in a manner that results in a feasible forest state (i.e., with species persisting indefinitely) when the system is unperturbed; and (2) to predict system changes during and after perturbation.

The following sections outline the interactions for each model community. Complete derivations are available in the Supplementary material. Model communities are described with a system of coupled ordinary differential equations, which are simulated and solved using the Julia language. Code used to conduct simulations may be found at the Open Science Framework, URL osf.io/zvnw9. The modeling of physical parameters includes ocean temperature, day length and hydrodynamic forces, and is outlined in Supplementary material. Parameterization of biological parameters is also explained in Supplementary material. Definitions of all symbols in the following section are listed in Supplementary Table 1.

2.1.1. Species dynamics and biotic interactions

2.1.1.1. Kelp

M. pyrifera population dynamics are modeled as a function of an intrinsic rate of population increase, and predation by urchins and sea cows (when present). Population growth is dependent on both temperature and spatial competition with algae, where kelp are a superior competitor because of their greater access to light (Detmer et al., 2021). Population growth rate also accounts for a steady rate of extrinsic juvenile recruitment, but recruitment from external sources is dependent on spatial competition with resident adult kelp and understory algae, as well as available light at the substrate surface. Urchin predation is a source of direct kelp mortality, whereas sea cow browsing of fronds affects individual productivity and hence the population growth rate, but not mortality. Finally, storm surges are a seasonally stochastic source of mortality.

Complete kelp dynamics is expressed as

$\begin{array}{l}\begin{array}{lll}\frac{dM}{dt}\hfill & =\hfill & M\left[\right[\frac{r_{M}F}{20}(1-\frac{M}{2-C})]-[({\omega }_1{U}_E)+{\delta }_M]]\hfill \end{array}\\ \begin{array}{l}+\left[{\sigma }_M\right(1-\frac{M+C}{2})(1-0.03F)]\hfill \end{array}& (1)\end{array}$

where r_M is the rate of intrinsic increase of kelp, M and F are kelp and frond density, respectively, C is understory algal density, U_E is the fraction of the sea urchin population that is exposed and actively foraging, and ω₁ is the kelp-urchin prey-predator interaction. (1−M/(2−C)) is a logistic limitation of kelp density reflecting exploitative competition with understory algae C. σ_M is a temperature-dependent rate of external recruitment, where recruits are limited by both spatial competition with understory algae and adult kelp, and shading by adult kelp ((1 − 0.003F)). U_E is the density of actively foraging urchins (Supplementary Equation 14). ω₁ is the predator-prey interaction between kelp and urchins. This and all subsequent such interactions are ratio-dependent Arditi-Ginzburg-Contois (AGC) interactions (Arditi and Ginzburg, 2012) (Supplementary Section 2.1.4.1). δ_M is wave-induced mortality (see Supplementary Section 2.1.3). The constants present in the equation represent the maximum number of fronds per kelp, 20, and the shading rate of frond density, 0.003.

Temperature-dependent parameters, including rates of intrinsic increase and σ_M, have asymmetric, convex, relationships with temperature, with a single maximum at an optimum temperature (T_opt; Supplementary Table 1). Here, we modeled those relationships for individual species parameters using a Linex function (Supplementary Equations 1–3).

2.1.1.2. Kelp canopy

The canopy of kelp fronds at the water's surface plays several important roles in the community. The canopy absorbs a majority of incoming sunlight and inhibits primary production in the water column and benthos, including that of understory algae. Fronds are one of the major sources of algal macro-detritus, hence a food resource for urchins, and they were apparently the sole food for H. gigas.

Fronds are consumed by sea cows when they are present in the system, which in turn creates a negative feedback to kelp growth rate (dM/dt). Fronds also have a natural senescence rate at which they age and detach from the stalk. Thus,

$\begin{array}{ll}\frac{dF}{dt}=r_{F}M\left(1-\frac{F}{20M}\right)-F({\omega }_{5}H+{\delta }_F)& (2)\end{array}$

where r_F is frond production rate (implicitly assumed to equal 1), F is frond density, H is sea cow density, ω₅ is the interaction between fronds and sea cows of the AGC form (see Supplementary Equation 9), and δ_F is the frond senescence rate.

2.1.1.3. Drift detritus

Both kelp and the understory algae in the model contribute to drifting macrodetritus. Detritus is produced as a fixed fraction of frond and understory algal density, and becomes refractory (inedible) also at a fixed rate. Thus,

$\begin{array}{ll}\frac{dG}{dt}=0.8(F+C)-0.15G& (3)\end{array}$

where G is detrital abundance, the constant 0.8d⁻¹ is derived from empirical measures of detritial production by another kelp species, Laminaria hyperborea (Pedersen et al., 2020), and 0.15d⁻¹ is the offshore export and refraction rate. This is a donor-driven function in the sense that there is no impact of sea urchin consumption on detrital density in the model.

2.1.1.4. Understory algae

The red alga Chondrocanthus corymbiferus has been shown to be nutritious for urchins, readily consumed when available, and capable of sustaining urchins after kelp density has been reduced by large storm events (Foster et al., 2015). The alga contributes to the detrital pool, competes spatially with kelp on the substrate, is light-inhibited by the kelp, and is grazed by actively foraging urchins. Like kelp, its rate of intrinsic increase is temperature-dependent, with T_opt set here at 13^oC, slightly cooler than that of the kelp, based on the more limited southern extent of its range. Dynamics are given by

$\begin{array}{ll}\frac{dC}{dt}=C\left[\left[r_C(1-0.03F)\left(1-\frac{C}{2-M}\right)\right]-{\omega }_2{U}_E\right]& (4)\end{array}$

where r_C is the rate of intrinsic increase, and (0.03F) is the degree to which population growth is impacted by kelp frond density. ω₂ is the interaction coefficient with exposed, actively foraging urchins. Understory algal abundance is capable of doubling in the absence of kelp. While the carrying capacities of each species are equal to one in the model in the presence of the other species, those bounds can be exceeded if the competing species declines. Also, as a consequence, total kelp and understory algal abundances in the current model are always approximately equal to two (varying because of hydrodynamically-driven stochasticity).

2.1.1.5. Sea urchins

Growth of the total urchin population is determined by the foraging activity of exposed urchins, detritivory by both sheltered and exposed urchins, and predation by both major predators of the urchins: sea stars and sea otters.

$\begin{array}{ll}\frac{dU}{dt}={\epsilon }_1\left[{\omega }_{3}UG+U_E({\omega }_{1}M+{\omega }_{2}C)\right]-U({\omega }_{4}P+{\omega }_{5}E+{\delta }_1)& (5)\end{array}$

where ε₁ is a temperature-dependent ecological efficiency of the urchins of the same form as Equation S1, ω_n are predator-prey coefficients, P is sea star density, E is sea otter density, and δ₁ is the urchin mortality rate.

2.1.1.6. Sunflower stars

Pycnopodia helianthoides, like the urchins, is assumed to have temperature-dependent rates of prey interaction and ecological efficiency. Predation of sheltered vs. exposed urchins are tracked separately, and the interaction between sunflower sea stars and urchins is therefore expressed as

$\begin{array}{ll}\frac{dP}{dt}=\{\begin{array}{ll}P({ϵ}_2{\omega }_{4}U-{\delta }_2)\hfill & \text{if}t\in \text{yearsofnoperturbation}\hfill \\ -0.0013P\hfill & \text{otherwise}\hfill \end{array}& (6)\end{array}$

where ε₂ is the ecological efficiency, ω₄ is the predator-prey interaction term, and δ₂ is the mortality rate. ε₂ is the ecological efficiency, ω₄ is the predator-prey interaction term, and δ₂ is sea star mortality rate. Here, perturbation refers specifically to the effect of disease outbreak on sea star populations, which was implemented as an exponential decline of the P. helianthoides population, modeled from empirical data of population decline from central and northern California (Hamilton et al., 2021).

2.1.1.7. Steller's sea cow

Hydrodamalis gigas was apparently an obligate canopy feeder, based on observations of its feeding and its buoyancy. No instances of non-human predation were ever observed, although anecdotal behavioral reports suggest that the animals were wary of orcas and swam in formations that may have shielded juveniles (Marsh et al., 2012). It is also reasonable to speculate that large sharks, particularly great whites (Carcharodon carcharias), may have preyed on juveniles. Nevertheless, no predation is included in the model. H. gigas dynamics were therefore controlled by frond availability and a natural mortality rate.

$\begin{array}{ll}\frac{dH}{dt}=H({\epsilon }_3{\omega }_{5}F-{\delta }_3)& (7)\end{array}$

where ε₃ is the ecological efficiency with which consumed kelp was “converted” to new sea cow individuals. δ₃ is the natural mortality rate.

2.1.1.8. Sea otters

Sea otters, Enhyda lutris, are present only when kelp are sufficiently abundant and urchins are healthy in the model. Otter population dynamics are ultimately kelp-dependent, dependent on urchin nutritional state. Given the ephemerality of sea otters in our model community, we assumed that the otters are drawn from a fixed metapopulation pool. Further assuming that urchin gonadal index is a function of kelp density, and thus that the presence of otters preying on urchins is likewise dependent on kelp density, we used published data of the relationship between gonadal index and otter predation (Smith et al., 2021) to derive a logistic relationship between kelp density and the probability that otters are present and foraging (Supplementary Section 2.2).

$\begin{array}{ll}E=0.934{\left[1+e^{-0.2057(40M-11.31)}\right]}^{-1}& (8)\end{array}$

where E is sea otter density in the community, and M is kelp density.

2.2. Perturbation

Model communities were subjected to four types of perturbation: no perturbation; a MHW to simulate the warm Pacific Ocean blob; an outbreak of SSWD; and coincident warming and disease, SSWD+MHW. Model communities were constructed as described above and simulated for a burn-in period of 200,000 days, at which point the community was assessed to be in one of two states, forest or quasiperiodic, the latter state describing oscillations between the forest and barren states (see Results). Perturbations were then applied to the models.

In the case of no perturbation, the model was simulated for an additional 14,600 days (40 years) post burn-in. Perturbations were initiated precisely 10 years after the burn-in (day 203,650). The MHW consisted of an abrupt increase of temperature 3.0^oC in excess of the seasonal cycle, and lasted for 2 years (730 days), after which temperatures returned to the seasonal cycle. SSWD was also initiated 10 years after burn-in (see Equation 6), after which the sea star population was allowed to recover. Models were therefore simulated for 30 and 37 years after the termination of disease and warming, respectively, as the MHW persisted for 7 years less than SSWD.

A total of 100 simulations were performed for each of 16 total model-perturbation combinations. This included the Historical and Modern models (i.e., with and without sea cows) with either no perturbation, warming, disease, or warming plus disease, each with and without sea otters. This resulted in a total of 1,600 simulations.

2.2.1. Perturbed forests

Perturbed simulations were classified according to the scheme outlined for unperturbed simulations. Only simulations that were in the forest state at the end of the burn-in period were analyzed, as we are interested in the impact of the perturbations on forests. Those simulations could then either remain in the forest state after initiation of a perturbation, or transition to a barrens. The relative frequency of each type of response was recorded for each model-perturbation combination, of which there are 12.

2.2.2. Measuring resilience

Resilient simulations were defined as those that return to a forest state within 10 years after the end of a perturbation. The frequency of resilience was recorded for each set of simulations of a model-perturbation combination, along with the corresponding frequency of simulations that transitioned permanently away from the forest state.

Resilient simulations were further sub-categorized according to whether they exhibited no state transitions, or underwent one or more transitions but subsequently recovered permanently (within the time frame of the simulation) to the forest state. Simulations that did not undergo any transitions were categorized as resistant (Grimm and Wissel, 1997; Walker et al., 2004). Resilience was described with the time of onset of the first transition to barrens, and the time of the last transition back to forest. These two measures encompass one classic definition of ecological resilience, which is recovery time (Pimm and Lawton, 1977).

3. Results

3.1. Community states

Three states emerge from the unperturbed models, two of which have dense kelp forests. All species persist in the first forest state, corresponding to the conventional concept of a “healthy”, high diversity kelp forest (Figure 3A); we refer to this state hereon as simply “forest”. The other forested state is one in which S. purpuratus and P. helianthoides fail to persist, or do so at densities less than 1% relative to their initial values. We term this hereon specifically as “low diversity forest”. The third state is a quasiperiodic oscillation between healthy forest and urchin barrens, termed here “forest-barrens oscillation”, the latter condition being one in which algal densities are drastically reduced (Figure 3B). The barrens state is, however, unstable in our model because of continuous temperature-dependent recruitment of kelp and understory algae from external sources. Therefore, healthy forest may recover after transition to barrens, although the timescale varied among simulations.

FIGURE 3

Figure 3. Forest and forest-barrens states. Plot lines show species densities relative to density at the end of the simulation burn-in. Days 5,000–9,500 post burn-in are plotted. (A) Forest state, showing population trajectories of giant kelp (green), understory algae (blue), and purple sea urchin (purple). (B) Trajectories of the sunflower star (orange) and Steller's sea cow (brown) in forested state. (C,D) Species trajectories in the quasicyclic state. Color indications as in (A,B).

Forest, low diversity forest, and forest-barrens oscillation occurred at statistically indistinguishable frequencies (chi-squared test, χ² = 8.9215, p = 0.1780) across all model types, yielding frequencies of 66.5% forest, 25.5% forest-barrens oscillation, and 7.7% low diversity forest. Given that systems exhibiting quasiperiodic oscillation also exhibit intervals of forest, and that intervals spent in the forest condition during oscillations are at least as lengthy as barrens intervals, then under unperturbed conditions, areas with the potential to host M. pyrifera should be forested in at least 74.4% of occurrences.

Sea urchin ecological efficiency differed significantly among the states (Figure 4) and predicts when an unperturbed simulation will produce either the forest or forest-barrens state (Supplementary Table 2). The parameter varies stochastically among simulations with the addition of a small error (Supplementary Section 3). This is reflected by the significantly different distributions of the parameter among the states (ANOVA, F = 55.32, p < 0.0001, Scheffe's multicomparison test; Figure 4A). The dominance of urchins in the system thus increases with their ecological efficiency.

FIGURE 4

Figure 4. (A) Distributions of sea urchin ecological efficiency per the three states of unperturbed systems. Unstable forest-barrens systems have the highest values. (B) Logistic regression of the high diversity forest and forest-barrens simulations against sea urchin ecological efficiency.

We evaluated the dependence of transitions between high diversity forest and quasiperiodicity on ecological efficiency with a logistic regression (χ² = 47.65, R² = 0.441, p < 0.00001; Figure 4B). Quasiperiodicity emerges at an efficiency of approximately 0.102, with an increasing probability of emergence as efficiency increases. The healthy forest state is still possible at those values of efficiency though, and the system is therefore capable of existing in one of the two alternative states at higher efficiency values.

3.1.1. Forest types

A comparison of population sizes of species common to all the models (M. pyrifer, C. corymbiferus, S. purpuratus, and P. helianthoides) across unperturbed models, during the final year of each simulation shows that the forest state itself differs significantly among the models (MANOVA, Wilks' λ = 0.4043, F = 9286.86, p < 0.00001). The differences are driven by an inverse relationship between giant kelp and understory algal abundances (Figure 5 and Supplementary Figure 2). Models or systems without sea otter predation have greater abundances of M. pyrifera relative to C. corymbiferus. The Historical model with sea otters has the lowest relative abundance of giant kelp, whereas the Historical model without sea otters has the highest relative abundance. Thus the forest state across all models is itself heterogeneous.

FIGURE 5

Figure 5. Forest types produced by each model type. Box plots are of average species population density during the final year of 100 unperturbed simulations. Left boxes (green)—giant kelp; right boxes (blue)—understory algae.

3.2. Response to perturbation

3.2.1. State transitions

Responses to perturbation were classified as either being in a forest or forest-barrens state at the end of the simulation (Figures 6A–D). Frequencies (Supplementary Table 3) were compared among model-perturbation combinations. Model types did not differ in the frequency with which they produce the forest state when unperturbed (previous section), therefore any differences when perturbed are caused by perturbations only, and do not reflect model variance.

FIGURE 6

Figure 6. Resistance, transition and resilience in response to a combined perturbation by warming and sea star wasting disease. (A) Resistant system. The sea urchin population (purple) responds positively to both warming and release from sunflower star predation, but giant kelp (green) and understory algal (blue) populations do not decline. (B) Sunflower star (orange) and Steller's sea cow (brown) populations corresponding to those in (A). Note the slight increase of the sunflower star population in response to warming, before the steep decline caused by the wasting disease. The population recovers quickly after cessation of the epidemic. The cyclicity of the sea urchin and sunflower star populations is annual, being driven by seasonal temperature variation. (C,D) A system that undergoes transition to the forest-barrens state in response to the perturbation. All populations become quasiperiodic after the transition, but the periodicity does not correspond to seasonal temperature variation, occurring at a lower frequency. (E,F) A resilient system, which transitions to quasiperiodicity, but recovers after only two cycles, a duration of approximately 4 years.

Almost all models exhibited similar frequencies of transition to the forest-barrens state, with the exception of the Historical model when subjected to SSWD+MHW. This model produced a significantly greater transition frequency when compared to other models (χ² comparisons, α = 0.05). The pooled frequency of all other model-perturbation combinations is 3.29%, whereas it differs significantly at 32.23% when SSWD+MHW is applied to the Historical model (χ² = 36.718, p < 0.00001).

3.2.2. Resilience and resistance

Resilient or resistant simulations began in the forest state prior to perturbation and ended in that state. Simulations were classified as resistant if barrens never appeared (Figures 6A,B), or resilient if barrens occurred but the system subsequently returned to stable forest (Figures 6E,F). Resilient systems always underwent multiple state switches between forest and barrens before returning permanently to the forest state.

The relative frequencies of resistant and resilient simulations varied significantly among model types and perturbations (Supplementary Table 4), and models fall into several groups (Figure 7). The group with the greatest number of resistant simulations was that perturbed by warming only, regardless of model type. In contrast, the group with the fewest resistant simulations comprised both the Historical and Modern models without sea otters, when subjected to SSWD+MHW. Intermediate were the Historical and Modern models with sea otters, when subjected to either SSWD, or SSWD+MHW. There is no significant difference between the Modern system and the Historical one when the latter is subjected to SSWD+MHW (compared to Modern SSWD and SSWD+MHW; χ² = 3.24 and 3.34; p=0.072 and 0.067, respectively). The Historical model subjected to disease only, however, is significantly more resistant than the Modern models (compared to Modern SSWD and SSWD+MHW; χ² = 12.66 and 12.78; p=0.0004 and 0.0003, respectively).

FIGURE 7

Figure 7. Relative proportions of resistant and resilient simulations per model and perturbation type. Circles—Historical models, squares—Modern models, open symbols—systems without sea otter predation, closed symbols—systems with sea otter predation, red—warming perturbation, blue—sea star wasting disease perturbation, purple—both warming and disease. For example, the closed purple square are Modern systems, with sea otter predation, subjected to both warming and disease. Systems without sea otter predation, and subjected to disease or warming plus disease are the least resistant, whereas all systems subjected to warming only are the most resistant. Systems with sea otter predation, and subjected to either disease only, or warming plus disease, are intermediate, with the Historical system subjected to disease only being the most resistant.

3.2.3. Onset and recovery from transitions

Resilient simulations were characterized by comparing the times of onset of the first transition to barrens, and the time of permanent recovery to the forest state. Models perturbed by MHW were excluded from the analysis. There were no significant differences of either onset or recovery times among those models perturbed by SSWD (ANOVA, n = 144; onset, F = 1.06, p = 0.368; recovery, F = 0.68, p = 0.567), or SSWD+MHW (ANOVA, n = 131; onset, F = 0.70, p = 0.5; recovery, F = 0.88, p = 0.418), nor did models differ in recovery time if the perturbations were pooled (ANOVA, n = 293; recovery, F = 0.61, p = 0.747). Onset time however, differed among the pooled Historical models (ANOVA, n = 293; onset, F = 2.06, p = 0.048), with earlier onset times when sea otters are absent (Supplementary Figure 2). This is also true for the Modern model when the perturbation is SSWD+MHW (Table 1). Sea otters thus delay the onset of transitions to urchin barrens when Steller's sea cow is present, but fail to play a similar role in the Modern system when perturbed by SSWD only.

TABLE 1

Table 1. Mean onset times of transition to urchin barrens, dependent on model, perturbation, and the absence or presence of sea otters.

3.3. Predicting transition

Pre-perturbation conditions might determine whether simulations are resistant or resilient, potentially serving as predictive indicators of how communities would respond to perturbation. We tested this with a principal components analysis (PCA) of population sizes of the four species common to all models: M. pyrifera, C. corymbiferus, S. purpuratus and P. helianthoides. We analyzed population sizes during the final year before the onset of perturbation (SSWD, or SSWD+MHW). The PCA differentiates between resistant and resilient simulations (Figure 8A). Loadings of the variables indicate that the relative abundances of M. pyrifera and C. corymbiferus are the most important determinants of whether a response will be resistant or resilient (Supplementary Tables 5, 6).

FIGURE 8

Figure 8. Principal components analysis of resistant Modern simulations. (A) Principal components scores of community states during the final year before onset of perturbations (SSWD, and SSW disease plus MHW warming). Blue points are resistant simulations, and orange points are resilient but not resistant. Vector arrows show principal component alignment of variable (population densities), and arrow lengths correspond to relative eigenvector loadings. (B) Box plots of distributions of understory algae:kelp population density ratios under model-perturbation combinations, and correspond to community states in (A). Color code as in (A).

We compared ratios of understory algae to kelp population sizes during the year before perturbation, demonstrating that those ratios differ significantly between resistant and resilient responses among all the Modern model types (Figure 8B). In each model type, that is with or without sea otters, and subjected to either type of perturbation, systems that were subsequently resistant to any transient transitions to quasiperiodicity had significantly greater ratios of understory algae to kelp. Furthermore, understory algae have larger population sizes in the presence of sea otters.

3.4. Sensitivity analyses

Variance of the intrinsic rates of increase for M. pyrifera and C. corymbiferous showed that forest states in the sensitivity simulations were the same as those generated in the standard unperturbed model. Sensitivity to consumer mortality rates indicate that while model types vary in sensitivity to this parameter, none showed sensitivity at or below 1% variation, and thus all values in the standard model are consistent with results reported here (see Supplementary material).

4. Discussion

The model generated several broad conclusions, summarized here and discussed in more detail below. First, the M. pyrifera forest state is metastable, capable of seemingly spontaneous transitions to urchin barrens, or dramatic declines of kelp density, although such transitions are transient. Second, there are multiple types of forest, defined by the relative abundances of M. pyrifera and C. corymbiferus. As expected, kelp dominates the understory algae, although the extent of dominance varies according to the presence or absence of both sea otters and Steller's sea cow. Exceptionally, historical systems with both mammals present would have had nearly equal dominance of kelp and understory algae. Third, both SSWD, and SSWD+MWH could drive transitions from forest to a forest-barrens state where forest alternates with urchin barrens. The probability of transition depends critically on the composition of the system, as well as the nature of the perturbation. In general, M. pyrifera forests are highly resistant against transition, and are otherwise resilient, with resilience being greater in the presence of Steller's sea cow.

4.1. Alternative states without perturbation

Unperturbed simulations of all model types (Historical, Historical without otters, Modern, and Modern without otters) revealed three possible alternative states of M. pyrifera giant kelp communities (Figure 3), termed here “low diversity forest”, “high diversity forest” (or simply forest), and “forest-barrens”, the latter of which transitions quasiperiodically between a forested community and urchin barrens. Why are the models capable of generating different states if unperturbed? Sea urchin ecological efficiency emerged as the single explanatory parameter (Supplementary Table 2). This parameter varies to reflect uncertainty and true variation of its value in vivo. The generation of three discrete states represents a probabilistic system bifurcation, in which state depends on a noisy parameter (Mirzakhalili and Epureanu, 2019). Higher values of urchin efficiency increase the likelihood of quasiperiodic variation and hence urchin barrens. The ecological implication is that any factors leading to an increase of that efficiency, or to its greater influence on population growth, are also likely to give rise to quasiperiodic transitions between kelp forest and urchin barrens. Sea urchin ecological efficiency is temperature dependent (see Supplementary Equation 1). Therefore, ocean warming or cooling toward optimum temperature increases the probability of quasiperiodicity.

4.2. Multiple types of forest

The models also reveal several high diversity forest types, differentiated by the relative densities of giant kelp and understory algae. There is an inverse relationship between kelp and understory algal abundances in three of the model types, with the “Historical with sea otters” model being an exception relative to the others (Figure 5). The models therefore predict that kelp forests today should vary in the relative abundances of giant kelp and understory algae that are favored by urchins, and that the variance is dependent on the presence or absence of predatory otters. We are not aware of any prior theoretical or empirical studies that have suggested this variance, but it is amenable to empirical verification. The role of the sea otter is not surprising, but conventional thought is that their positive impact is the result of controlling sea urchins (Estes and Palmisano, 1974; Estes and Duggins, 1995). Here, we suggest an additional route whereby decreased grazing on understory algae by urchins in the presence of sea otters allows the algae to attain greater densities, being limited then solely by shading and competition with kelp. Understory algae thereby become a more important alternative food resource for urchins. Thus, the sea otter as a keystone predator may be important for persistence of the forest state and increased primary producer diversity via both direct and indirect interactions.

4.2.1. A no-analog community

Our results further predict a novel type of forest that would have existed when Steller's sea cow was extant. The sirenian would have had a positive effect on understory algal abundance, further amplified by the sea otter, resulting in kelp forests where the abundance of understory algae could match or even exceed that of the kelp. Today there are no megaherbivores that feed on kelp forest canopies, and herbivory by invertebrates and herbivorous fish are unlikely to open the canopy to the extent that Steller's sea cow was able to. The decline and extinction of Hydrodamalis species would therefore have marked the demise of a type of kelp forest for which there is no analog today. Results suggest that this would have had significant effects on both forest resilience and the probability of transitions to barrens when perturbed by MHW and SSWD.

4.3. Resilience and transition

Several parameters describing the ecologies of kelp, understory algae, and the invertebrate echinoderms are temperature sensitive in the model and decline as temperature increases. These include the intrinsic rates of population increase of the producers, and the predator-prey interaction coefficients and ecological efficiencies of consumers. For example, increases of water temperature decreases P. helianthoides growth rates (Fernández et al., 2021). Furthermore, SSWD reduces the interaction strength between sea urchins and sunflower sea stars, and rapid declines of P. helianthoides have been documented to lead to increases of the fraction of urchin populations that are exposed and foraging, and declines of the urchins' nutritional state (Smith et al., 2021). In the Modern system, urchin populations would then be controlled by sea otters only, and this is apparently sufficient to maintain the system in a forested state. In the Historical system, however, Steller's sea cow continues to exert a negative impact on the kelp both through consumption of fronds, and subsequent release of understory algae from shading. The combined effects of urchin release from sea star predation, depression of kelp growth rates, and sustained sea cow grazing, is apparently sufficient to tip the Historical system into the forest-barrens state.

The Historical system when subjected to SSWD only does not have a high probability of transition because the kelp rate of intrinsic increase is not depressed by warming. Paradoxically, when the system lacks sea otter predation, it is also highly resistant to SSWD+MHW, which would appear to be inconsistent with the above explanation because urchins are now under no predation pressure. We propose a type of intermediate disturbance hypothesis (IDH) (Connell, 1978) to explain these results. The classic IDH seeks to explain the maximization of biodiversity at intermediate levels of disturbance by a minimization of competition. Here, we propose that the probability of transition to quasiperiodicity of the Historical system is maximized at an intermediate level of disturbance. The lowest level of disturbance occurs when sea otters are present and the system is perturbed by disease only; resistance is maximized because kelp and understory algae are not affected by the disease (Figure 9). The intermediate disturbance and greatest probability of transition arises when systems are subjected to SSWD+MHW, and urchins are released from sea star predation but not predation by sea otters. Maximum disturbance occurs under SSWD+MHW if sea otters are absent. In this case, we propose that the urchin population is so tightly coupled to the dynamics of its prey (kelp and understory algae) prior to any disturbance, that warming weakens this producer-herbivore coupling and the impact of the urchins on their prey sufficiently to prevent the system from transitioning easily into quasiperiodicity.

FIGURE 9

Figure 9. (A) Schematic of intermediate disturbance hypothesis of Historical model-perturbation combination responses to perturbation. Hd, Historical community with sea otters, perturbed by SSWD; Hwd, same model perturbed by SSWD and MHW warming; Hwd, no otters, Historical model without sea otters, perturbed by SSWD and MHW warming. The y-axis expresses the hypothesis that transitions to a permanent forest-barrens state are most probable under the Hwd circumstances. Network diagrams show each corresponding system. Arrow thickness indicates relative strength of interaction. (B) Hd; (C) Hwd; (D) Hwd, no otters.

4.4. Resistance and resilience

Resistant systems are of great interest because although sea urchin populations may increase during outbreaks of SSWD or MHW (Figure 6A), resistant communities will not transition to barrens. Three conclusions can be drawn from the frequency of resistant responses. First, results fall neatly into three groups: Responses to MHW were overwhelmingly resistant, whereas those perturbed by SSWD, or SSWD+MHW, were more likely to be resistant if sea otters were present (Figure 7). This suggests that warming during 2013-2015 by itself would not have caused widespread transitions of forests to urchin barrens, but that it did act synergistically with SSWD, perhaps exacerbating the effect of the disease (Harvell et al., 2019; McPherson et al., 2021) or accelerating spread of the disease. Second, both Historical and Modern forests subjected to SSWD or SSWD+MHW, are more likely to be resistant when sea otters are present (Figure 7), confirming the keystone role of this species in maintaining giant kelp communities in the forest state. Third, the presence of sea cows confers additional resistance to the system, which we refer to as the “sea cow effect” (Figure 10).

FIGURE 10

Figure 10. Relative frequencies of resistant (maximally resilient) simulations under all perturbation-model combinations. Bars in each perturbation group are in the following order: 1 (blue)—Historical with sea otters; 2 (light blue)—Historical without sea otters; 3 (brown)—Modern with sea otters; 4 (light brown)—Modern without sea otters. Red bars represent the greater resistance of Historical vs. Modern systems when sea cows are present, termed the sea cow effect.

4.4.1. Onset of transitions

The best of all possible worlds for maintaining high density, ecologically resistant or resilient M. pyrifera forests when perturbed by SSWD, or SSWD+MHW, occurs in the presence of both Steller's sea cow and sea otters. In the present day, without the sea cow, the outcome depends on the perturbation, with the system faring best with sea otters present, and if perturbed by SSWD+MHW. The disease by itself affects the system similarly to SSWD+MHW in the absence of the otter. Although this is counterintuitive, with a combined perturbation sometimes having a lesser effect than a single one, and appearing to run counter to suggestions that warming alone can be a major driver of the transition to urchin barrens, there is a synergistic explanation (McPherson et al., 2021).

Onset time of quasiperiodicity is influenced most by the removal of sunflower sea stars (SSWD). Urchins thrive under that condition, even in the presence of sea otters. This is not the case when the Modern system is subjected to SSWD+MHW, because urchins too are affected negatively by warming and thus do not benefit from predation release to the same extent as when warming is absent. This distinction does not unfold in the Historic systems because Steller's sea cow moderates the system regardless of the presence or absence of sea otters, or the type of perturbation. Thus, the delay of a transition to urchin barrens when a Modern forest is disturbed is controlled by interactions of the sea urchins with their predators, but the delay of a transition would be enhanced in the presence of the extinct megaherbivore. Sea otters also had a delaying effect on the onset of transitions to quasiperiodicity in the Historic system. That effect is absent from the Modern system, highlighting an indirect synergistic interaction between the two mammals that was lost with the extinction of H. gigas.

4.4.2. Recovery

Recovery from the final transient barrens to a permanent forest state is invariant among the models, nor does it differ according to whether the perturbation is SSWD only, or SSWD+MHW. The mean time to recovery is approximately 9.6 years after the onset of SSWD. The disease itself persists for 10 years in the model, and P. helianthoides populations recover rapidly after that. Kelp forests are beginning to recover in 2022 at various locations in California, 9 years after the first appearance of SSWD (Bland, 2021). Many of these areas are the sites of conservation and restoration efforts that consist of large scale removal of urchins, and in some locations of northern California the re-planting of bull kelp, N. luetkeana (Williams et al., 2021; Eger et al., 2022). The congruent timing of these efforts with the average recovery time predicted by our modeling might be coincidental, but urchin removal does mimic the recovery of sunflower sea star populations (Williams et al., 2021), and the model supports the timing as being optimal for forest recovery. Filbee-Dexter and Scheibling (2014) have suggested that the barrens state may be unstable because of predation, a contention supported by the positive impact of recovering P. helianthoides populations in our model. Unfortunately, evidence for initial sea star recovery from SSWD is currently equivocal at best (Hamilton et al., 2021).

4.5. Restoring functionality and regenerating forests

Predictions of Steller's sea cow's impact on kelp forest dynamics depend on the type of perturbation, as well as organismal traits. Models predict that restoration of sea cow functionality in the form of reducing canopy cover could be key to both the resistance and resilience of M. pyrifera forests when sea otters are absent. Could a framework of enhanced resilience be implemented by operationalizing model results? Artificial restoration would depend on extensive experimental testing prior to any such actions. Given the complicated nature of the system's dynamics, the circumstances under which frond removal is done would depend on the perturbation, as well as the presence/absence of sea otters. The dynamics of unperturbed central Californian kelp forests are not steady state, as shown by multi-decadal records of canopy cover (Bell et al., 2015). In fact it could be argued that quasiperiodicity on a timescale longer than predicted by our model might be common. Both our model and canopy cover records from northern California (McPherson et al., 2021) suggest that MHW alone may be insufficient to trigger the transitions observed in the past 10 years. Instead, a decline in predation pressure and/or a synergy with warming are more likely culprits.

Given uncertainty about H. gigas consumption rates, the intensity of artificial browsing by canopy reduction would have to be established experimentally. Experiments would also be necessary to assess the roles of other species; for example, the impacts of urchin predators such as the sheephead Semicossyphus pulcher (Eisaguirre et al., 2020). Finally, experiments could facilitate the re-direction of kelp productivity back into the forest community. Today, significant amounts of kelp productivity are exported offshore in the form of kelp detritus, where it is believed to comprise a major nutrient supply to the offshore water column (Smale et al., 2021). That export would have been less in the presence of the sea cow, and the egestion and excretion of the megaherbivore would have created a significant nutrient feedback directly into the forest community that is no longer present. Harvesting of the kelp canopy with local retention of degraded detritus could therefore recreate some ecological functionality of the sea cow, increasing forest resilience under certain circumstances while providing an additional nutrient source to the community.

4.6. Proposing a new approach

What is the current status (e.g., diversity, resilience, integrity) of an ecosystem? What did this ecosystem look like in the past? Given current trends of anthropogenic change, what kinds of interventions, in what order and over what time period, are most likely to generate positive results? Answering these questions requires a PPF lens that extends the standard time frame of conservation, and utilizes multiple tools and types of data. Conservation paleobiology is the temporal framework for exploring paleoecological snapshots (Dietl and Flessa, 2011) that extend at least as far back as the Pleistocene or early Holocene. Mathematical modeling, which is at the core of this study, has potential to interweave diverse data streams (Scheffer, 2020) and, in doing so, provide answers to the above questions. The methodology advocated here can be summarized as follows. First, establish basic ecological models of a present-day ecosystem, including both environmental conditions and biological components. Second, integrate historical and/or paleontological data to explore putative past states of that ecosystem, emphasizing key elements that no longer exist. Third, integrate the first and second steps to create mathematical models of possible future ecosystem states. Finally, generate predictions about the likely impacts of specific interventions in order to inform experimentation, policies, and action.

5. Summary

As anthropogenically-driven environmental changes multiply, intensify and become more permanent features of the natural world, we will need an increasingly diverse and powerful set of tools to address the conservation and regeneration of natural systems. Similarly, as species decline, and extirpations and extinctions increase, historical and paleontological studies must be brought to bear to understand the altered states of modern ecosystems. We believe that the PPF modeling approach advocated in this paper has considerable potential to address these concerns. Here, we combined paleontological, historical and ecological data into a mathematical model of how an extinct marine megaherbivore might have affected the dynamics and resilience of an important coastal marine ecosystem. Those effects were generally non-trivial and often counterintuitive. On one hand, the sea cow may have increased the probability that giant kelp forests would switch to an alternative state exhibiting periodic urchin barrens when subjected to marine heat waves and sea star wasting disease. On the other hand, forests that did not transition permanently to an alternative state would have been significantly more resilient compared to the present. When both sea cows and sea otters were present in historical forests, they would have interacted indirectly to increase resilience, a function that the sea otters today apparently do not fulfill.

There is a growing appreciation of the roles of extant and extinct marine megafauna as ecosystem engineers (Pyenson and Vermeij, 2016). Direct impacts of large vertebrates include physical habitat restructuring, consumption of large quantities of biomass at or near the base of food chains, facilitation of processes vital to the survival of other species, and flux rates and transport of nutrients. To this may be added the browsing of giant kelp canopies by Steller's sea cow. Important functions of megafauna may extend well beyond such direct impacts, however, and include the indirect interactions and impacts of these powerful consumers that affect community resilience and alternative state dynamics. Efforts to conserve and regenerate natural systems can, in our opinion, benefit from consideration of novel, mathematically rigorous interventions like the one proposed here.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://osf.io/zvnw9/.

Author contributions

RB, PR, and SS conceived of the study. PR developed the mathematical model and wrote simulation code. RB and PR conducted statistical analyses. All authors contributed to the final manuscript and approved the final submission.

Acknowledgments

The authors thank the helpful comments and suggestions of three reviewers. The authors acknowledge that the California Academy of Sciences occupies land that is the ancestral homeland of the Ramaytush Ohlone, the vast majority of whom were forcibly removed from their homelands. We acknowledge this painful history and honor their legacy and continuing contributions to better understanding the natural history of the United States.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.983558/full#supplementary-material

References

Arditi, R., and Ginzburg, L. R. (2012). How Species Interact: Altering the Standard View on Trophic Ecology. New York, NY: Oxford University Press. doi: 10.1093/acprof:osobl/9780199913831.001.0001

CrossRef Full Text | Google Scholar

Bell, T. W., Cavanaugh, K. C., Reed, D. C., and Siegel, D. A. (2015). Geographical variability in the controls of giant kelp biomass dynamics. J. Biogeogr. 42, 2010–2021. doi: 10.1111/jbi.12550

CrossRef Full Text | Google Scholar

Bland, A. (2021). Kelp Forests Surge Back on Parts of the North Coast, With a Lesson About Environmental Stability. Technical report, Bay Nature. Available online at: https://baynature.org/2021/09/13/kelp-forests-surge-back-on-the-north-coast-with-a-lesson-about-stable-environments/

Google Scholar

Bullen, C. D., Campos, A. A., Gregr, E. J., McKechnie, I., and Chan, K. M. (2021). The ghost of a giant-Six hypotheses for how an extinct megaherbivore structured kelp forests across the North Pacific Rim. Glob. Ecol. Biogeogr. 30, 2101–2118. doi: 10.1111/geb.13370

CrossRef Full Text | Google Scholar

Burt, J. M., Tinker, M. T., Okamoto, D. K., Demes, K. W., Holmes, K., and Salomon, A. K. (2018). Sudden collapse of a mesopredator reveals its complementary role in mediating rocky reef regime shifts. Proc. R. Soc. B Biol. Sci. 285, 20180553. doi: 10.1098/rspb.2018.0553

PubMed Abstract | CrossRef Full Text | Google Scholar

Castorani, M. C., Reed, D. C., and Miller, R. J. (2018). Loss of foundation species: disturbance frequency outweighs severity in structuring kelp forest communities. Ecology 99, 2442–2454. doi: 10.1002/ecy.2485

PubMed Abstract | CrossRef Full Text | Google Scholar

Connell, J. H. (1978). Diversity in tropical rain forests and coral reefs: high diversity of trees and corals is maintained only in a nonequilibrium state. Science 199, 1302–1310. doi: 10.1126/science.199.4335.1302

PubMed Abstract | CrossRef Full Text | Google Scholar

Dayton, P. K., Tegner, M. J., Edwards, P. B., and Riser, K. L. (1998). Sliding baselines, ghosts, and reduced expectations in kelp forest communities. Ecol. Appl. 8, 309–322. doi: 10.1890/1051-0761(1998)008[0309:SBGARE]2.0.CO;2

CrossRef Full Text | Google Scholar

Detmer, A. R., Miller, R. J., Reed, D. C., Bell, T. W., Stier, A. C., and Moeller, H. V. (2021). Variation in disturbance to a foundation species structures the dynamics of a benthic reef community. Ecology 102, e03304. doi: 10.1002/ecy.3304

PubMed Abstract | CrossRef Full Text | Google Scholar

Dietl, G. P., and Flessa, K. W. (2011). Conservation paleobiology: putting the dead to work. Trends Ecol. Evol. 26, 30–37. doi: 10.1016/j.tree.2010.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Domning, D. P. (1976). An ecological model for late tertiary sirenian evolution in the North Pacific Ocean. Syst. Zool. 25, 352–362. doi: 10.2307/2412510

CrossRef Full Text | Google Scholar

Eger, A. M., Marzinelli, E. M., Christie, H., Fagerli, C. W., Fujita, D., Gonzalez, A. P., et al. (2022). Global kelp forest restoration: past lessons, present status, and future directions. Biol. Rev. Camb. Philos. Soc. 97, 1449–1475. doi: 10.1111/brv.12850

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisaguirre, J. H., Eisaguirre, J. M., Davis, K., Carlson, P. M., Gaines, S. D., and Caselle, J. E. (2020). Trophic redundancy and predator size class structure drive differences in kelp forest ecosystem dynamics. Ecology 101, e02993. doi: 10.1002/ecy.2993

PubMed Abstract | CrossRef Full Text | Google Scholar

Estes, J. A., Burdin, A., and Doak, D. F. (2016). Sea otters, kelp forests, and the extinction of Steller's sea cow. Proc. Natl. Acad. Sci. U.S.A. 113, 880–885. doi: 10.1073/pnas.1502552112

PubMed Abstract | CrossRef Full Text | Google Scholar

Estes, J. A., and Duggins, D. O. (1995). Sea otters and kelp forests in Alaska: generality and variation in a community ecological paradigm. Ecol. Monogr. 65, 75–100. doi: 10.2307/2937159

CrossRef Full Text | Google Scholar

Estes, J. A., and Palmisano, J. F. (1974). Sea otters: their role in structuring nearshore communities. Science 185, 1058–1060. doi: 10.1126/science.185.4156.1058

PubMed Abstract | CrossRef Full Text | Google Scholar

Feehan, C. J., Grauman-Boss, B. C., Strathmann, R. R., Dethier, M. N., and Duggins, D. O. (2018). Kelp detritus provides high-quality food for sea urchin larvae. Limnol. Oceanogr. 63, S299–S306. doi: 10.1002/lno.10740

CrossRef Full Text | Google Scholar

Fernández, P. A., Navarro, J. M., Camus, C., Torres, R., and Buschmann, A. H. (2021). Effect of environmental history on the habitat-forming kelp Macrocystis pyrifera responses to ocean acidification and warming: a physiological and molecular approach. Sci. Rep. 11, 1–15. doi: 10.1038/s41598-021-82094-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Filbee-Dexter, K., and Scheibling, R. E. (2014). Sea urchin barrens as alternative stable states of collapsed kelp ecosystems. Mar. Ecol. Prog. Ser. 495, 1–25. doi: 10.3354/meps10573

CrossRef Full Text | Google Scholar

Foster, M. C., Byrnes, J. E., and Reed, D. C. (2015). Effects of five southern California macroalgal diets on consumption, growth, and gonad weight, in the purple sea urchin Strongylocentrotus purpuratus. PeerJ 3, e719. doi: 10.7717/peerj.719

PubMed Abstract | CrossRef Full Text | Google Scholar

Foster, M. S., Reed, D. C., Carr, M. H., Dayton, P. K., Malone, D. P., Pearse, J. S., et al. (2013). “Kelp forests in California,” in Research and Discoveries: The Revolution of Science through Scuba, Vol. 39 of Smithsonian Contributions to the Marine Sciences. Smithsonian Institution Press.

Gill, J. L. (2014). Ecological impacts of the late quaternary megaherbivore extinctions. N. Phytol. 201, 1163–1169. doi: 10.1111/nph.12576

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, M. H. (2004). Effects of local deforestation on the diversity and structure of southern California giant kelp forest food webs. Ecosystems 7, 341–357. doi: 10.1007/s10021-003-0245-6

CrossRef Full Text | Google Scholar

Graham, M. H., Harrold, C., Lisin, S., Light, K., Watanabe, J. M., and Foster, M. S. (1997). Population dynamics of giant kelp Macrocystis pyrifera along a wave exposure gradient. Mar. Ecol. Prog. Ser. 148, 269–279. doi: 10.3354/meps148269

CrossRef Full Text | Google Scholar

Graham, M. H., Kinlan, B. P., Druehl, L. D., Garske, L. E., and Banks, S. (2007). Deep-water kelp refugia as potential hotspots of tropical marine diversity and productivity. Proc. Natl. Acad. Sci. U. S. A. 104, 16576–16580. doi: 10.1073/pnas.0704778104

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimm, V., and Wissel, C. (1997). Babel, or the ecological stability discussions: an inventory and analysis of terminology and a guide for avoiding confusion. Oecologia 109, 323–334. doi: 10.1007/s004420050090

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamilton, S. L., Saccomanno, V., Heady, W., Gehman, A., Lonhart, S., Beas-Luna, R., et al. (2021). Disease-driven mass mortality event leads to widespread extirpation and variable recovery potential of a marine predator across the eastern Pacific. Proc. R. Soc. B 288, 20211195. doi: 10.1098/rspb.2021.1195

PubMed Abstract | CrossRef Full Text | Google Scholar

Harvell, C. D., Montecino-Latorre, D., Caldwell, J. M., Burt, J. M., Bosley, K., Keller, A., et al. (2019). Disease epidemic and a marine heat wave are associated with the continental-scale collapse of a pivotal predator (Pycnopodia helianthoides). Sci. Adv. 5, eaau7042. doi: 10.1126/sciadv.aau7042

PubMed Abstract | CrossRef Full Text | Google Scholar

Hyvarinen, O., Te Beest, M., le Roux, E., Kerley, G., de Groot, E., Vinita, R., et al. (2021). Megaherbivore impacts on ecosystem and earth system functioning: the current state of the science. Ecography 44, 1579–1594. doi: 10.1111/ecog.05703

CrossRef Full Text | Google Scholar

Karatayev, V. A., Baskett, M. L., Kushner, D. J., Shears, N. T., Caselle, J. E., and Boettiger, C. (2021). Grazer behaviour can regulate large-scale patterning of community states. Ecol. Lett. 24, 1917–1929. doi: 10.1111/ele.13828

PubMed Abstract | CrossRef Full Text | Google Scholar

Kriegisch, N., Reeves, S., Flukes, E., Johnson, C., and Ling, S. (2019). Drift-kelp suppresses foraging movement of overgrazing sea urchins. Oecologia 190, 665–677. doi: 10.1007/s00442-019-04445-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Mann, K. (1977). Destruction of kelp-beds by sea-urchins: a cyclical phenomenon or irreversible degradation? Helgoländer Wissenschaftliche Meeresuntersuchungen 30, 455–467. doi: 10.1007/BF02207854

CrossRef Full Text | Google Scholar

Marsh, H., O'Shea, T. J., and Reynolds, J. E. III (2012). Ecology and Conservation of the Sirenia: Dugongs and Manatees. New York, NY: Cambridge University Press. doi: 10.1017/CBO9781139013277

CrossRef Full Text | Google Scholar

McPherson, M. L., Finger, D. J., Houskeeper, H. F., Bell, T. W., Carr, M. H., Rogers-Bennett, L., et al. (2021). Large-scale shift in the structure of a kelp forest ecosystem co-occurs with an epizootic and marine heatwave. Commun. Biol. 4, 1–9. doi: 10.1038/s42003-021-01827-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirzakhalili, E., and Epureanu, B. I. (2019). Probabilistic analysis of bifurcations in stochastic nonlinear dynamical systems. J. Comput. Nonlinear Dyn. 14, 081009. doi: 10.1115/1.4043669

PubMed Abstract | CrossRef Full Text | Google Scholar

Mogen, S. C., Lovenduski, N. S., Dallmann, A. R., Gregor, L., Sutton, A. J., Bograd, S. J., et al. (2022). Ocean biogeochemical signatures of the North Pacific Blob. Geophys. Res. Lett. 49, e2021GL096938. doi: 10.1029/2021GL096938

CrossRef Full Text | Google Scholar

Owen-Smith, R. N. (1988). Megaherbivores: The Influence of Very Large Body Size on Ecology. New York, NY: Cambridge University Press. doi: 10.1017/CBO9780511565441

CrossRef Full Text | Google Scholar

Pedersen, M. F., Filbee-Dexter, K., Norderhaug, K. M., Fredriksen, S., Frisk, N. L., Wernberg, T., et al. (2020). Detrital carbon production and export in high latitude kelp forests. Oecologia 192, 227–239. doi: 10.1007/s00442-019-04573-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Pimm, S., and Lawton, J. (1977). Number of trophic levels in ecological communities. Nature 268, 329–331. doi: 10.1038/268329a0

CrossRef Full Text | Google Scholar

Pyenson, N. D., and Vermeij, G. J. (2016). The rise of ocean giants: maximum body size in Cenozoic marine mammals as an indicator for productivity in the Pacific and Atlantic Oceans. Biol. Lett. 12, 20160186. doi: 10.1098/rsbl.2016.0186

PubMed Abstract | CrossRef Full Text | Google Scholar

Reed, D. C., and Foster, M. S. (1984). The effects of canopy shadings on algal recruitment and growth in a giant kelp forest. Ecology 65, 937–948. doi: 10.2307/1938066

CrossRef Full Text | Google Scholar

Reed, D. C., Rassweiler, A., Carr, M. H., Cavanaugh, K. C., Malone, D. P., and Siegel, D. A. (2011). Wave disturbance overwhelms top-down and bottom-up control of primary production in California kelp forests. Ecology 92, 2108–2116. doi: 10.1890/11-0377.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Rogers-Bennett, L., and Catton, C. (2019). Marine heat wave and multiple stressors tip bull kelp forest to sea urchin barrens. Sci. Rep. 9, 1–9. doi: 10.1038/s41598-019-51114-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheffer, M. (2020). Critical Transitions in Nature and Society, Vol. 16. Princeton, NJ: Princeton University Press. doi: 10.2307/j.ctv173f1g1

CrossRef Full Text | Google Scholar

Smale, D. A., Pessarrodona, A., King, N., and Moore, P. J. (2021). Examining the production, export, and immediate fate of kelp detritus on open-coast subtidal reefs in the Northeast Atlantic. Limnol. Oceanogr. doi: 10.1002/lno.11970. [Epub ahead of print].

Smith, J. G., and Tinker, M. T. (2022). Alternations in the foraging behaviour of a primary consumer drive patch transition dynamics in a temperate rocky reef ecosystem. Ecol. Lett. 25, 1827–1838. doi: 10.1111/ele.14064

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, J. G., Tomoleoni, J., Staedler, M., Lyon, S., Fujii, J., and Tinker, M. T. (2021). Behavioral responses across a mosaic of ecosystem states restructure a sea otter-urchin trophic cascade. Proc. Natl. Acad. Sci. U.S.A. 118, e2012493118. doi: 10.1073/pnas.2012493118

PubMed Abstract | CrossRef Full Text | Google Scholar

Spindel, N. B., Lee, L. C., and Okamoto, D. K. (2021). Zombies of the nearshore: metabolic depression in sea urchin barrens associated with food deprivation. Bull. Ecol. Soc. Am. 102, e01926. doi: 10.1002/bes2.1926

CrossRef Full Text | Google Scholar

Steller, G. W. (1751). De Bestiis Marinis, or, the Beasts of the Sea (1751). St. Petersburg, FL: Press of the Academy of Sciences.

Steneck, R. S., Graham, M. H., Bourque, B. J., Corbett, D., Erlandson, J. M., Estes, J. A., et al. (2002). Kelp forest ecosystems: biodiversity, stability, resilience and future. Environ. Conserv. 29, 436–459. doi: 10.1017/S0376892902000322

CrossRef Full Text | Google Scholar

Walker, B., Holling, C. S., Carpenter, S. R., and Kinzig, A. (2004). Resilience, adaptability and transformability in social-ecological systems. Ecol. Soc. 9, 5. doi: 10.5751/ES-00650-090205

CrossRef Full Text | Google Scholar

Williams, J. P., Claisse, J. T., Pondella, D. J. II., Williams, C. M., Robart, M. J., et al. (2021). Sea urchin mass mortality rapidly restores kelp forest communities. Mar. Ecol. Prog. Ser. 664, 117–131. doi: 10.3354/meps13680

CrossRef Full Text | Google Scholar

Wing, S. R., Leichter, J. J., and Denny, M. W. (1993). A dynamic model for wave-induced light fluctuations in a kelp forest. Limnol. Oceanogr. 38, 396–407. doi: 10.4319/lo.1993.38.2.0396

CrossRef Full Text | Google Scholar