# Understanding the depth limit of the seagrass Cymodocea nodosa as a critical transition: Field and modeling evidence

Elvira Mayol, Jordi Boada, Marta Pérez, Neus Sanmartí, Mario

Minguito-Frutos, Rohan Arthur, Teresa Alcoverro, David Alonso, Javier
Romero

► To cite this version:

Elvira Mayol, Jordi Boada, Marta Pérez, Neus Sanmartí, Mario Minguito-Frutos, et al. Understanding the depth limit of the seagrass Cymodocea nodosa as a critical transition: Field and modeling evidence. Marine Environmental Research, 2022, 182, pp.105765. 10.1016/j.marenvres.2022.105765. hal-03956334

HAL Id: hal-03956334

https://hal.science/hal-03956334

Submitted on 26 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Understanding the depth limit of the seagrass *Cymodocea nodosa* as a critical transition: Field and modeling evidence

Elvira Mayol<sup>a, *</sup>, Jordi Boada<sup>b</sup>, Marta Pérez<sup>c</sup>, Neus Sanmartí<sup>c</sup>, Mario Minguito-Frutos<sup>d</sup>, Rohan Arthur<sup>d, e</sup>, Teresa Alcoverro<sup>d, e</sup>, David Alonso<sup>d</sup>, Javier Romero<sup>c</sup>

## ARTICLE INFO

**Keywords:**

Seagrass

Light attenuation

Plant traits

Tipping points

Alternative stable states

Resilience

## ABSTRACT

Changes in light and sediment conditions can sometimes trigger abrupt regime shifts in seagrass meadows resulting in dramatic and unexpected die-offs of seagrass. Light attenuates rapidly with depth, and in seagrass systems with non-linear behaviours, can serve as a sharp boundary beyond which the meadow transitions to bare sand. Determining system behaviour is therefore essential to ensuring resilience is maintained and to prevent stubborn critical ecosystem transitions caused by declines in water quality. Here we combined field and modelling studies to explore the transition from meadow to bare sand in the seagrass *Cymodocea nodosa* at the limit of its depth distribution in a shallow, light-limited bay. We first describe the relationship between light availability and seagrass density along a depth gradient in an extensive unfragmented meadow (Alfacis bay, NE Spain). We then develop a simple mechanistic model to characterise system behaviour. In the field, we identified sharp decline in shoot density beyond a threshold of ~1.9 m depth, shifting from a vegetated state to bare sand. The dynamic population model we developed assumes light-dependent growth and an inverse density-dependent mortality due to facilitation between shoots (mortality rate decreases as shoot density increases). The model closely tracked our empirical observations, and both the model and the field data showed signs of bistability. This strongly suggests that the depth limit of *C. nodosa* is a critical transition driven by photosynthetic light requirements. While the mechanisms still need to be confirmed with experimental evidence, recognizing the non-linear behaviour of *C. nodosa* meadows is vital not only in improving our understanding of light effects on seagrass dynamics, but also in managing shallow-water meadows. Given the shallow threshold (<2m), light-limited systems may experience significant and recalcitrant meadow retractions with even small changes in sediment and light conditions. Understanding the processes underlying meadow resilience can inform the maintenance and restoration of meadows worldwide.

## Introduction

Ecosystems emerge as a product of species interactions – both with each other, as well as with their environment. These interactions often result in dynamic systems that show complex non-linear behaviours in response to changing conditions (Conversi et al., 2015). Many complex ecosystems are inherently surprising: relatively minor increments of a stressor can lead to abrupt and catastrophic shifts to another state (Scheffer et al. 2001). The shifted state (forward shift) involves a

complete reorganization of species-environment interactions, is usually simpler and less diverse, and is often functionally poorer. Non-linear system behaviours are usually characterised by stabilizing feedback mechanisms that confer a degree of resilience against collapse (Scheffer et al., 2001). When the simpler shifted state also reorganizes around a set of stabilizing mechanisms, the system displays hysteresis and remains resilient to recovery even when stressful conditions abate (Scheffer et al., 2001; Nyström et al., 2012). Such bistable systems can therefore exist in either ecosystem state at exactly the same level of

<footnote><sup>a</sup> Institut Mediterrani d'Estudis Avançats (IMEDEA-CSIC), Carrer Miquel Marqués 21, 07190, Esporles, Spain</footnote>

<footnote><sup>b</sup> Laboratoire d'Océanographie de Villefranche, Sorbonne Université, Villefranche-sur-Mer, France</footnote>

<footnote><sup>c</sup> Departament de Biologia Evolutiva, Ecologia i Ciències Ambientals, Facultat de Biologia, Universitat de Barcelona, Avinguda Diagonal 643, 08028, Barcelona, Spain</footnote>

<footnote><sup>d</sup> Centre d'Estudis Avançats de Blanes (CEAB-CSIC), Carrer d'Accés a la cala Sant Francesc 14, 17300, Blanes, Spain</footnote>

<footnote><sup>e</sup> Nature Conservation Foundation, 1311 Amritha, 12th Cross, Vijayanagara 1st Stage, Mysore, 570017, India</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: elvira.mayol@imedea.uib-csic.es (E. Mayol).</footnote>
stressor (Scheffer and Carpenter, 2003; Conversi et al., 2015), and shifts between alternative states occur through “critical transitions” (Scheffer et al., 2001).

Instances of critical transitions in both terrestrial or aquatic ecosystems are widespread, from savannas to lakes and coral reefs, and there is growing interest in understanding their underlying dynamics (Scheffer and Carpenter, 2003). They present a formidable challenge to our capacity to predict large changes in ecosystem state, and studies have focused on identifying an impending collapse and on understanding the feedback mechanisms that promote ecosystem resistance (Barbier et al., 2008; Duarte et al., 2009; Maxwell et al., 2017; Berdugo et al., 2021). However, we are still far from fully understanding the interacting drivers of abrupt changes in most ecosystems well enough to effectively predict (and, following predictions and taking adequate actions, mitigate and/or prevent) their collapse (Holling, 1973; Scheffer et al., 2001).

Critical transitions are relatively common in nearshore ecosystems, and reports of their collapse are growing increasingly frequent with the combined anthropogenic stressors of overfishing, pollution and global warming (Ling et al., 2015; Wernberg et al., 2016; Melis et al., 2019; Wang et al., 2020; Verdura et al., 2021). Examples of coastal ecosystem collapse include, for instance, the overharvesting of tropical herbivore fishes, that has resulted in shifts from coral to macroalgal-dominated bottom in Caribbean reefs (Kennedy et al., 2013), or the herbivorous sea urchin outbreaks in temperate areas caused by overfishing of their fish predators, and resulting in the expansion of barrens (Estes and Palmisano, 1974; Filbee-Dexter and Scheibling, 2014; Ling et al., 2015) or shifts in nutrient regimes (Tett et al., 2007). Critical transitions have also been observed in seagrass meadows, shifting from communities with lush vegetation to bare sand associated to an increase in sediment and nutrients supply (McGlathery et al., 2013).

Seagrass meadows are among the most structurally important and biodiverse habitats in coastal waters, and support coastal societies with a host of essential goods and services. Seagrass meadow collapses can represent significant losses to dependent human communities (McGlathery et al., 2007; Fourqurean et al., 2012; Christianen et al., 2013). Like other coastal systems, seagrass meadows are exposed to multiple anthropogenic stressors which act together to erode their resilience (Unsworth et al., 2015). These include changes in sediment biogeochemistry (Holmer and Hasler-Sheetal, 2014), increases in sediment resuspension (McGlathery et al., 2013), warming (Arias-Ortiz et al., 2018; Ontoria et al., 2019), eutrophication (Burkholder et al., 2007) and herbivory (Prado et al., 2007; Arnold et al., 2008) among others. While the decline of meadows to these stressors is relatively well documented, there is still need to understand if seagrass responses involve non-linear dynamics and critical transitions (but see Carr et al., 2012; McGlathery et al., 2013) and the mechanisms that produce these behaviors. The implications are non-trivial. Managing systems subject to critical transitions requires completely different strategies, since stabilizing feedbacks may mask an approaching threshold until it is too late for ameliorative action. Additionally, bistable systems are remarkably difficult to restore post-collapse, and it may require more than a mere reduction of stressor levels for them to shift back to a more socio-ecologically desirable state.

Large-scale meadow shifts to bare sand have been described linked to eutrophication (Connel et al., 2017) increased sediment resuspension (Carr et al., 2010; Connell et al., 2017) on light attenuation (Dennison, 1987). As stated earlier, the abruptness of these shifts suggests that the system has several feedback mechanisms that contribute, implicitly or explicitly, to the persistence of the vegetated state (Maxwell et al., 2017). For instance, the capacity of seagrasses to absorb nutrients and reduce their concentrations from the water column is a crucial mechanism preventing phytoplankton or epiphyte growth and promoting seagrass persistence (Nielsen et al., 2004). The presence of seagrass itself could indirectly facilitate light penetration supporting further seagrass growth (Carr et al., 2012). Moreover, seagrass canopies attenuate water flow and trap fine sediment particles, enhancing deposition and

reducing bottom susceptibility to resuspension, which feeds back to increased water transparency (de Boer, 2007). Seagrasses also alleviate the mechanical impact and erosion by hydrodynamic forces through the protection provided by the leaf canopy. In particular, local erosion should be greater at the edge of the meadow, compromising local plant growth (Peterson et al., 2004; Maxwell et al., 2017). Seagrasses are clonal species, and integration between ramets can be a powerful facilitative mechanism of meadow resilience (Kendrick et al., 2005). Because of positive ramet interactions, the rate of expansion of seagrass patches accelerates with size (Vidondo et al., 1997; Marbà and Duarte, 1998; Kendrick et al., 1999; Sintes et al., 2005). Feedbacks have been proposed between seagrass growth and the physical environment due to clonal integration (Duarte and Sand-Jensen 1990; Olesen and Sand-Jensen, 1994; Vidondo et al., 1997; Ramage and Schiel 1999). Taken together, clonal integration, physical protection and enhancement of water transparency, hint that seagrass meadows have several facilitative mechanisms that promote their resilience and can result in density-dependent survival in the meadow.

Seagrass abundance (shoot density and/or cover and/or biomass) usually decreases with depth; beyond a point the plant disappears, strongly linked to declining light availability (Duarte, 1991). This transition is often not a gradual thinning out, but, as seen in Fig. 1, can often be a sharp boundary (Dennison, 1987; Duarte, 1991; Mayot et al., 2005; Collier et al., 2007; Boudouresque et al., 2009; Enríquez et al., 2019), hinting at a non-linear relationship between seagrass abundance and light. If the depth limit of meadows represents a critical transition, it has far-reaching implications, not merely for advancing our understanding of light effects on seagrass dynamics, but also for the management of the ecosystem. The contraction of seagrass depth limits with declining environmental condition can lead to large losses in the areal extent of meadows (Dennison, 1987; Duarte, 1991; Nielsen et al., 2002; Duarte et al., 2007). In fact, the seagrass depth limit is commonly used in biomonitoring programmes related to eutrophication (López y Royo et al., 2010). For shallow-water seagrasses with narrow depth ranges, even a small reduction in depth limits can be potentially devastating. In this context, the regime shift framework represents not only an adequate tool to better understand the depth limit of seagrass distribution, but is also an opportunity for increasing our understanding of seagrass-light interactions and, to help identify mechanisms underlying critical transitions in seagrasses. The regime shift perspective can be combined with water clarity monitoring programs to detect conditions that endanger the ecosystem and, in the absence of any action, would lead to its collapse (Saunders et al., 2013). This understanding would contribute to reducing the unpredictability of massive seagrass losses as a consequence of relatively small changes in the light regime, caused either by natural or, increasingly, man-made disturbances, giving a chance to possible management actions aimed at preventing these losses.

Our working hypothesis here is that the transition from a seagrass meadow state to a bare sand state in the deep limit of meadow distribution can be best described by an abrupt change caused by a reduction in the available light beyond a given threshold. This abrupt change is based on the breakdown of resilience mechanisms linked to clonal integration. To test this hypothesis, we use complementary empirical and modelling approaches. First, we explore the type of transition between seagrass abundance (shoot density) and sand along the depth gradient in a dense and extensive meadow of *Cymodocea nodosa* in a shallow embayment (Fig. 1). Then, we develop a mechanistic model to approximately reproduce the observed behaviour of the system. The model explores the relation between plant abundance and light along a depth gradient. In the model we assume the existence of clonal integration between ramets, as well as other protective effects (see previous paragraphs and Nielsen et al., 2004; de Boer 2007; McGlathery et al., 2013) that result in a density-dependent avoidance of mortality, acting as a stabilizing feedback that contributes to the resilience of the vegetated state.
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Upper panel: Edge of Cymodocea nodosa meadow on its deep limit in Alfacs bay, Ebro Delta, showing a sharp transition between the meadow and the bare sand. Photo: Marta Pérez. Lower panel: Study area and sampling transects on a depth gradient (black bars) and detail of the sampled points (in yellow, n = 119).</figcaption></figure>

# Materials and methods

## 2.1. *Empirical approach*

### 2.1.1. *Study site and sampling methodology*

Our field observations were conducted in the Alfacs Bay (Ebro River Delta), a shallow embayment (up to 6 m deep) with an estimated area of 50 km², located in the Northwestern Mediterranean. The bay has extensive seagrass meadows of *C. nodosa* that grow from the surface to 2–3 m (Pérez and Camp, 1986). The embayment receives seasonal freshwater inputs as runoff from rice paddy fields, causing high nutrient and organic matter concentrations and suspended materials increasing water turbidity (Pérez et al., 1994).

Sampling took place in two consecutive years. To avoid seasonal variability, sampling was performed over a short period (3 days) in June 2018 and June 2019, when the seasonal peak in seagrass biomass and shoot density occurs (Mascaró et al., 2014). In order to empirically assess the relationship between plant abundance and available light, we

extensively sampled the meadows, to characterise the depth gradient (cross-shore, from 0.4 to 3.6 m) and potential horizontal variability (long-shore), with a total of 119 sampling points for the two years (Fig. 1). Sampling points were randomly selected at each depth range, and distributed along approximately 10 km of seagrass meadow (Fig. 1). At each sampling site, we collected all the plant material within a hand-held corer (15 cm internal diameter, pushed 30 cm into the sediment), according to Pérez et al. (2001). In 2018, each sample was thoroughly rinsed *in situ* with seawater to eliminate the sediment, and plant material sealed in plastic bags, immediately frozen and carried to the laboratory where it was conserved at –25 °C until processing, which consisted in counting the number of shoots. In 2019, in order to optimize sampling effort, the shoots were counted in the field, and no samples were transported to the laboratory. The number of shoots obtained at each sample was transformed to shoot density (shoots m⁻²) and then presented as a percentage relative to the maximum value of shoot density found (5375 shoots m⁻²).

At each sampling point and in both years, incident light as PAR
(photosynthetically active radiation) was measured using a LI-COR: LI-1000 Data Logger, immediately below the surface and at the top of the plant canopy (or 20 cm above the sediment in the absence of plants). Values were obtained integrating 1 min of readings, and percent light extinction (relative to subsurface irradiance) was computed for each point as a proxy for light deprivation stress. These light measurements are merely indicative, as they were taken close to the yearly maximum light availability, and do not incorporate the strong seasonality of the site (Pérez and Romero 1992).

### 2.1.2. Statistical analyses

To statistically assess the type of relationship between seagrass abundance (as shoot density) and depth (or light attenuation), and to identify a potential threshold of abrupt transitions, we used the R package, strucchange (Zeileis et al., 2001; Andersen et al., 2009; R Core Team, 2021). We computed maximum values of shoot density associated to a given depth (or light) as input to the algorithm. Based on different fitting to linear regressions, the algorithm divides our dataset into groups. We used this to identify the potential change points between groups of data. To assess their significance, an F-statistic (Chow test) was computed as has been done for time-series and other data types (Zeileis et al., 2001; Boada et al., 2017). The F-statistics is calculated for all potential change points in the dataset interval. The null hypothesis of no structural change or unique linear regression for the whole dataset can be rejected if the F-statistic is large and crosses the boundary of the 95% level of probability (Zeileis et al., 2001).

## 2.2. Modeling approach

We developed a mechanistic plant abundance-light attenuation model (1) based on the interplay of two dynamic components, plant growth and plant mortality.

For growth, we assumed that C. nodosa abundance (as shoot density, B) follows a logistic pattern driven by the available light F(I), limited by the carrying capacity Bmax. Respiration (Rd) was included as a constant and negative contribution to growth.

Mortality was assumed to be a function of shoot density (δ(B)).

Growth and mortality were expressed on a per capita basis, resulting in

$$ \frac{dB}{dt} = \left[ F(I) \left( 1 - \frac{B}{B_{max}} \right) - R_d \right] B - \left[ \delta(B) \right] B \quad (1) $$

The relationship between growth rate and light is described by a Michaelis-Menten function (Pérez and Romero, 1992), with linear increase until saturating irradiance Km leading to a growth asymptote Fmax.

$$ F(I) = \frac{F_{max}I}{I + K_m} \quad (2) $$

Mortality follows an asymptotic decreasing function of density, assuming that, as density increases, either clonal integration or protective effects (or both) among shoots result in intraspecific facilitation, reducing shoot mortality and acting as a reinforcing feedback. We consider m0 the basal mortality, while the facilitation effect is ε, which begins to be relevant beyond a threshold of shoot density Bo, which is expressed as:

$$ \delta(B) = m_0 \frac{1}{1 + \exp^{(B-B_0)}} \quad (3) $$

Thus, combining (1), (2) and (3) we obtain:

$$ \frac{dB}{dt} = \left[ \frac{F_{max} I}{I + K_m} \left( 1 - \frac{B}{B_{max}} \right) - Rd \right] B - \left[ m_0 \frac{1}{1 + \exp^{(B-B_0)}} \right] B \quad (4) $$

Population growth rate adopted a quadratic form, with null values at zero density and at the carrying capacity density (Bmax), while population mortality rate was null at zero density, increasing sharply until a

maximum and then decreasing smoothly (Fig. 2a and b). The interaction between our two dynamic components (growth and mortality rates, equation solves) results in three points (A, B and C in Fig. 2c) where the net change of the system is zero (dB/dt = 0). To explore these three equilibrium points, we numerically solved the ODE (Ordinary Differential Equation) (4) using the deSolve package in R software (Soetaert et al., 2010; R Core Team, 2021). We qualitatively assessed the interaction between both rates and evaluated the local stability of the model using a Jacobian matrix and numerically solved equilibrium, implemented by the jacobian full function in R package rootSolve (Soetaert, 2015). Moreover, to appraise the effect of light deprivation stress, we analyzed growth and mortality rates, and their balance, for different depths. Depth has an effect to incident irradiance following an exponential and negative relation through a light extinction coefficient (a), according to the following expression (Ryther, 1956):

$$ I_h = I_0 e^{-ah} $$

where I0 is the surface irradiance and Ih the irradiance at depth h. The range of simulated depths (0–6 m) was in accordance with the real depth range in our study site. Finally, to identify the bistable region predicted by our model over the depth gradient, we conducted a bifurcation analysis using the rootSolve package in R software (Soetaert, 2014). This allowed us to identify different equilibrium points, stable or unstable, for a range of depths.

To parameterize the model and to keep it as close as possible to the real world, the values of the parameters were chosen based on the empirical data obtained here or available from literature (Table 1, Pérez and Romero, 1992; Mascaró et al., 2014).

# Results

## 3.1. Empirical approach

Shoot density of C. nodosa decreased with depth displaying a substantial variability, particularly between 0.5 and 2.5 m (Fig. 3a). The highest values (100% relative to the maximum value found, 5375 shoots m⁻²) appear between 0.5 and 1 m deep, with an abrupt decrease close to 2 m depth to values as low as 7–8% (relative to the maximum mentioned, Fig. 3a). The same data plotted against light (expressed as a percentage relative to subsurface light) show the same trend, but with even more dispersion (Fig. 3b). The Chow test indicated a significant threshold at 1.9 m depth (Fig. 3a), and 61% light (as extinction relative to subsurface irradiance), beyond which the meadow collapses to almost bare sand (Fig. 3b; and see also Fig. 3c and d for the statistics).

Shoot density before and after the threshold (mean values represented by green lines in Fig. 3a and b) shows an abrupt decrease, with a drop to near zero in the interval between 1.7 and 2.1 m deep.

## 3.2. Modeling approach

The model was numerically analyzed using the values of the parameters defined in Table 1 for a depth in accordance with a realistic depth range at our study site (0–6 m).

At 1 m depth, the interplay of these two curves results in three equilibrium points A, B and C, which correspond to values on the X-axis for C. nodosa shoot density: (BA) when both rates have null values at shoot density of 0% (eigenvalue of the Jacobian matrix, λ = −0.03); (BB) at shoot density of 24% (relative to the maximum density found, eigenvalue, λ = 0.03); and (BC) at shoot density of 88% relative to maximum density found (eigenvalue, λ = −0.06). According to these results, when seagrass meadow density is between BA and BB, density will tend to decrease to BA (zero, bare sand), given that mortality is higher than growth. In contrast, when seagrass meadow density is between BB and BC, density will tend to increase to BC (dense and stable meadow), since growth is higher than mortality. Finally, when density is
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** (a) Population growth and (b) mortality rate of *Cymodocea nodosa* as a function of shoot density at default parameters of the model (eq. (4); Table 1), and 1 m depth. (c) Equilibrium points of system (A, B and C) due to growth and mortality rates interaction, where dB/dt = 0. Bmax is the maximum shoot density, i.e. the carrying capacity, see eqs. (1) and (4).</figcaption></figure>

between B~C~ and B~max~, density will tend to decrease to B~C~, since mortality is again higher than growth. Therefore, A and C are stable points, while B is an unstable point (Fig. 4, bottom).

Simulations at different depths modify the position of equilibrium points B and C as a consequence of changes in light that diminish seagrass growth (Fig. 5a). The basin of attraction of each stable state is also modified with depth (Fig. 5b). Solving the model and determining the equilibrium points (both stable and unstable) at a variety of depths indicates three dynamic properties: (i) the vegetated stable equilibrium (Fig. 5a, filled points C) sustains lower values of shoot density at greater depths; (ii) at depths greater than 3 m, the system has only one state possible (bare sand, Fig. 5a, filled points A, Fig. 5c above dashed red line); (iii) depth increase results in an increase in the distance between the unvegetated stable state (Fig. 5a, filled points A) and the unstable point (Fig. 5a, empty points B), which corresponds to a greater basin of attraction of the bare sand state (Fig. 5b).

Based on the model results, at shallow depths (0.5–1 m) the

probability of a state shift from vegetated to bare sand is low (and, therefore, the system is more resilient), as indicated by the distance in the graph between vegetated stable point C and the unstable equilibrium point B (Fig. 5a). Seagrass resilience decreases as depth increases (Fig. 5b and c) as indicated by the increasing proximity between points B and C (Fig. 5a). According to the model, the system is bi-stable from 0.5 to 3 m, while, from 3 m depth onwards, the only possible state is the bare sand (Fig. 5c). Differences in values aside, the same bi-stable dynamic coincides with the pattern observed in Fig. 3, including (roughly) the critical depth beyond which vegetation disappears.

Discussion

Both field and modelling approaches confirm that *Cymodocea nodosa* seagrass meadows have strongly non-linear behaviours at their depth limits linked to light deprivation. Light-dependent growth and inverse density-dependent mortality supports the resilience of this system, but
<table><caption>**Table 1**
Setting of default values and dimension for model variables and parameters.</caption><thead><tr><th></th><th>Default</th><th>Unit</th><th>Description</th><th>References</th></tr></thead><tbody><tr><td>Variables</td><td></td><td></td><td></td><td></td></tr><tr><td>B</td><td></td><td>shoot m<sup>-2</sup></td><td>C. nodosa shoot density (relative value (%))</td><td></td></tr><tr><td>I</td><td></td><td>µmol s<sup>-1</sup> m<sup>-2</sup></td><td>Incident irradiance</td><td></td></tr><tr><td>Parameters</td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>0.6</td><td></td><td>Light extinction coefficient</td><td>Pérez and Romero (1992)</td></tr><tr><td>F<sub>max</sub></td><td>0.13</td><td>month<sup>-1</sup></td><td>Net shoot recruitment rate at saturating irradiance</td><td>Mascaró et al. (2014)</td></tr><tr><td>K<sub>m</sub></td><td>25</td><td>µmol s<sup>-1</sup> m<sup>-2</sup></td><td>Saturating irradiance</td><td>Pérez and Romero (1992)</td></tr><tr><td>B<sub>max</sub></td><td>100</td><td>shoot m<sup>-2</sup></td><td>Carrying capacity</td><td>Mascaró et al. (2014)</td></tr><tr><td>m<sub>0</sub></td><td>0.19</td><td>month<sup>-1</sup></td><td>Mortality rate</td><td></td></tr><tr><td>ε</td><td>0.05</td><td></td><td>Density dependent facilitation coefficient</td><td></td></tr><tr><td>B<sub>0</sub></td><td>10</td><td>shoot m<sup>-2</sup></td><td>Shoot density threshold for facilitation effect</td><td></td></tr><tr><td>Rb</td><td>0.007</td><td>month<sup>-1</sup></td><td>Respiration rate</td><td>Pérez and Romero (1992)</td></tr></tbody></table>

beyond a (remarkably shallow) depth, these stabilizing feedbacks no longer operate, and the dense seagrass meadow gives way abruptly to bare sand. The complementary approaches used in this work highlights the benefit of mixing modeling and observational studies to identify the mechanisms enhancing resilience of natural ecosystems (Scheffer and Carpenter, 2003). In our case study, *Cymodocea nodosa* meadows change abruptly from a vegetated state with high to moderate shoot densities to

an unvegetated state (bare sand) at ca. 2 m deep. At shallower depths, seagrass density displayed a high variability, from high to near zero shoots, suggesting a region of bistability typical of systems with alternative stable states. Although our sampling design did not allow us to specifically test bistability, these observations suggest bimodality and coexistence of the two alternative, potentially stable states. On the other hand, our mechanistic model, based on a simple interplay between light-dependent growth and density-dependent mortality, is consistent with these empirical observations, showing a large region of bistability. Although neither empirical field observations nor the mechanistic model are conclusive evidence on their own, together they strongly suggest the existence of a critical transition at the depth limit of the meadow.

Note that, given the strong seasonality and fluctuations of water transparency, the light extinction observations (Fig. 3b) are merely comparative between sampling points, and no great importance should be given to the 60% extinction found at the 2 m depth threshold. According to previous results on light extinction in the area (Pérez and Romero 1992), minimum light attenuation coefficients were in early summer (the time of our sampling), while the highest ones were in winter. Making a very rough estimate, we conclude that in winter only 20% of incident light (80% of light extinction) reaches the meadow at 2 m depth, and this would probably be the actual limiting factor to the meadow development. Where the depth limit of seagrass meadows lies is of considerable importance, as small retractions of this limit can imply large areal meadow losses, with subsequent losses of the associated ecosystem services. In the case of *C. nodosa* in Alfacs Bay, with a depth limit of ~2m, the meadow does not have much to retreat before it completely disappears. This may be true of many other meadows with naturally shallow vertical distributions, which are particularly prone to depth contractions due to declines in water quality. The location of the

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** *Cymodocea nodosa* shoot density (% relative to the maximum value found, 5375 shoots m⁻²) measured in each sampling site (both years, n = 119). Data are represented as a function of depth (a) and as a function of light extinction (b). In the depth gradient figure (a), bubble size indicates frequency of combined shoot density values per a given depth. Dashed lines indicating the position of the thresholds and their respective confidence intervals (95%, in grey). Green horizontal lines represent mean shoot density before and after the threshold. (c) and (d) represent the F statistic value (Chow test) obtained from the shoot density trajectory over the depth gradient and the light extinction gradient, respectively, and shows a major peak where the most significant break point occurs.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Above: Population growth (green curve) and mortality rates (red curve) of *C. nodosa* as a function of shoot density at default parameters of the model, and 1 m depth. Closed dots represent stable equilibrium points (A, C), whereas the open dot represents the unstable equilibrium point (B). $B_A$, $B_B$, $B_C$ are the values of shoot density at each equilibrium point (A, B, C, respectively), and $B_{max}$ is the maximum shoot density, i.e. the carrying capacity. The vertical red line indicates a threshold, for which, to the left (densities < $B_B$), the meadow density declines (mortality > population growth), while, to the right (densities > $B_B$) it increases (mortality < population growth). Below: representation of the basin of attraction for the two ecosystems stable states as an approximation to the meadow status resilience.</figcaption></figure>

depth limit has been used as a reliable bioindicator of water quality (Krause-Jensen et al., 2005), and identifying the factors that control the limits of distribution of seagrasses and the dynamics shaping them are highly relevant for conservation and management purposes. Although the ultimate driver of seagrass depth limit is clearly light (Duarte, 1991), our results suggest that the interplay between light and seagrasses can be more complex than previously thought. The balance between photosynthesis and respiration and light dependence has been often invoked to determine the position of the depth limit of seagrass distribution (Giddings and Eddlemon, 1978; del Barrio et al., 2014). However, if the decline in seagrass abundance with depth was only dependent on a simple and direct response to light attenuation, one would expect a gradual decrease, a much more even distribution of density values between maximum and zero, and a smoother depth limit. This is not the case in either our empirical data or in other species and/or areas (Mayot et al., 2005; Collier et al., 2007; Boudouresque et al., 2009; Enríquez

et al., 2019). Far from following a steady decline, seagrass density often abruptly shifts to bare sand after a light attenuation threshold is exceeded, suggesting the breaking of stabilizing feedback mechanisms (Maxwell et al., 2017), eventually triggering the collapse that determines the establishment of the seagrass depth limit (McGlathery et al., 2013). Our results align with other studies in which they suggest a similar mechanistic understanding of the neat shape shift typical of edges (not necessarily in the depth limit) between seagrass and bare sand habitats. These so-called spatially dependent mechanisms are based on small-scale (20–30 m) self-facilitating and self-competitive interactions that shape the landscape at kilometer-wide scales: from the formation of bare sand circles interrupting an extensive meadow, to strips and patches of vegetation over an expanse of bare sand. (Ruiz-Reynés et al., 2017).

Several processes have been described as reinforcing feedbacks to contribute to the stability of seagrasses related to light. Internally plant-driven mechanisms such as physiological, morphological and structural photoacclimation capacities can substantially increase the resilience of plants to light deprivation (Ruiz and Romero, 2001; Marín-Guirao et al., 2022). However, an important set of additional feedbacks are dependent on a minimum density to maintain the meadow environment. Specifically, when seagrasses attain a certain density, they can modify their environment, establishing the conditions for further seagrass growth by enhancing self-protection (Folkard, 2005), improving water clarity (Abal and Dennison, 1996; Carr et al., 2012), increasing nutrient retention (Maxwell et al., 2017) and, in particular, preventing sediment resuspension and erosion (Carr et al., 2012; Hansen and Reidenbach, 2013). In addition, clonal integration in these species also facilitates shoot growth (Vidondo et al. 1997; Kendrick et al. 1999, 2005; Sintes et al., 2005). We introduced this feedback in our model in the form of mortality attenuation as shoot density increases and it tracks our empirical results quite closely. Independent from its biological basis, this mechanism would result in density-dependent decrease in seagrass mortality, and, incorporated in our model, confers bistability to the system, as observed in the field.

Our mathematical model reproduces the balance between plant population growth (light dependent) and mortality (density dependent, as explained above) to reflect the collapse of the seagrass meadow at depth. Modeling results show then that the higher the distance between the unstable equilibrium point ($B_B$) and stable seagrass ($B_C$), the higher the resilience of the meadow against disturbances (see Fig. 4). For instance, if a disturbance occurs that causes mortality and then decreases density, two possibilities arise: (i) if the new density is above $B_B$, the density will recover; (ii) if the new density is below $B_B$, the meadow will shift to a bare sand state. At higher depths, this distance decreases meaning that the capacity of the meadow to maintain a high-density state is eroded. This has two implications. First, the equilibrium density ($B_C$) decreases with depth; second, the shoot density interval at which the meadow can persist also decreases (Fig. 5a), becoming less resilient (Fig. 5b). Finally, beyond a depth threshold (about 3 m), mortality is higher than growth for any density, and the dynamic system collapses to a single equilibrium point (bare sand).

The model is able to robustly reproduce the behavior of the system along a depth gradient. We acknowledge that, like all models, ours is a simplification of reality. Importantly, it does not incorporate most of the internal stabilizing feedbacks outlined in previous paragraphs, nor does it account for other potential feedbacks affecting the bare sand state, as hydrodynamics implications of the edge of the seagrass meadow (Scheffer et al., 2001; van der Heide et al., 2007; Carr et al., 2010; Maxwell et al., 2017), that would most probably increase the hysteretic behavior in the system. However, it is interesting to note that realistic complex behaviors can be predicted with relatively simple assumptions, and equally simple models (Boada et al., 2017).

Non-linearities make ecosystems predictability complex (Solé and Bascompte, 2012), rending management inherently more difficult. There have been several attempts to identify early warning indicators of
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

(caption on next column)

**Fig. 5.** (a) Population growth (green lines) at default parameters and different depths (0.5–6 m), and mortality rates (red line, note that mortality is depth-independent), both as a function of shoot density, following the model. A, B, C and $B_{max}$ as in Fig. 4. Closed dots represent stable equilibrium, and the open dot represent the unstable one. (b) Representation of the basins of attraction (green scaled) and states resilience (red and white bars) for bare sand and meadow, respectively, through depth gradient (pink arrow, on the right side). Again, stable and unstable equilibrium represented by dotted and opened dots, respectively. (c) Density values of stable and unstable equilibrium points obtained by the model solved for a range of depths (0–4m). Closed black line represents stable equilibriums of bare (A) and seagrass (B), and opened line represents the unstable equilibrium (B). Red-grey shaded area represents basin of attraction of the bare sand state, whereas blue-grey shaded one represents basin of attraction of the meadow state. Red dotted line indicates theoretical threshold for the shift from meadow to bare sand.

collapse helping to forecast at what specific level of stress transition are likely to occur (Scheffer et al., 2009; Kéfi et al., 2014; Gsell et al., 2016). However, they are generally based on time series, do not comprise a mechanistic (ecological) understanding of the processes involved, and signaling may not be exclusive to non-linear dynamics (Dakos et al., 2015). As a result, their application to management or conservation is sometimes difficult (Wouters et al., 2015; Kéfi et al., 2013; Rietkerk et al., 2021). These phenomenological approaches would benefit greatly from engaging with the underlying mechanisms that foster the resilience and trigger abrupt transitions. As a practical tool, interdisciplinary approaches to dynamic systems (combining field and modelling approaches) are much more powerful in helping determine the thresholds of stability, and in identifying if these changes represent critical transitions or not.

## 4.1. Conclusions

In summary, combining field observations with model simulations allowed us to conclude that the depth limit of *C. nodosa* distribution is an abrupt transition to bare sand, compatible with an alternative stable state dynamic, and is potentially explained by density-dependent avoidance of mortality. Given the generally clear-cut edges observed in most seagrass depth limits, it is reasonable to conclude that our results can be generalized beyond our specific case-study, although this remains to be tested. Our findings highlight the interaction between internal ecosystem feedbacks and abiotic factors, and can be a step forward in reducing unpredictability in benthic ecosystems subjected to anthropogenic stress. Current climate change scenarios paint a worrying picture of rising sea levels (Carr et al., 2012) and disrupted nutrient cycles (Havens et al., 2001; Sardans et al., 2012). These projections make the future of seagrass meadows increasingly bleak, and we should expect further decreases in depth limits and a contraction of meadows from the world's coasts. Managing meadows into the future will require a much better understanding of what drives critical transitions in seagrass ecosystems, so we can design management actions to prevent them from occurring.

## Authors contributions

**Elvira Mayol:** formal analysis (lead), data curation (lead), writing - original draft (lead), visualization (lead), investigation (equal), methodology (equal). **Jordi Boada:** methodology (equal), investigation (equal), data curation (equal), conceptualization (supporting), software (supporting), formal analysis (supporting), writing - original draft (supporting). **Marta Pérez:** project administration (lead). methodology (equal), investigation (equal), conceptualization (supporting), data curation (supporting), writing - original draft (supporting), supervision (supporting). **Neus Sanmartí:** resources (lead), methodology (equal), data curation (equal), investigation (equal), formal analysis (supporting), writing - review & editing (supporting), project administration
(supporting). **Mario Minguito-Frutos:** formal analysis (supporting), writing – review & editing (supporting). **Rohan Arthur:** conceptualization (supporting), writing – review & editing, supervision (supporting). **Teresa Alcoverro:** conceptualization (supporting), writing – review & editing (supporting), supervision (supporting), funding acquisition (supporting). **David Alonso:** software (lead), writing – review & editing (supporting). **Javier Romero:** conceptualization (lead), supervision (lead), funding acquisition (lead), methodology (equal), investigation (equal), data curation (supporting), writing – original draft (supporting), project administration (supporting).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgements

Especially grateful to the whole UMBRAL Project (CTM2013-48027-C3-1-R and CTM2013-48027-C3-3-R). We thank Guim Aguadé-Gorgorió for his theoretical advice and fruitful discussion, and Marina Almendro for her support during the laboratory procedures and the discussion of results. Biol Obrador provided the LI-COR device. Expressed gratitude to Embarcaciones Badia Alfacs for facilitation of boat renting during sampling fieldwork. DA acknowledges funding from the Spanish Ministerio de Ciencia, Innovación y Universidades and the European Regional Development Fund (ERDF) under the project CRISIS (PGC2018-096577-B-100). JB acknowledges the support received by the Spanish Ministry of Science and Innovation under the JdC fellowship (FJC2018-035566-I) and the European Commission – European Union’s Horizon 2020 MSCA – SHIFT2SOLVE-1030591.

References

Abal, E.G., Dennison, W.C., 1996. Seagrass depth range and water quality in southern Moreton Bay, Queensland, Australia. Mar. Freshw. Res. 47 (6), 763-771.

Andersen, T., Carstensen, J., Hernandez-García, E., Duarte, C.M., 2009. Ecological thresholds and regime shifts: approaches to identification. Trends Ecol. Evol. 24 (1), 49-57.

Arias-Ortiz, A., Serrano, O., Masqué, P., Lavy, P.S., Mueller, U., Kendrick, G.A., Rozaimi, M., Esteban, A., Fourqurean, J.W., Marbà, N., Mateo, M.A., Murray, K., Rule, M.J., Duarte, C.M., 2018. A marine heatwave drives massive losses from the world's largest seagrass carbon stocks. Nat. Clim. Change 8 (4), 338.

Arnold, T.M., Tanner, C.E., Rothen, M., Bullington, J., 2008. Wound-induced accumulations of condensed tannins in turtlegrass, Thalassia testudinum. Aquat. Bot. 89 (1), 27-33.

Barbier, E.B., Koch, E.W., Silliman, B.R., Hacker, S.D., Wolanski, E., Primavera, J., Granek, E.F., Polasky, S., Aswani, S., Cramer, L.A., Stoms, D.M., Kennedy, C.J., Bael, D., Kappel, C.V., Perillo, G.M.E., Reed, D.J., 2008. Coastal ecosystem-based management with nonlinear ecological functions and values. Science 319, 321-323.

del Barrio, P., Ganju, N.K., Aretxabaleta, A.L., Hayn, M., García, A., Howarth, R.W., 2014. Modeling future scenarios of light attenuation and potential seagrass success in a eutrophic estuary. Estuar. Coast Shelf Sci. 149, 13-23.

Berdugo, M., Rocamura, B.V., Solé, R., Maestre, F.T., 2021. Ecological Mechanisms Underlying Aridity Thresholds in Global Drylands.

Boada, J., Arthur, R., Alonso, D., Pagès, J.F., Pessarrodona, A., Oliva, S., Ceccherelli, G., Piazzi, L., Romero, J., Alcoverro, T., 2017. Immanent conditions determine imminent collapses: nutrient regimes define the resilience of macroalgal communities. Proc. Biol. Sci. 284 (1851), 20162814.

De Boer, W.F., 2007. Seagrass–sediment interactions, positive feedbacks and critical thresholds for occurrence: a review. Hydrobiologia 591 (1), 5-24.

Boudouresque, C.F., Bernard, G., Pergent, G., Shili, A., Verlaque, M., 2009. Regression of Mediterranean Seagrasses Caused by Natural Processes and Anthropogenic Disturbances and Stress: a Critical Review.

Burkholder, J.M., Tomasko, D.A., Touchette, B.W., 2007. Seagrasses and eutrophication. J. Exp. Mar. Biol. Ecol. 350 (1-2), 46-72.

Carr, J., D'odorico, P., McGlathery, K., Wiberg, P.L., 2010. Stability and bistability of seagrass ecosystems in shallow coastal lagoons: role of feedbacks with sediment resuspension and light attenuation. J. Geophys. Res.: Biogeosciences 115 (G3).

Carr, J.A., D'Odorico, P., McGlathery, K.J., Wiberg, P.L., 2012. Stability and resilience of seagrass meadows to seasonal and interannual dynamics and environmental stress. J. Geophys. Res.: Biogeosciences 117 (G1).

Christianian, M.J., van Belzen, J., Herman, P.M., van Katwijk, M.M., Lamers, L.P., van Leent, P.J., Bouma, T.J., 2013. Low-canopy seagrass beds still provide important coastal protection services. PLoS One 8 (5), e62413.

Collier, C.J., Lavy, P.S., Masini, R.J., Ralph, P.J., 2007. Morphological growth and meadow characteristics of the seagrass Posidonia sinuosa along a depth-related gradient of light availability. Mar. Ecol. Prog. Ser. 337, 103-115.

Connell, S.D., Fernandes, M., Burnell, O.W., Doubleday, Z.A., Griffin, K.J., Irving, A.D., Leung, J.Y.S., Owen, S., Russell, B.D., Falkenberg, L.J., 2017. Testing for thresholds of ecosystem collapse in seagrass meadows. Conserv. Biol. 31 (5), 1196-1201.

Conversi, A., Dakos, V., Gårdmark, A., Ling, S., Folke, C., Mumby, P.J., Greene, C., Edwards, M., Blenckner, T., Casini, M., Pershing, A., Möllmann, C., 2015. A Holistic View of Marine Regime Shifts, 1659, vol. 370. Philosophical Transactions of the Royal Society B: Biological Sciences, 20130279. Dynamics and Environmental Stress'. Journal of Geophysical Research: Biogeosciences 117 (G1).

Dakos, V., Carpenter, S.R., van Nes, E.H., Scheffer, M., 2015. Resilience indicators: prospects and limitations for early warnings of regime shifts. Phil. Trans. Biol. Sci. 370 (1659), 20130263.

Dennison, W.C., 1987. Effects of light on seagrass photosynthesis, growth and depth distribution. Aquat. Bot. 27 (1), 15-26.

Duarte, C.M., Sand-Jensen, B., 1990. Seagrass colonization: biomass development and shoot demography in Cymodocea nodosa patches. Marine ecology progress series. Oldendorf 67 (1), 93-103.

Duarte, C.M., 1991. Seagrass depth limits. Aquat. Bot. 40 (4), 363-377.

Duarte, C.M., Marbà, N., Krause-Jensen, D., Sánchez-Camacho, M., 2007. Testing the predictive power of seagrass depth limit models. Estuar. Coast 30 (4), 652.

Duarte, C.M., Conley, D.J., Carstensen, J., Sánchez-Camacho, M., 2009. Return to neverland: shifting baselines affect eutrophication restoration targets. Estuar. Coast 32, 29-36.

Enriquez, S., Olivé, I., Cayabyab, N., Hedley, J.D., 2019. Structural complexity governs seagrass acclimatization to depth with relevant consequences for meadow production, macrophyte diversity and habitat carbon storage capacity. Sci. Rep. 9 (1), 1-14.

Estes, J.A., Palmisano, J.F., 1974. Sea otters: their role in structuring nearshore communities. Science 185 (4156), 1058-1060.

Filbee-Dexter, K., Scheibling, R.E., 2014. Sea urchin barrens as alternative stable states of collapsed kelp ecosystems. Mar. Ecol. Prog. Ser. 495, 1-25.

Folkard, A.M., 2005. Hydrodynamics of model Posidonia oceanica patches in shallow water. Limnol. Oceanogr. 50 (5), 1592-1600.

Fourqurean, J.W., Duarte, C.M., Kennedy, H., Marbà, N., Holmer, M., Mateo, M.A., Apostolaki, E.T., Kendrick, G.A., Krause-Jensen, D., McGlathery, K.J., Serrano, O., 2012. Seagrass ecosystems as a globally significant carbon stock. Nat. Geosci. 5, 505-509.

Giddings, J., Eddleman, G.K., 1978. Photosynthesis/respiration ratios in aquatic microcosms under arsenic stress. Water Air Soil Pollut. 9 (2), 207-212.

Gsell, A.S., Scharfenberger, U., Özkundakci, D., Walters, A., Hansson, L.A., Janssen, A.B., Nöges, P., Reid, P.C., Schindler, D.E., Van Donk, E., Dakos, V., Adrian, R., 2016. Evaluating early-warning indicators of critical transitions in natural aquatic ecosystems. Proc. Natl. Acad. Sci. USA 113 (50), E8089-E8095.

Hansen, J.C., Reidenbach, M.A., 2013. Seasonal growth and senescence of a Zostera marina seagrass meadow alters wave-dominated flow and sediment suspension within a coastal bay. Estuar. Coast 36 (6), 1099-1114.

Havens, K.E., Hauxwell, J., Tyler, A.C., Thomas, S., McGlathery, K.J., Cebrian, J., Valiela, I., Steinman, A.D., Hwang, S., 2001. Complex interactions between autotrophs in shallow marine and freshwater ecosystems: implications for community responses to nutrient stress. Environ. Pollut. 113 (1), 95-107.

van der Heide, T., van Nes, E.H., Geerling, G.W., Smolders, A.J., Bouma, T.J., van Katwijk, M.M., 2007. Positive feedbacks in seagrass ecosystems: implications for success in conservation and restoration. Ecosystems 10 (8), 1311-1322.

Holling, C.S., 1973. Resilience and stability of ecological systems. Annu. Rev. Ecol. Systemat. 4 (1), 1-23.

Holmer, M., Hasler-Sheetal, H., 2014. Sulfide intrusion in seagrasses assessed by stable sulfur isotopes—a synthesis of current results. Front. Mar. Sci. 1, 64.

Kéfi, S., Dakos, V., Scheffer, M., Van Nes, E.H., Rietkerk, M., 2013. Early warning signals also precede non-catastrophic transitions. Oikos 122 (5), 641-648.

Kendrick, G.A., Eckersley, J., Walker, D.I., 1999. Landscape-scale changes in seagrass distribution over time: a case study from Success Bank, Western Australia. Aquat. Bot. 65 (1-4), 293-309.

Kéfi, S., Guttal, V., Broc, W.A., Carpenter, S.R., Ellison, A.M., Livina, V.N., Seekell, D.A., Scheffer, M., van Nes, E.H., Dakos, V., 2014. Early warning signals of ecological transitions: methods for spatial patterns. PLoS One 9 (3), e92097.

Kendrick, G.A., Duarte, C.M., Marbà, N., 2005. Clonality in seagrasses, emergent properties and seagrass landscapes. Mar. Ecol. Prog. Ser. 290, 291-296.

Kennedy, E.V., Perry, C.T., Halloran, P.R., Iglesias-Prieto, R., Schönberg, C.H., Wisshak, M., Form, A.U., Carricart-Ganivet, J.P., Fine, M., Eakin, C.M., Mumby, P.J., 2013. Avoiding coral reef functional collapse requires local and global action. Curr. Biol. 23 (10), 912-918.

Krause-Jensen, D., Greve, T.M., Nielsen, K., 2005. Eelgrass as a bioindicator under the European water framework directive. Water Resour. Manag. 19 (1), 63-75.

Ling, S.D., Scheibling, R.E., Rassweiler, A., Johnson, C.R., Shears, N., Connell, S.D., Salomon, A.K., Norderhaug, K.M., Pérez-Matus, A., Hernández, J.C., Clemente, S., Blamey, L.K., Hereu, B., Ballesteros, E., Sala, E., Garrabou, J., Cebrian, E., Zabala, M., Fujita, D., Johnson, L.E., 2015. Global regime shift dynamics of catastrophic sea urchin overgrazing. Phil. Trans. Biol. Sci. 370 (1659), 20130269.
López y Royo, C.L., Casazza, G., Pergent-Martini, C., Pergent, G., 2010. A biotic index using the seagrass Posidonia oceanica (BiPo), to evaluate ecological status of coastal waters. Ecol. Indicat. 10 (2), 380–389.

Marbà, N., Duarte, C.M., 1998. Rhizome elongation and seagrass clonal growth. Mar. Ecol. Prog. Ser. 174, 269–280.

Marín-Guirao, Lazaro, Bernardeau-Esteller, Jaime, Belando, María Dolores, Ramos-Segura, Aranzazu, Alcoverro, Teresa, Minguito-Frutos, Mario, Ruiz-Fernández, Juan M., 2022. Photo-acclimatory thresholds anticipate sudden shifts in seagrass ecosystem state under reduced light conditions. Mar. Environ. Res. 177, 105636 https://doi.org/10.1016/j.marenvres.2022.105636.

Mascaró, O., Romero, J., Pérez, M., 2014. Seasonal uncoupling of demographic processes in a marine clonal plant. Estuar. Coast Shelf Sci. 142, 23–31.

Mayot, N., Boudouresque, C.F., Leriche, A., 2005. Unexpected response of the seagrass
posidonia oceanica to a warm-water episode in the north Western Mediterranean
sea. Comptes Rendus Biol. 328 (3), 291–296.

Maxwell, P.S., Eklöf, J.S., van Katwijk, M.M., O'Brien, K.R., de la Torre-Castro, M., Boström, C., Bouma, T.J., Krause-Jensen, D., Unsworth, R.K.F., van Tussenbroek, B. L., van der Heide, T., 2017. The fundamental role of ecological feedback mechanisms for the adaptive management of seagrass ecosystems-a review. Biol. Rev. 92 (3), 1521-1538.

McGlathery, K.J., Matthew, A.R., d'Odorico, P., Fagherazzi, S., Pace, M.L., Porter, J.H., 2013. Nonlinear dynamics and alternative stable states in shallow coastal systems. Oceanography 26 (3), 220–231.

McGlathery, K.J., Sundbäck, K., Anderson, I.C., 2007. Eutrophication in shallow coastal bays and lagoons: the role of plants in the coastal filter. Mar. Ecol. Prog. Ser. 348, 1–18.

Melis, R., Ceccherelli, G., Piazzi, L., Rustici, M., 2019. Macroalgal forests and sea urchin barrens: structural complexity loss, fisheries exploitation and catastrophic regime shifts. Ecol. Complex. 37, 32–37.

Nielsen, S.L., Sand-Jensen, K., Borum, J., Geertz-Hansen, O., 2002. Depth colonization of eelgrass (Zostera marina) and Macroalgae as determined by water transparency in Danish coastal waters. Estuaries 25 (5), 1025–1032.

Nielsen, S.L., Banta, G.T., Pedersen, M.F., 2004. Decomposition of marine primary producers: consequences for nutrient recycling and retention in coastal ecosystems. In: Estuarine Nutrient Cycling: the Influence of Primary Producers. Springer, Dordrecht, pp. 187–216.

Nyström, M., Norström, A.V., Blenckner, T., de la Torre-Castro, M., Eklöf, J.S., Folke, C., Österblom, H., Steneck, R.S., Thyresson, M., Troell, M., 2012. Confronting feedbacks of degraded marine ecosystems. Ecosystems 15 (5), 695–710.

Olesen, B., Sand-Jensen, K., 1994. Patch dynamics of eelgrass Zostera marina. Mar. Ecol. Prog. Ser. 106, 147, 147.

Ontoria, Y., Gonzalez-Guedes, E., Sanmartín, N., Bernardeau-Esteller, J., Ruiz, J.M., Romero, J., Pérez, M., 2019. Interactive effects of global warming and eutrophication on a fast-growing Mediterranean seagrass. Mar. Environ. Res. 145, 27–38.

Pérez, M., Mateo, M.A., Alcoverro, T., Romero, J., 2001. Variability in detritus stocks in beds of the seagrass Cymodocea nodosa. Bot. Mar. 44 (6), 523–531.

Pérez, M., Duarte, C.M., Romero, J., Sand-Jensen, K., Alcoverro, T., 1994. Growth plasticity in Cymodocea nodosa stands: the importance of nutrient supply. Aquat. Bot. 47 (3–4), 249–264.

Pérez, M., Romero, J., 1992. Photosynthetic response to light and temperature of the seagrass Cymodocea nodosa and the prediction of its seasonality. Aquat. Bot. 43 (1), 51–62.

Pérez, M., Camp, J., 1986. Distribución Espacial Y Biomasa de Las Fanerógamas de Las Bahías Del Delta Del Ebro.

Peterson, C.H., Luettich Jr., R.A., Micheli, F., Skilleter, G.A., 2004. Attenuation of water flow inside seagrass canopies of differing structure. Mar. Ecol. Prog. Ser. 268, 81–92.

Prado, P., Tomas, F., Alcoverro, T., Romero, J., 2007. Extensive direct measurements of Posidonia oceanica defoliation confirm the importance of herbivory in temperate seagrass meadows. Mar. Ecol. Prog. Ser. 340, 63–71.

R Core Team, 2021. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL. https://www.R-project. org/.

Ramage, D.L., Schiel, D.R., 1999. Patch dynamics and response to disturbance of the seagrass Zostera novazelandica on intertidal platforms in southern New Zealand. Mar. Ecol. Prog. Ser. 189, 275–288.

Rietkerk, M., Bastiaansen, R., Banerjee, S., van de Koppel, J., Baudena, M., Doelman, A., 2021. Evasion of tipping in complex systems through spatial pattern formation. Science 374 (6564), eabj0359.

Ruiz, J.M., Romero, J., 2001. Effects of in situ experimental shading on the Mediterranean seagrass Posidonia oceanica. Mar. Ecol. Prog. Ser. 215, 107–120.

Ruiz-Reynés, D., Gomila, D., Sintes, T., Hernández-García, E., Marbà, N., Duarte, C.M., 2017. Fairy circle landscapes under the sea. Sci. Adv. 3 (8), e1603262.

Ryther, J.H., 1956. Photosynthesis in the ocean as a function of light intensity 1. Limnol. Oceanogr. 1 (1), 61–70.

Sardans, J., Rivas-Ubach, A., Peñuelas, J., 2012. The C: N: P stoichiometry of organisms and ecosystems in a changing world: a review and perspectives. Perspect. Plant Ecol. Evol. Systemat. 14 (1), 33–47.

Saunders, M.I., Leon, J., Phinn, S.R., Callaghan, D.P., O'Brien, K.R., Roelfsema, C.M., Lovelock, C.E., Lyons, M.B., Mumby, P.J., 2013. Coastal retreat and improved water quality mitigate losses of seagrass from sea level rise. Global Change Biol. 19 (8), 2569–2583.

Scheffer, M., Carpenter, S., Foley, J.A., Folke, C., Walker, B., 2001. Catastrophic shifts in ecosystems. Nature 413 (6856), 591–596.

Scheffer, M., Carpenter, S.R., 2003. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol. 18 (12), 648–656.

Scheffer, M., Bascompte, J., Brock, W.A., Brovkin, V., Carpenter, S.R., Dakos, V., Held, H., Van Nes, E.H., Rietkerk, M., Sugihara, G., 2009. Early-warning signals for critical transitions. Nature 461, 53e59'.

Sintes, T., Marba, N., Duarte, C.M., Kendrick, G.A., 2005. Nonlinear processes in seagrass colonisation explained by simple clonal growth rules. Oikos 108 (1), 165–175.

Soetaert, K., Petzoldt, T., Setzer, R.W., 2010. Solving differential equations in R: package deSolve. J. Stat. Software 33, 1–25.

Soetaert, K., 2014. Package rootSolve: Roots, Gradients and Steady-States in R. Google Scholar.

Soetaert, K., 2015. rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations [Software]. R-package version, p. 1.

Solé, R., Bascompte, J., 2012. MPB-42). In: Self-Organization in Complex Ecosystems. Princeton University Press.

Tett, P., Gowen, R., Mills, D., Fernandes, T., Gilpin, L., Huxham, M., Kennington, K., Read, P., Service, M., Wilkinson, M., Malcolm, S., 2007. Defining and detecting undesirable disturbance in the context of marine eutrophication. Mar. Pollut. Bull. 55 (1–6), 282–297.

Unsworth, R.K., Collier, C.J., Waycott, M., Mckenzie, L.J., Cullen-Unsworth, L.C., 2015. A framework for the resilience of seagrass ecosystems. Mar. Pollut. Bull. 100 (1), 34–46.

Verdura, J., Santamaría, J., Ballesteros, E., Smale, D.A., Cefalì, M.E., Golo, R., de Caralt, S., Vergés, A., Cebrian, E., 2021. Local-scale climatic refugia offer sanctuary for a habitat-forming species during a marine heatwave. J. Ecol. 109 (4), 1758–1773.

Vidondo, B., Duarte, C.M., Middelboe, A.L., Stefansen, K., Lützen, T., Nielsen, S.L., 1997. Dynamics of a landscape mosaic: size and age distributions, growth and demography of seagrass Cymodocea nodosa patches. Mar. Ecol. Prog. Ser. 158, 131–138.

Wang, C., Wang, Z.H., Sun, L., 2020. Early-warning signals for critical temperature transitions. Geophys. Res. Lett. 47 (14), e2020GL088503.

Wernberg, T., Bennett, S., Babcock, R.C., De Bettignies, T., Cure, K., Depczynski, M., Dufois, F., Fromont, J., Fulton, C.J., Hovey, R.K., Harvey, E.S., Holmes, T.H., Kendrick, G.A., Radford, B., Santana Garcon, J., Saunders, B.J., Smale, D.A., Thomsen, M.S., Tuckett, C.A., Tuya, F., Vanderklift, M.A., Wilson, S., 2016. Climate-driven regime shift of a temperate marine ecosystem. Science 353 (6295), 169–172.

Wouters, N., Dakos, V., Edwards, M., Serafim, M.P., Valayer, P.J., Cabral, H.N., 2015. Evidencing a regime shift in the North Sea using early-warning signals as indicators of critical transitions. Estuar. Coast Shelf Sci. 152, 65–72.

Zeileis, Achim, Leisch, Friedrich, Kurt, Hornik, Kleiber, Christian, 2001. Strucchange. An R Package for Testing for Structural Change in Linear Regression Models.