High Dimensionality of the Stability of a Marine Benthic Ecosystem

Stability is a central property of complex systems and encompasses multiple dimensions such as resistance, resilience, recovery, and invariability. How these dimensions correlate among them is focus of recent ecological research, but empirical evidence at regional scales, at which conservation decisions are usually made, remains absent. Using a field-based manipulative experiment conducted in two marine intertidal regions, we analyze the correlations among different aspects of stability in functioning (community cover) and composition of local communities facing a press disturbance. The experiment involved the removal of the local space-dominant species for 35 months in eight sites under different environmental regimes in northern- and southern-central Chile (ca. 30 and 40°S, respectively). After the disturbance, the magnitude of the initial responses and the recovery patterns were similar among communities dominated by different species, but varied between the functional and compositional response variables, and among four dimensions of stability. The recovery trajectories in function and composition remained mostly uncorrelated across the system. Yet, larger initial functional responses were associated with faster recovery trajectories—high functional resilience, in turn, was associated with both, high and low variability in the pattern of recovery. Finally, the compositional stability dimensions were independent from each other. The results suggest that varying community compositions can perform similar levels of functioning, which might be the result of strong compensatory dynamics among species competing for space in these communities. Knowledge of several, and sometimes independent, aspects of stability is mandatory to fully describe the stability of complex ecological systems.

Introduction

Anthropogenic environmental change is disrupting the functioning and composition of ecosystems over multiple spatiotemporal scales (De Laender et al., 2016; Isbell et al., 2017). Understanding the capacity of ecosystems to absorb and respond to such disturbances is challenging, because stability is a multifaceted property characterized by several dimensions (Pimm, 1984; Ives and Carpenter, 2007; Donohue et al., 2013). Recent theoretical and empirical advances provide the framework to assess this dimensionality through correlations between multiple aspects of stability such as resistance, resilience, recovery, and invariability (Donohue et al., 2013; Hillebrand et al., 2018; Radchuk et al., 2019). Most scientific research to date has focused on either one or two dimensions, which has impaired our capacity to fully understand the stability of complex ecological systems (Donohue et al., 2016; Kéfi et al., 2019). In addition, the functional consequences of biodiversity loss can be exacerbated at regional spatial scales, which are relevant for conservation or policy (Isbell et al., 2017). To our best knowledge, however, the correlations of multiples dimensions of stability have not been assessed at large spatial scales yet (Gonzalez et al., 2020).

Each dimension of stability can be assessed either for aggregate ecosystem properties performed by the whole community, such as productivity, community biomass, or cover (functional stability), or for the combination of species identities and abundances in the community (compositional stability; Micheli et al., 1999). The relationship between both domains of stability sheds light into the contribution of particular combinations of species to the aggregate function. For example, if the community is characterized by asynchronous species dynamics, in which the decay of some species is compensated by proportional increases of others, then the community could maintain relatively stable aggregate properties over time (Yachi and Loreau, 1999; Gonzalez and Loreau, 2009). This implies that a low compositional stability can sustain a high functional stability (Tilman, 1996; Micheli et al., 1999). In these cases, negative or null correlations between functional and compositional aspects of stability can be observed (e.g., Hillebrand et al., 2018; Hillebrand and Kunze, 2020). The stabilizing role of such asynchronous or compensatory species dynamics is supported by several studies (Gonzalez and Loreau, 2009; Valencia et al., 2020). However, it is less clear how functional and compositional stabilities will be correlated in communities facing the widespread chronic disturbances experienced by ecosystems (Guelzow et al., 2017).

Within the functional and compositional domains, the dimensions of stability include, but are not limited to, resistance, or the degree to which an aggregate function or species composition remains unchanged after a disturbance, resilience, or the speed of community recovery after the disturbance, recovery, or the degree to which the community recovers from the disturbance, and invariability, or the degree of constancy of the function or composition around the temporal recovery trend (Grimm and Wissel, 1997; Hillebrand et al., 2018). These aspects of stability can be independent of each other or can be correlated. Several models suggest that most stability dimensions are positively correlated (e.g., Scheffer et al., 2009; Radchuk et al., 2019). For instance, communities with low invariability (i.e., highly variable functioning or composition over time) can show low resistance against disturbances, because of the increased risk of extinctions of species with small population sizes (Lande et al., 2003). In those cases, stability is unidimensional and conservation decisions could maximize both aspects of stability, i.e., enhancing resistance would also enhance invariability. However, negative correlations among dimensions can be also observed. Resistance, for example, can be negatively related with resilience when disturbance has a strong impact on the community and intrinsic growth rates are high: Larger initial disturbances can allow for faster recovery trajectories in function or composition (Harrison, 1979; Hillebrand et al., 2018). Negatively correlated stability dimensions imply that policy trade-offs would be needed to maintain either one or another aspect of stability. Finally, chronic disturbances, like the loss of dominant species, can lead to uncorrelated dimensions of stability (Donohue et al., 2013). The lack of correlation among stability dimensions (i.e., a high dimensionality) indicates that multiple aspects need to be quantified simultaneously to fully characterize the overall stability of the ecosystem.

In this study, we analyze the results of a field-based manipulative experiment in which a press disturbance, i.e., the chronic extinction of the locally dominant species, was simulated to assess the dimensionality of stability within the functional and compositional domains of rocky intertidal sessile communities. To this end, we used total community cover as a relevant metric of ecosystem functioning. For communities of sessile species, this function expresses the degree to which competing species use the space as a limiting resource and thus, its relevant to understand both functional and compositional stability (Paine, 1980; Dunstan and Johnson, 2006; Stachowicz et al., 2008). Multiple sites spanning a large environmental gradient (Figure 1) were used to assess the dimensionality of stability at an ecologically relevant spatial scale (Isbell et al., 2017). Our experimental design involved that the locally dominant species varied among sites, which allowed us to enhance the generality of our findings (e.g., Arnillas and Cadotte, 2019). These dominant species were functionally and evolutionary distinct and included the red corticated alga Mazzaella laminarioides, chthamalid barnacles (mixture of Jehlius cirratus and Notochthamalus scabrosus), and the purple mussel Perumytilus purpuratus (hereafter referred to as Mazzaella, barnacles, and Perumytilus, respectively; Table 1). This setup allowed us to test the following three hypotheses.

FIGURE 1

Figure 1. Location of study sites in the northern and southern regions in the SE Pacific shore of Chile. Color scales denote long-term temporal mean (2003–2018) of (A) sea surface temperature and (B) log₁₀ Chlorophyll-a concentration using L3 MODIS-aqua data (Lara et al., 2019). Site codes are in Table 1. Mazzaella-dominated communities: LIMA and PTAL; barnacle-dominated communities: TEMB, GUAN, and CHAI: Perumytilus-dominated communities: CHEU, CALF, and PUCA. Sea surface temperature and Chlorophyll-a data were not included in the statistical analyses and are shown here to describe the environmental conditions across the province.

TABLE 1

Table 1. Identity of the dominant species experimentally removed at each rocky shore site in each region.

H1: The magnitude of the stability response will depend on the identity of the removed dominant species.

The three dominant species use the settlement substratum in different ways: Perumytilus is able to monopolize almost the 100% of primary substratum, while the barnacle complex occurs in combination with bare rock and filamentous algae (Table 1; Broitman et al., 2001; Moreno, 2001; Lagos et al., 2005). Mazzaella, on the other hand, uses a smaller proportion of the primary substratum, but generates a large canopy-like structure above the rock (Hoffmann and Santelices, 1997). Thus, the removal of these species should trigger different responses of the ensuing assemblages (see also next paragraph, below).

H2: Across sites, the removal of the dominant species will elicit negative or neutral correlations between functional and compositional stability dimensions due to prevailing compensatory dynamics.

The removal of Perumytilus triggers the compensation between subordinate taxa: fast-growing ephemeral macroalgae such as Ulothrix sp. and Pyropia spp. dominate the assemblage during the first 2 months after the removal, then they are replaced by barnacles and Ulva sp., and eventually by Mazzaella (Aguilera and Navarrete, 2012; Tejada-Martinez et al., 2016). Barnacles facilitate small grazers; therefore, the removal of barnacles leads to increased abundances of microalgae but decreased abundances of perennial macroalgae (Harley, 2006; Mrowicki et al., 2014; Whalen et al., 2016). Finally, the removal of Mazzaella is associated with the decay of the red alga Gymnogongrus furcellatus (Asterfilopsis furcellata), but the increase in the abundance of Ulva rigida, Lithophyllum spp., and the ephemerals Petalonia fascia and Scytosiphon lomentaria (Jara and Moreno, 1984)—these dynamics are ascribed to shading effects of Mazzaella on the understorey (Moreno and Jaramillo, 1983; Nielsen and Navarrete, 2004).

H3: Across sites and within the stability domains (functional and compositional), the stability dimensions will be correlated with each other.

According to previous modeling work (Radchuk et al., 2019), the stability dimensions will show positive pairwise correlations. The exception will be negative correlations of resistance-resilience, because stronger impacts can allow faster recovery rates, and resilience-invariability, because a faster recovery is expected to be associated with more variable temporal trajectories in function and composition (Hillebrand et al., 2018).

Materials and Methods

Study Province

The study was conducted at multiple field sites on the extremes of a large transitional biogeographic province that extends between the Peruvian and Magellan provinces along the Chilean sector of the South East Pacific (Camus, 2001; Thiel et al., 2007; Lara et al., 2019). We used eight study sites that where evenly split between the northern and southern extremes of the province (around 30.23°S and 39.9°, respectively)—each region encompassing ca. 200 km of shoreline (see Figure 1). Therefore, the experimental design considered four spatial scales: The province encompassed two regions; each region included four sites. Then, each site contained several replicated plots. A plot was a fixed location where the abundance of macrobenthic sessile species was estimated over time. There were two types of plots, either control or disturbed in which the dominant species was removed. Within each site, each stability metric was calculated for randomly paired plots (one of each type of plot). Depending on the tested hypothesis, the analyses were conducted at the scale of individual plots, paired plots, or after aggregating these data according to the identity of the removed species, which was a site-level attribute.

The northern region is characterized by semi-permanent coastal upwelling activity, which coincides with a significant transition in abiotic (Hormazabal et al., 2004; Garreaud et al., 2011; Tapia et al., 2014) and biotic conditions (Valdivia et al., 2015). Compared with neighbor locations, this upwelling point has been associated with more persistent abiotic environmental conditions, more synchronous species fluctuations, and lower invariability in total cover (Valdivia et al., 2013).

The southern region presents significant correlations between Chlorophyll-a concentration, riverine inputs, and runoff from forestry and agricultural activities (Van Holt et al., 2012; Lara et al., 2019). In this region, upwelling activity is weaker and more intermittent than the sector equatorward from Concepción, which is dominated by strong seasonal upwelling (ca. 37 S; Atkinson et al., 2002; Letelier et al., 2009; Pinochet et al., 2019).

Rocky-shore communities inhabiting the mid-intertidal zone across the province are composed of a diverse set of mobile and sessile species, where human harvesting and environmental context modulate ecological interactions (Broitman et al., 2001, 2011; Lagos et al., 2005; Caro et al., 2010). The sessile communities are dominated by green macroalgae like Ulva compressa and U. rigida, red macroalgae like Mazzaella laminarioides and Gelidium chilense, or filter feeders like the barnacles Jehlius cirratus, Notochthamalus scabrosus, and Perumytilus purpuratus (reviewed in Aguilera et al., 2019).

Experimental Design and Setup

At each site, we conducted a manipulative experiment with dominant-species removal as a fixed factor. The experiments were conducted between October 2014 and August 2017 in the southern region, and between December 2014 and September 2017 in the northern region—in both cases, the prolonged time span (ca. 3 years) allowed us to fully capture the seasonal and intra-annual variation in abiotic conditions and species composition. Prior to the start of the experiment, we estimated relative species abundances in ten 30 × 30 cm plots haphazardly located on the mid-intertidal zone to determine locally dominant species as those with the highest relative abundances at each site (see also Valone and Balaban-Feld, 2018). At each site, quadrats were located within an alongshore transect (ca. 3 m wide) covering ca. 50 m of the coastline. The red corticated alga Mazzaella laminarioides and chthamalid barnacles (a mixture of Jehlius cirratus and Notochthamalus scabrosus) were the dominant taxa in the northern region; chthamalid barnacles and the purple mussel Perumytilus purpuratus were the dominant species in the southern region (hereafter referred as Mazzaella, barnacles, and Perumytilus, respectively; Table 1). Before the manipulations, mean percentage cover of dominant species ranged from ca. 30% for Mazzaella in the northern region to 100% for Perumytilus in the southern region. Simpson’s concentration index varied between 0.4 and 0.9 in the northern and southern region, respectively, picturing a higher evenness in species abundance distribution in the former than the latter region (Table 1).

Following site characterization, we started identical field experiments. We used 30 × 30 cm experimental plots on the mid-intertidal zone of each site, which were permanently marked with stainless-steel bolts. The size of the experimental unit was selected to minimize the impact on the local intertidal communities. For logistic reasons, we deployed 10 replicate plots in the sites located on the northern sector and 20 replicates in the southern sites. Experimental plots were haphazardly located within areas of high abundance of the locally dominant species and relatively flat and horizontal rocky surfaces. Thus, we avoided crevices, tide-pools, and vertical surfaces to reduce biotic variation related to local spatial heterogeneity.

At the onset of experimental manipulations (i.e., October and December 2014 in the southern and northern regions, respectively), we haphazardly selected half of the plots at each site and completely removed the dominant species with scrapers and chisels. The remaining half of the experimental units remained un-manipulated and assigned to the control level. Settlers and recruits of the dominant species were further removed approximately every 3 months, i.e., the press disturbance treatment (see also Bulleri et al., 2012). Prior field studies on recruitment patterns in the province suggest that the 3-month window between manipulations allows to functionally exclude dominant species (Aguilera et al., 2013; Valdivia et al., 2013, 2015).

All plots were sampled immediately before the removal of local dominants, 1–2 months after the manipulation, and then approximately every 3 months until the end of the experiment. All macroalgae and sessile macro-invertebrates (>5 mm) occurring in each plot were identified in situ and classified to the lowest possible taxonomic level (usually species). For each plot, we used a 0.09-m² frame, divided in 25 equal fields with a monofilament line, to estimate sessile species abundances as percentage cover (1% resolution). We quantified the organisms growing on both primary and secondary substrata (i.e., rock and other organisms, respectively), which allowed total community cover to increase beyond 100%. This protocol has been used in several studies of benthic diversity along Chile and elsewhere (see Broitman et al., 2011; Bulleri et al., 2012). The data obtained during the 11 successive surveys, conducted after the first removal, were used to estimate four dimensions of stability, i.e., resistance, resilience, recovery, and invariability. Since we were interested in the response of the remaining community to the local extinction of dominants, the latter were removed from the dataset when we assessed community stability.

Species Compensation

The degree of species compensation was measured by assessing community-wide synchrony in species abundances. To this aim, we used the φ statistic (Loreau and de Mazancourt, 2008):

$\begin{array}{c}\hfill \mathrm{\phi }=\frac{{\mathrm{\sigma }}_F^2}{{\left({\sum }_{i=1}^N{\mathrm{\sigma }}_{f_i}\right)}^2}\hfill \end{array}$

$\begin{array}{c}\hfill F=\sum _{i=1}^N{f}_i,\hfill \end{array}$

where f_i is the percentage cover of species i, F is the equivalent aggregate community-level function (total community cover), and ${\mathrm{\sigma }}_{f_i}^2$ and ${\mathrm{\sigma }}_F^2$ are their respective temporal variances. The index is standardized between 0 (perfect asynchrony, i.e., most species’ abundances are negatively correlated over time) and 1 (perfect synchrony, i.e., most species’ abundances are positively correlated). Since φ is independent of the magnitude and distribution of species abundances and variances, it allows us to quantitatively compare communities with different numbers of species (Loreau and de Mazancourt, 2008).

Stability Dimensions

Stability was assessed according to the metrics for resistance, resilience, recovery, and invariability described by Hillebrand et al. (2018). Each metric was estimated by randomly pairing one control and one disturbed plot in each site. In this way, we obtained replicate estimations within each site. The stability metrics in the functional domain were calculated on the basis of total community cover (F), and the stability metrics in the compositional domain were calculated on the basis of Bray-Curtis similarities (sim) between the disturbed and control plots (C).

Resistance (a) was estimated for function as the log-response ratio (LRR) calculated after the disturbance took place (1–2 months after the removal of local dominants):

$\begin{array}{c}\hfill a=ln\left(\frac{F_{dist}}{F_{con}}\right),\hfill \end{array}$

where F_dist is total community cover in the dominant-removed plot, and F_con is the same property in the control plot. A value of a = 0 indicates maximum resistance of the community to the disturbance, a < 0 indicates low resistance due to underperformance, and a > 0 low resistance due to overperformance relative to controls, respectively (Hillebrand et al., 2018). This means that both initially increased and decreased functioning will reflect low functional resistance. For that reason, we used the absolute value of resistance in some hypothesis tests as described below (see section “Statistical analyses”: H1 and H2).

The resistance of composition (a) was calculated on the basis of Bray-Curtis similarities between the disturbed and the control plots:

$\begin{array}{c}\hfill a=\left(\frac{C_{dist}}{C_{con}}\right),\hfill \end{array}$

Where C_dist and C_con are the species-abundance matrices in the dominant-removed and control plots, respectively. Maximum similarities of a = 1 indicate 100% of similarity between disturbed and control plots and thus maximum resistance. Minimum similarities of a = 0 indicate 0% of similarity between both conditions and thus extremely low resistance.

Resilience (b) was calculated for each plot as the slope of the linear regression of LRR calculated for each sampling event (Hillebrand et al., 2018):

$\begin{array}{c}\hfill ln\left(\frac{F_{dist}}{F_{con}}\right)=i+b\times t,\hfill \end{array}$

$\begin{array}{c}\hfill dis\left(\frac{C_{dist}}{C_{con}}\right)=i+b\times t\hfill \end{array}$

where i, b, and t are the intercept, slope, and time, respectively. A b-value equaling zero indicates no recuperation over time, b > 0 indicates faster recuperation in case of a negative initial effect of disturbance on the function or composition, and b < 0 indicates further deviation from control in case of negative initial effects. When the disturbance generates positive effects on functioning (i.e., a > 0), b > 0 is interpreted as further deviation from control and vice versa.

Recovery (c) was estimated as the LRR and similarity calculated at the end of the experiment. c-values equal to zero denote maximum recovery, c < 0 incomplete recovery, and c > 0 overcompensation relative to control. Like for resistance, both c < 0 and c > 0 reflect low functional recovery; thus, the absolute value of functional c was analyzed in the test of hypotheses H1 and H3 (see section “Statistical analyses”).

Invariability (d) was calculated as the multiplicative inverse of the standard deviation of residuals (resid) from the regression slope b as:

$\begin{array}{c}\hfill d=\left(\frac{1}{sd\left(resi{d}_b\right)}\right)\hfill \end{array}$

Larger values de d denote smaller temporal variations around the trend of recovery in function or composition (Hillebrand et al., 2018). Note that invariability is termed as Temporal stability by Hillebrand et al. (2018). However, we have preferred to use Invariability to avoid confusions with the term Stability (see also Radchuk et al., 2019).

Statistical Analyses

For simplicity, hereafter the community in each site will be referred to according to the identity of the locally dominant species that was removed (i.e., Mazzaella-, barnacle-, or Perumytilus-dominated communities). This means that the plots and sites will be grouped into each of the three community types in some of the statistical analyses.

H1: The magnitude of the stability response will depend on the identity of the removed dominant species

To test this hypothesis, we first used t-tests for each community and within each stability domain (either functional or compositional) to determine if resistance, resilience, and recovery significantly differed from the reference values. Secondly, we used separate general lineal models (LMM) to compare the magnitude of each stability dimension among the three communities. Accordingly, we used the absolute value of functional resistance and recovery to focus on their magnitudes. The models included the identity of the removed dominant species as fixed effect, and site as random effect. Nevertheless, site was removed from the models of the compositional domain due to singularity (i.e., variance nearly zero).

H2: Across sites, the removal of the dominant species will elicit negative or neutral correlations between functional and compositional stability dimensions due to prevailing compensatory dynamics

First, we used one-tailed t-tests to determine whether φ was significantly lower than 1 (i.e., perfect synchrony in species fluctuations) for each community. In addition, we used a generalized linear mixed model (GLMM) to determine if the removal of the locally dominant species decreased φ (i.e., increased the degree of species compensation within the community). For the model we used a Gamma structure of errors and a log link function; the model included the identity of the dominant species in each site and the removal treatment (removal or control) as crossed and fixed factors, and the site (eight levels) as random factor. After initial the fit, however, site was removed due to singularity.

Second, Pearson product-moment correlations were calculated between functional and compositional dimensions of stability. Correlations were calculated separately for each community and stability dimension. In these correlations, we used the absolute value of every functional stability dimension to assess the magnitude, rather than the direction, of these responses. The absence of statistically significant correlations, in addition to near-zero values of φ were used as evidence of species compensation after the removal of dominant species.

Third, to identify the species driving compensation (or the lack thereof), we used Similarity Percentage Routines (SIMPER; Clarke, 1993) that partitioned the contribution of each species to the between-treatment dissimilarity in each community.

H3: Across sites and within the stability domains (functional and compositional), the stability dimensions will be correlated with each other

Pearson product-moment correlations were used to test this hypothesis. Within the functional and compositional domains, analyses were conducted combining all sites and communities to test for a tendency between these dimensions to correlate independently of the ecological context and species identities.

All analyses, calculations, and plots were done in R programming environment version 4.0.2 and the R-packages lme4, MuMIn, sjPlot, tidyverse, and vegan (Bates et al., 2015; Oksanen et al., 2019; Wickham et al., 2019; Bartoń, 2020; Lüdecke, 2020; R Core Team, 2020), except for Figure 1, which was made using MatLab R2014a and data from MODIS-Aqua.

Results

H1: The magnitude of the stability response will depend on the identity of the removed dominant species

The different communities significantly responded to the experimental disturbances within both functional and compositional domains (Figures 2A–H, respectively; Supplementary Appendix Table 1). The exception was the functional resilience of the barnacle- and Perumytilus-dominated communities, which did not differ from zero (Supplementary Appendix Table 1 and Figure 2C). When we compared the magnitude of the responses among communities, the three communities exhibited remarkably similar stability responses to the experimental disturbance for both domains (Figure 2). Accordingly, we did not detect statistically significant differences among the communities, with the exception of functional resilience (GLMM contrasts in Supplementary Appendix Table 2: Mazzaella > barnacles; Figure 2C).

FIGURE 2

Figure 2. Four dimensions of stability in the (A,C,E,G) functional and (B,D,F,H) compositional domains. Note that (A) functional resistance and (E) recovery are shown as absolute values in order to compare magnitudes, rather than directions of effects. Values are mean ± 95% confidence intervals. Dashed lines represent maximum resistance and recovery; solid gray lines represent minimum resilience. All stability dimensions statistically differed from reference values, excepting (C) functional resilience of barnacle- and Perumytilus-dominated communities (t-tests: t₁₆ = −1.24, P > 0.05; t₂₈ = 1.28, P > 0.05, respectively). Accordingly, the three types of communities exhibited similar stability responses, excepting functional resilience (GLMM contrasts: Mazzaella > barnacles; Mazzaella > Perumytilus).

H2: Across sites, the removal of the dominant species will elicit negative or neutral correlations between functional and compositional stability dimensions due to prevailing compensatory dynamics

The low values of the proxy for species synchrony, φ, indicated that species compensation was common across communities and removal treatments (Figure 3). The parameter was on average 0.11 (0.06–0.16), 0.30 (0.21–0.40), and 0.17 (0.12–0.21) for Mazzaella, barnacles, and Perumytilus communities, respectively (mean [95% CI]). These values were significantly lower than 1 in each community (one-tailed t-tests: t₉ = −42.2, t₁₆ = −17.5, and t₂₈ = −38.2 for Mazzaella, barnacles, and Perumytilus, respectively; P < 0.001). However, the response of these dynamics to the dominant removal varied among communities. The press disturbance increased the strength of species compensation in the Perumytilus communities, but not in the others (GLM: R² = 0.45, P < 0.001 for an 80-% decrease in φ in the Perumytilus communities, see Figure 3). In these communities, the strengthening of species compensation was accounted for by the overshoot in the abundance of Ulvoids (mixture of Ulva rigida and U. compressa) during the first half of the experiment (late 2014 – early 2016; see Supplementary Appendix Figure 1); this taxon was replaced by chthamalid barnacles over time, which remained relatively stable until the end of the experiment in mid-2017 (Supplementary Appendix Figure 1).

FIGURE 3

Figure 3. Community-wide synchrony among subordinate species in each community type (Mazzaella-, barnacle-, or Perumytilus-dominated community). The range of φ is [0, 1], with 1 = perfectly synchronous species fluctuations and 0 = perfectly asynchronous, compensatory species fluctuations. Values are mean ± 95% confidence intervals. All communities statistically differed from 1 (one-tailed t-tests: t₉ = −42.2, P < 0.001; t₁₆ = −17.5, P < 0.001; t₂₈ = −38.2, P < 0.001, respectively), indicating prevalent compensatory dynamics across the region. The removal of the dominant species systematically increased the strength of species compensation in the Perumytilus communities, but not in the others (Gamma GLM, R² = 0.45, P < 0.001 for Perumytilus).

About seven (TEMB and GUAN) and ten (PTAL) taxa accounted for ca. 90% of the dissimilarities between treatments in the northern region; about three (CALF) and five (CHAI) taxa explained the 90% of the multivariate differences between treatments in the southern metacommunity (Supplementary Appendix Table 3). Across both regions, these species included encrusting, filamentous, and corticated macroalgae, in addition to barnacles (e.g., Austromegabalanus psittacus; Supplementary Appendix Table 3). The lichen Verrucaria sp. accounted for a relatively important percentage of between-treatment dissimilarities in the northern region (TEMB, GUAN, and PTAL; Supplementary Appendix Table 3).

According to the high prevalence of species compensation, we observed almost no statistically significant correlation between functional and compositional stabilities across communities (Figure 4). Separately for the three community types, we observed non-statistically significant, correlations for resistance (Figure 4A and Supplementary Appendix Table 4). For resilience, only the communities in which barnacles were removed showed a significant, and negative, correlation between the functional (absolute value) and compositional domains (green symbols in Figure 4B and Supplementary Appendix Table 4). Similarly, the barnacle-dominated communities exhibited a negative correlation between functional (absolute value) and compositional recoveries (green symbols in Figure 4C and Supplementary Appendix Table 4). Finally, functional invariability was uncorrelated from compositional invariability (Figure 4D and Supplementary Appendix Table 4).

FIGURE 4

Figure 4. Pairwise associations between functional and compositional (A) resistance, (B) resilience, (C) recovery, and (D) invariability. In these analyses, all functional dimensions were used as absolute values in order to correlated magnitudes, rather that directions of effects. Dashed lines represent maximum resistance and recovery; solid gray lines represent minimum resilience. Only (B) resilience and (C) recovery of the barnacle-dominated communities (green symbols and trend lines ± 95% confidence intervals) showed significant function-composition correlations (r = −0.52 and r = −0.57, respectively; P < 0.05).

H3: Across sites and within the stability domains (functional and compositional), the stability dimensions will be correlated with each other

Our results did not support this hypothesis, yet we detected some significant correlations (Figure 5 and Supplementary Appendix Table 5). In the functional domain, statistically significant correlations were found between resistance and resilience (negative; Figure 5A and Supplementary Appendix Table 5), resistance and invariability (positive; Figure 5D and Supplementary Appendix Table 5), and resilience and invariability (negative; Figure 5E and Supplementary Appendix Table 3). In the compositional domain, we did not detect statistically significant correlations among stability dimensions (Figure 6 and Supplementary Appendix Table 5).

FIGURE 5

Figure 5. Associations among four dimensions of functional stability. All pairwise combination were analyzed: resilience against (A) resistance; recovery against (B) resistance and (C) resilience; and invariability against (D) resistance, (E) resilience, and (F) recovery. Dashed lines represent reference values for maximum resistance and recovery; solid gray lines represent reference values for minimum resilience. Across communities, statistically significant correlations were detected for the resistance-resilience, resistance-invariability, and resilience-invariability pairs (trend lines ± 95% confidence intervals; r = −0.38, r = 0.39, and r = −0.36, respectively; P < 0.01).

FIGURE 6

Figure 6. Associations among four dimensions of compositional stability. All pairwise combination were analyzed: resilience against (A) resistance; recovery against (B) resistance and (C) resilience; and invariability against (D) resistance, (E) resilience, and (F) recovery. Dashed lines represent maximum resistance and recovery; solid gray lines represent minimum resilience. No statistically significant correlation was detected across communities.

Discussion

This study showed that the stability of rocky intertidal communities in a broad biogeographical province of the SE Pacific shore is a multidimensional property. After the removal of locally dominant species, the magnitude of the effect of the disturbance and the temporal pattern of recovery were surprisingly similar among communities but varied between functional and compositional response variables, and among the four dimensions of stability examined here; namely, resistance, resilience, recovery, and invariability. Only in the functional domain we still observed significant correlations between resistance and resilience (negative), resistance and invariability (positive), and resilience and invariability (negative). The high prevalence of species compensation and asynchronous dynamics could explain the general lack of correlation between functional and compositional stability domains, indicating that functional recovery may not involve a fully compositional recovery in this study system. Our results are in line with recent evidence of the multifaceted nature of ecological stability (Donohue et al., 2016; Hillebrand et al., 2018; Pennekamp et al., 2018; Kéfi et al., 2019; Radchuk et al., 2019), which we further discuss below.

Different Communities Exhibited Similar Initial and Long-Term Stability Responses to Local Extinction

The disturbance elicited significant stability responses within both the functional and compositional domain regardless of the dominant species identity, with the exception of functional resilience, which was indistinguishable from zero in the case of the invertebrate-dominated communities. Despite these differences, no community fully recovered in the end of the experiment, 35 months after the start of the disturbance. Importantly, the consistently negative compositional resilience provides further support for a persistent departure of these assemblages from the reference state (i.e., control plots), as shown in phytoplankton mesocosm experiments (Hillebrand et al., 2018). This result may indicate that the recovery, and restauration, of this ecosystem may need a longer time span than here considered. Long-term (i.e., 20 years) monitoring of rocky intertidal communities suggests the existence of cyclic succession after severe disturbances, with the alternation of dominant species every ca. 2 years (Benincà et al., 2015). These non-equilibrium dynamics can well be amplified by environmental stochasticity, as shown by experimental removals of dominant mussels (Wootton, 2010). Thus, the non-equilibrium and near-chaotic nature of coastal systems may obstruct our ability to restore disturbed coastal ecosystems. Indeed, a major review of restoration strategies indicates that functional and compositional recovery of coastal systems may reach on average a 34% in 12 years, even when active restoration is enforced (Moreno-Mateos et al., 2020).

Between-site environmental differences could have generated “unexplained” variability in the stability responses, precluding our ability to detect differences among the dominant species. Nevertheless, our analyses suggest that the factors (other than dominant’s identity) associated to between-site differences may have little influence on the stability measurements (compare the pseudo-R² between fixed factors and full models in Supplementary Appendix Table 2). In the case of functional stability, between-site differences accounted for small, but non-zero, amount of variation in the stability responses. Dampened environmental variability in sites under the influence of persistent upwelling could have influenced the invariability in community cover (Valdivia et al., 2013). Sea water temperature patterns define the phenology of most intertidal organisms (Sanford, 2002; Baldanzi et al., 2018), which, mediated through larval dispersal processes, can influence overall population growth rates, and thus functional resilience and recovery (Svensson et al., 2005). In the case of compositional stability, on the contrary, site accounted for almost no variation in Bray-Curtis dissimilarity. If site differences seemly accounted for a small variation in stability responses, then, to what degree the stability observed at the local scales can scale up to regional and seascape scales?

The area-invariability relationship indicates that stability should increase from local to regional scales due to asynchronous dynamics among sites within the landscape, and also among species within sites (Wang and Loreau, 2014; Zhang et al., 2018). Moreover, increased beta diversity (i.e., between-site differences in species composition) is expected to increase the spatial asynchrony in functioning among sites, and thus the stability of the whole landscape (Wang and Loreau, 2016). According to the high prevalence of asynchronous species dynamics (this study; Valdivia et al., 2013) and high beta diversity (Broitman et al., 2001; Valdivia et al., 2015, 2017), predictions based on the area-invariability relationship should hold in our study system. Albeit the extension of this theory to other stability dimensions and the compositional domain still needs further development, we suggest that maximizing protected areas, rather than focusing on patches that contain unique assemblages, should be a conservation goal to maintain the functioning in the analyzed ecosystem (see also Economo, 2011; Ospina-Alvarez et al., 2020).

Compensatory Dynamics Were Common Within the System, Leading to Null or Negative Function-Composition Correlations

In agreement with our prediction, asynchronous and compensatory species dynamics were strong in the study region. The strength of these dynamics, however, varied among the three analyzed community types. In particular, the removal of Perumytilus reduced the synchrony in species abundances following compensatory dynamics between ulvoid algae and barnacles, according to the facilitation model of ecological succession (Connell and Slatyer, 1977). In addition to competitive species replacement during succession, asynchronous species dynamics can also reflect differential species environmental responses, population compensatory cycles, species-wise variations in dispersal, and neutral, zero-sum dynamics (e.g., Lehman and Tilman, 2000; Chesson et al., 2001; Loreau et al., 2003; Loreau and de Mazancourt, 2008). At global scales, asynchronous dynamics can have stronger effects on stability than species richness (Valencia et al., 2020). Accordingly, evidence of the central role of species synchrony in the stability of ecosystem functions is supported by studies across ecosystems and taxonomic groups (reviewer by Gonzalez and Loreau, 2009).

The lack of correlation between functional and compositional measurements of stability could well reflect the effects of species compensation. For instance, high resistance and recovery (i.e., near-zero log-response ratios) in the functional domain were observed along ample ranges of both metrics in the compositional domain, and the negative function-composition correlations detected for resilience and recovery further reinforce the notion of functional compensation within these communities (Guelzow et al., 2017; Hillebrand et al., 2018; Hillebrand and Kunze, 2020). This indicates that a given level of functioning could be performed by communities with different compositions, which may have important consequences for ecosystem restauration (Figure 7A; Borja et al., 2010).

FIGURE 7

Figure 7. Summary of scenarios based on correlations between (A) functional and compositional stability domains and (B) stability dimensions. Panel (A) Systems exhibiting positive functional-compositional stability correlations may need a full recovery in species composition to achieve a full recovery in functioning (yellow quadrant); if not, no functional recovery would be possible (orange quadrant). Contrarily, negative functional-compositional correlations (e.g., barnacle communities in this study) imply that the altered community may be able to compensate the functioning of the pre-disturbance community, but also that restored communities are unable to recover pre-disturbance functional values. Panel (B) A positive correlation between two aspects of stability (e.g., functional resistance and invariability; this study) indicates that and conservation decisions could improve both aspects of stability; a negative correlation (e.g., functional resistance and resilience; this study) implies that policy trade-offs will be crucial to maintain either one or another aspect of stability. The lack of correlation among stability dimensions (e.g., all compositional dimensions in this study) indicates that, due to high uncertainty, multiple aspects are mandatory to characterize the stability of the ecosystem.

The Stability Dimensions Were Mostly Uncorrelated Across Community Types and Stability Domains

In the functional domain, the correlations among the different dimensions of stability ranged from negative to null and to positive. We were able to detect negative correlations between functional resistance and resilience and between resilience and invariability. Negative functional resistance-resilience and resilience-invariability correlations have been also observed in modeling and experiments conducted on experimental phytoplanktonic mesocosm communities (Hillebrand et al., 2018), and have been associated to comparatively high growth rates that allow faster recovery rates with stronger initial impacts (Harrison, 1979). Similarly to aquatic ciliate communities (Pennekamp et al., 2018), larger initial functional responses (i.e., lower resistance) were associated with more temporally stable functioning, which could respond to a positive relationship between species’ invariability and relative abundance (Grman et al., 2010; Lamy et al., 2020).

In the compositional domain, on the other hand, no statistically significant correlation between stability dimensions was detected. This high dimensionality implies that the four aspects of compositional stability analyzed in this study are fully orthogonal and are mandatory to completely describe the response of these communities to an ecologically realistic disturbance akin to our press removal of the dominant species. Uncorrelated stability dimensions have been highlighted by previous work when ecological communities face severe environmental disturbances (Allison, 2004; van Ruijven and Berendse, 2010; Donohue et al., 2013; Hillebrand et al., 2018). Contrarily, most empirical research has been based on one or two stability dimensions; the temporal coefficient of variation has been widely used to assess the stability in aggregate ecosystem functions in the absence of major disturbances (reviewed by Kéfi et al., 2019). Yet, restoration assessment after pulse and press disturbances would largely benefit from considering multiple components of stability, because, as we have shown, these components may well be decoupled during the process of recovery (Figure 7B; Donohue et al., 2016; Moreno-Mateos et al., 2020).

Perhaps, percentage cover would have been too coarse to detect differences stability among community types and positive correlations among the stability dimensions. Albeit we considered the multi-layered nature of macrobenthic communities in our in situ estimations, it would be the case that between-plot variations in the tri-dimensional projection of the communities would have generated unwanted variation in the stability responses. However, our results are well supported by previous work assessing other functional proxy, such as consumption rates and biomass (Ghedini et al., 2015; Hillebrand et al., 2018; Valencia et al., 2020). Percentage cover in sessile communities is the results of the complex interaction between preemptive competition, disturbance, and recolonization (Paine, 1980), which makes it a good proxy for functioning.

Conclusion and Future Directions

Our results allow us to conclude that ecological stability of the analyzed intertidal communities is a multidimensional and complex phenomenon. The experimental removal of locally dominant species elicited significant, and similar in magnitude, stability responses of communities in a large biogeographical province of the SE Pacific shore. These responses varied between functional and compositional domains, likely due to pervasive and strong asynchronous species dynamics.

Still is necessary to determine how different types and aspects of disturbances can influence stability dimensionality in these ecosystems. Albeit recent modeling work indicates that disturbance types such as random, species-specific, or local disturbances can have similar effects on stability dimensionality (Radchuk et al., 2019), contrasting the effects of pulse and press disturbances can enhance our mechanistic understanding of the ecological consequences of current climate crisis (Kéfi et al., 2019)—anthropogenic climate change involves the interdependent effects of both types of disturbances over multiple spatiotemporal scales (Harris et al., 2018). Regional manipulative experiments considering a more ample range of environmental filters and dispersal limitations would allow us model the role of spatial aspects of disturbance (Zelnik et al., 2018, 2019) in how stability dimensionality scales from sites to entire regions; to our best knowledge, the scale dependence of stability has been investigated only for single aspects of stability (Wang et al., 2017). In summary, realizing that stability is indeed a multidimensional concept will enhance our ability to predict the response of communities and metacommunities to future and ongoing anthropogenic disturbances.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

Author Contributions

NV and BB conceived the study. NV analyzed the data. BB and MA contributed with field work, reagents, and analytical material. NV wrote the manuscript. BB and MA contributed significantly to the writing of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was financially supported by the FONDECYT #1190529, FONDECYT #1141037, and FONDAP #15150003 (IDEAL) to NV. MA was supported by FONDECYT #1160223 and PAI-CONICYT #79150002-2016. BB acknowledges the support of FONDECYT #1181300 and the MUSEL Millennium Nucleus, funded by Iniciativa Científica Milenio.

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.

Acknowledgments

We thank all the enthusiastic undergraduate interns (Biología Marina; UACH), members of the Coastal Ecology Lab (Lafkenchelab, UACH), and the Changolab (CEAZA-UAI) who provided invaluable help during the fieldwork. Melissa González, José Pantoja, and Oscar Beytía were fundamental to complete the sampling schedule in the northern region. Constructive criticism by Fabio Bulleri, Matt Bracken and three reviewers greatly improved an early version of this manuscript.

Supplementary Material

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

References

Aguilera, M. A., Aburto, J. A., Bravo, L., Broitman, B. R., García, R. A., Gaymer, C. F., et al. (2019). “Chapter 29 - Chile: environmental status and future perspectives,” in World Seas: an Environmental Evaluation, Second Edn, ed. C. Sheppard (Cambridge, MA: Academic Press), 673–702. doi: 10.1016/b978-0-12-805068-2.00046-2

CrossRef Full Text | Google Scholar

Aguilera, M. A., and Navarrete, S. A. (2012). Functional identity and functional structure change through succession in a rocky intertidal marine herbivore assemblage. Ecology 93, 75–89. doi: 10.1890/11-0434.1

CrossRef Full Text | Google Scholar

Aguilera, M. A., Navarrete, S. A., and Broitman, B. R. (2013). Differential effects of grazer species on periphyton of a temperate rocky shore. Mar. Ecol. Prog. Ser. 484, 63–78. doi: 10.3354/meps10297

CrossRef Full Text | Google Scholar

Allison, G. (2004). The influence of species diversity and stress intensity on community resistance and resilience. Ecol. Monographs 74, 117–134. doi: 10.1890/02-0681

CrossRef Full Text | Google Scholar

Arnillas, C. A., and Cadotte, M. W. (2019). Experimental dominant plant removal results in contrasting assembly for dominant and non-dominant plants. Ecol. Lett. 22, 1233–1242.

Google Scholar

Atkinson, L. P., Valle-Levinson, A., Figueroa, D., De Pol-Holz, R., Gallardo, V. A., Schneider, W., et al. (2002). Oceanographic observations in chilean coastal waters between valdivia and concepcion. J. Geophys. Res. Oceans 107, 18–11.

Google Scholar

Baldanzi, S., Storch, D., Navarrete, S. A., Graeve, M., and Fernandez, M. (2018). Latitudinal variation in maternal investment traits of the kelp crab Taliepus dentatus along the coast of Chile. Mar. Biol. 165:12.

Google Scholar

Bartoń, K. (2020). MuMIn: Multi-Model Inference. (R package version 1.43.17. Available online at: https://CRAN.R-project.org/package=MuMIn).

Google Scholar

Bates, D., Machler, M., Bolker, B. M., and Walker, S. C. (2015). Fitting linear mixed-effects models using lme4. J. Stat. Software 67, 1–48.

Google Scholar

Benincà, E., Ballantine, B., Ellner, S. P., and Huisman, J. (2015). Species fluctuations sustained by a cyclic succession at the edge of chaos. Proc. Natl. Acad. Sci. U S A. 112, 6389–6394. doi: 10.1073/pnas.1421968112

PubMed Abstract | CrossRef Full Text | Google Scholar

Borja, Á, Dauer, D. M., Elliott, M., and Simenstad, C. A. (2010). Medium- and Long-term recovery of estuarine and coastal ecosystems: patterns, rates and restoration effectiveness. Estuaries Coasts 33, 1249–1260. doi: 10.1007/s12237-010-9347-5

CrossRef Full Text | Google Scholar

Broitman, B. R., Navarrete, S. A., Smith, F., and Gaines, S. D. (2001). Geographic variation of southeastern Pacific intertidal communities. Mar. Ecol. Prog. Ser. 224, 21–34. doi: 10.3354/meps224021

CrossRef Full Text | Google Scholar

Broitman, B. R., Véliz, F., Manzur, T., Wieters, E. A., Finke, G. R., Fornes, P. A., et al. (2011). Geographic variation in diversity of wave exposed rocky intertidal communities along central Chile. Revista Chilena de Historia Natural 84, 143–154. doi: 10.4067/s0716-078x2011000100011

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulleri, F., Benedetti-Cecchi, L., Cusson, M., Maggi, E., Arenas, F., Aspden, R., et al. (2012). Temporal stability of European rocky shore assemblages: variation across a latitudinal gradient and the role of habitat-formers. Oikos 121, 1801–1809. doi: 10.1111/j.1600-0706.2011.19967.x

CrossRef Full Text | Google Scholar

Camus, P. A. (2001). Marine biogeography of continental Chile. Revista Chilena de Historia Natural 74, 587–617.

Google Scholar

Caro, A. U., Navarrete, S. A., and Castilla, J. C. (2010). Ecological convergence in a rocky intertidal shore metacommunity despite high spatial variability in recruitment regimes. Proc. Natl. Acad. Sci. U S A. 107, 18528–18532. doi: 10.1073/pnas.1007077107

PubMed Abstract | CrossRef Full Text | Google Scholar

Chesson, P., Pacala, S., and Neuhauser, C. (2001). “Environmental niches and ecosystem functioning,” in The Functional Consequences of Biodiversity, eds A. Kinzig, S. Pacala, and D. Tilman (Princeton, NJ: Princeton University Press).

Google Scholar

Clarke, K. R. (1993). Non-parametric multivariate analyses of changes in community structure. Aus. J. Ecol. 18, 117–143. doi: 10.1111/j.1442-9993.1993.tb00438.x

CrossRef Full Text | Google Scholar

Connell, J. H., and Slatyer, R. O. (1977). Mechanisms of succession in natural communities and their role in community stability and organization. Am. Nat. 111, 1119–1144. doi: 10.1086/283241

CrossRef Full Text | Google Scholar

De Laender, F., Rohr, J. R., Ashauer, R., Baird, D. J., Berger, U., Eisenhauer, N., et al. (2016). Reintroducing environmental change drivers in biodiversity-ecosystem functioning research. Trends Ecol. Evol. 31, 905–915. doi: 10.1016/j.tree.2016.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Donohue, I., Hillebrand, H., Montoya, J. M., Petchey, O. L., Pimm, S. L., Fowler, M. S., et al. (2016). Navigating the complexity of ecological stability. Ecol. Lett. 19, 1172–1185.

Google Scholar

Donohue, I., Petchey, O. L., Montoya, J. M., Jackson, A. L., Mcnally, L., Viana, M., et al. (2013). On the dimensionality of ecological stability. Ecol. Lett. 16, 421–429.

Google Scholar

Dunstan, P. K., and Johnson, C. R. (2006). Linking richness, community variability, and invasion resistance with patch size. Ecology 87, 2842–2850. doi: 10.1890/0012-9658(2006)87[2842:lrcvai]2.0.co;2

CrossRef Full Text | Google Scholar

Economo, E. P. (2011). Biodiversity conservation in metacommunity networks: linking pattern and persistence. Am. Nat. 177, E167–E180.

Google Scholar

Garreaud, R. D., Rutllant, J. A., Muñoz, R. C., Rahn, D. A., Ramos, M., and Figueroa, D. (2011). VOCALS-CUpEx: the chilean upwelling experiment. Atmospheric Chem. Phys. 11, 2015–2029. doi: 10.5194/acp-11-2015-2011

CrossRef Full Text | Google Scholar

Ghedini, G., Russell, B. D., and Connell, S. D. (2015). Trophic compensation reinforces resistance: herbivory absorbs the increasing effects of multiple disturbances. Ecol. Lett. 18, 182–187. doi: 10.1111/ele.12405

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez, A., and Loreau, M. (2009). The causes and consequences of compensatory dynamics in ecological communities. Ann. Rev. Ecol. Evol. Systematics 40, 393–414. doi: 10.1146/annurev.ecolsys.39.110707.173349

CrossRef Full Text | Google Scholar

Gonzalez, A., Germain, R. M., Srivastava, D. S., Filotas, E., Dee, L. E., Gravel, D., et al. (2020). Scaling-up biodiversity-ecosystem functioning research. Ecol. Lett. 23, 757–776.

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

Grman, E., Lau, J. A., Schoolmaster, D. R., and Gross, K. L. (2010). Mechanisms contributing to stability in ecosystem function depend on the environmental context. Ecol. Lett. 13, 1400–1410. doi: 10.1111/j.1461-0248.2010.01533.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Guelzow, N., Muijsers, F., Ptacnik, R., and Hillebrand, H. (2017). Functional and structural stability are linked in phytoplankton metacommunities of different connectivity. Ecography 40, 719–732. doi: 10.1111/ecog.02458

CrossRef Full Text | Google Scholar

Harley, C. D. G. (2006). Effects of physical ecosystem engineering and herbivory on intertidal community structure. Mar. Ecol. Prog. Ser. 317, 29–39. doi: 10.3354/meps317029

CrossRef Full Text | Google Scholar

Harris, R. M. B., Beaumont, L. J., Vance, T. R., Tozer, C. R., Remenyi, T. A., Perkins-Kirkpatrick, S. E., et al. (2018). Biological responses to the press and pulse of climate trends and extreme events. Nat. Climate Change 8, 579–587.

Google Scholar

Harrison, G. W. (1979). Stability under environmental stress: resistance. resilience, persistence, and variability. Am. Nat. 113, 659–669. doi: 10.1086/283424

CrossRef Full Text | Google Scholar

Hillebrand, H., and Kunze, C. (2020). Meta-analysis on pulse disturbances reveals differences in functional and compositional recovery across ecosystems. Ecol. Lett. 23, 575–585. doi: 10.1111/ele.13457

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillebrand, H., Langenheder, S., Lebret, K., Lindström, E., Östman, Ö, Striebel, M., et al. (2018). Decomposing multiple dimensions of stability in global change experiments. Ecol. Lett. 21, 21–30. doi: 10.1111/ele.12867

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffmann, A. J., and Santelices, B. (1997). Flora Marina de Chile Central. Santiago: Ediciones Universidad Católica de Chile.

Google Scholar

Hormazabal, S., Shaffer, G., and Leth, O. (2004). Coastal transition zone off Chile. J. Geophys. Res. Oceans 109:C01021.

Google Scholar

Isbell, F., Gonzalez, A., Loreau, M., Cowles, J., Diaz, S., Hector, A., et al. (2017). Linking the influence and dependence of people on biodiversity across scales. Nature 546, 65–72. doi: 10.1038/nature22899

PubMed Abstract | CrossRef Full Text | Google Scholar

Ives, A. R., and Carpenter, S. R. (2007). Stability and diversity of ecosystems. Science 317, 58–62.

Google Scholar

Jara, H. F., and Moreno, C. A. (1984). Herbivory and structure in a midlittoral rocky community: a case in southern Chile. Ecology 65, 28–38. doi: 10.2307/1939455

CrossRef Full Text | Google Scholar

Kéfi, S., Dominguez-Garcia, V., Donohue, I., Fontaine, C., Thebault, E., and Dakos, V. (2019). Advancing our understanding of ecological stability. Ecol. Lett. 22, 1349–1356. doi: 10.1111/ele.13340

PubMed Abstract | CrossRef Full Text | Google Scholar

Lagos, N. A., Navarrete, S. A., Veliz, F., Masuero, A., and Castilla, J. C. (2005). Meso-scale spatial variation in settlement and recruitment of intertidal barnacles along the coast of central Chile. Mar. Ecol. Prog. Ser. 290, 165–178. doi: 10.3354/meps290165

CrossRef Full Text | Google Scholar

Lamy, T., Koenigs, C., Holbrook, S. J., Miller, R. J., Stier, A. C., and Reed, D. C. (2020). Foundation species promote community stability by increasing diversity in a giant kelp forest. Ecology 101:e02987.

Google Scholar

Lande, R., Engen, S., and Saether, B.-E. (2003). Stochastic Population Dynamics in Ecology and Conservation. New York, NY: Oxford University Press.

Google Scholar

Lara, C., Saldías, G., Cazelles, B., Rivadeneira, M., Haye, P., and Broitman, B. (2019). Coastal biophysical processes and the biogeography of rocky intertidal species along the south-eastern Pacific. J. Biogeography 46, 420–431. doi: 10.1111/jbi.13492

CrossRef Full Text | Google Scholar

Lehman, C. L., and Tilman, D. (2000). Biodiversity, stability, and productivity in competitive communities. Am. Nat. 156, 534–552. doi: 10.2307/3079056

CrossRef Full Text | Google Scholar

Letelier, J., Pizarro, O., and Nunez, S. (2009). Seasonal variability of coastal upwelling and the upwelling front off central Chile. J. Geophys. Res. Oceans 114, 1–16. doi: 10.1016/j.pocean.2011.07.019

CrossRef Full Text | Google Scholar

Loreau, M., and de Mazancourt, C. (2008). Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments. Am. Nat. 172, E48–E66.

Google Scholar

Loreau, M., Mouquet, N., and Gonzalez, A. (2003). Biodiversity as spatial insurance in heterogeneous landscapes. Proc. Natl. Acad. Sci. U S A. 100, 12765–12770. doi: 10.1073/pnas.2235465100

PubMed Abstract | CrossRef Full Text | Google Scholar

Lüdecke, D. (2020). ”sjPlot: Data Visualization for Statistics in Social Science”. (R package version 2.8.4: Available online at: https://CRAN. (R)-project.org/package=sjPlo).

Google Scholar

Micheli, F., Cottingham, K. L., Bascompte, J., Bjornstad, O. N., Eckert, G. L., Fischer, J. M., et al. (1999). The dual nature of community variability. Oikos 85, 161–169. doi: 10.2307/3546802

CrossRef Full Text | Google Scholar

Moreno, C. A. (2001). Community patterns generated by human harvesting on Chilean shores: a review. Aquatic Conserv. Mar. Freshwater Ecosystems 11, 19–30. doi: 10.1002/aqc.430

CrossRef Full Text | Google Scholar

Moreno, C. A., and Jaramillo, E. (1983). The role of grazers in the zonation of intertidal macroalgae of the Chilean coast. Oikos 41, 73–76. doi: 10.2307/3544348

CrossRef Full Text | Google Scholar

Moreno-Mateos, D., Alberdi, A., Morrien, E., Van Der Putten, W. H., Rodriguez-Una, A., and Montoya, D. (2020). The long-term restoration of ecosystem complexity. Nat. Ecol. Evol. 4, 676–685. doi: 10.1038/s41559-020-1154-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Mrowicki, R. J., Maggs, C. A., and O’connor, N. E. (2014). Does wave exposure determine the interactive effects of losing key grazers and ecosystem engineers? J. Exp. Mar. Biol. Ecol. 461, 416–424. doi: 10.1016/j.jembe.2014.09.007

CrossRef Full Text | Google Scholar

Nielsen, K. J., and Navarrete, S. A. (2004). Mesoscale regulation comes from the bottom-up: intertidal interactions between consumers and upwelling. Ecol. Lett. 7, 31–41. doi: 10.1046/j.1461-0248.2003.00542.x

CrossRef Full Text | Google Scholar

Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., Mcglinn, D., et al. (2019). ”vegan: Community Ecology Package”. R package version 2.5-6 ed.).

Google Scholar

Ospina-Alvarez, A., De Juan, S., Davis, K. J., González, C., Fernández, M., and Navarrete, S. A. (2020). Integration of biophysical connectivity in the spatial optimization of coastal ecosystem services. Sci. Total Environ. 733:139367. doi: 10.1016/j.scitotenv.2020.139367

PubMed Abstract | CrossRef Full Text | Google Scholar

Paine, R. T. (1980). Food webs: linkage, interaction strength and community infrastructure. J. Animal Ecol. 49, 667–685.

Google Scholar

Pennekamp, F., Pontarp, M., Tabi, A., Altermatt, F., Alther, R., Choffat, Y., et al. (2018). Biodiversity increases and decreases ecosystem stability. Nature 563, 109–112. doi: 10.1038/s41586-018-0627-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Pimm, S. L. (1984). The complexity and stability of ecosystems. Nature 307, 321–326.

Google Scholar

Pinochet, A., Garces-Vargas, J., Lara, C., and Olguin, F. (2019). Seasonal Variability of Upwelling off Central-Southern Chile. Remote Sensing 11:1737. doi: 10.3390/rs11151737

CrossRef Full Text | Google Scholar

R Core Team (2020). R: a Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

Google Scholar

Radchuk, V., De Laender, F., Cabral, J. S., Boulangeat, I., Crawford, M., Bohn, F., et al. (2019). The dimensionality of stability depends on disturbance type. Ecol. Lett. 22, 674–684. doi: 10.1111/ele.13226

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanford, E. (2002). The feeding, growth, and energetics of two rocky intertidal predators (Pisaster ochraceus and Nucella canaliculata) under water temperatures simulating episodic upwelling. J. Exp. Mar. Biol. Ecol. 273, 199–218. doi: 10.1016/s0022-0981(02)00164-8

CrossRef Full Text | Google Scholar

Scheffer, M., Bascompte, J., Brock, W. A., Brovkin, V., Carpenter, S. R., Dakos, V., et al. (2009). Early-warning signals for critical transitions. Nature 461, 53–59.

Google Scholar

Stachowicz, J. J., Graham, M., Bracken, M. E. S., and Szoboszlai, A. I. (2008). Diversity enhances cover and stability of seaweed assemblages: the role of heterogeneity and time. Ecology 89, 3008–3019. doi: 10.1890/07-1873.1

CrossRef Full Text | Google Scholar

Svensson, C. J., Jenkins, S. R., Hawkins, S. J., and Åberg, P. J. O. (2005). Population resistance to climate change: modelling the effects of low recruitment in open populations. Oecologia 142, 117–126. doi: 10.1007/s00442-004-1703-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Tapia, F. J., Largier, J. L., Castillo, M., Wieters, E. A., and Navarrete, S. A. (2014). Latitudinal discontinuity in thermal conditions along the nearshore of central-northern Chile. PLoS One 9:e110841. doi: 10.1371/journal.pone.0110841

PubMed Abstract | CrossRef Full Text | Google Scholar

Tejada-Martinez, D., López, D. N., Bonta, C. C., Sepúlveda, R. D., and Valdivia, N. (2016). Positive and negative effects of mesograzers on early-colonizing species in an intertidal rocky-shore community. Ecol. Evol. 6, 2045–7758.

Google Scholar

Thiel, M., Macaya, E. C., Acuna, E., Arntz, W. E., Bastias, H., Brokordt, K., et al. (2007). The Humboldt Current System of northern and central Chile. Oceanography Mar. Biol. Ann. Rev. 45, 195–344.

Google Scholar

Tilman, D. (1996). Biodiversity: population versus ecosystem stability. Ecology 77, 350–363. doi: 10.2307/2265614

CrossRef Full Text | Google Scholar

Valdivia, N., Aguilera, M. A., Navarrete, S. A., and Broitman, B. R. (2015). Disentangling the effects of propagule supply and environmental filtering on the spatial structure of a rocky shore metacommunity. Mar. Ecol. Prog. Ser. 538, 67–79. doi: 10.3354/meps11493

CrossRef Full Text | Google Scholar

Valdivia, N., Gonzalez, A. E., Manzur, T., and Broitman, B. R. (2013). Mesoscale variation of mechanisms contributing to stability in rocky shore communities. PLoS One 8:e54159. doi: 10.1371/journal.pone.0054159

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdivia, N., Segovia-Rivera, V., Fica, E., Bonta, C. C., Aguilera, M. A., and Broitman, B. R. (2017). Context-dependent functional dispersion across similar ranges of trait space covered by intertidal rocky shore communities. Ecol. Evol. 7, 1882–1891. doi: 10.1002/ece3.2762

PubMed Abstract | CrossRef Full Text | Google Scholar

Valencia, E., De Bello, F., Galland, T., Adler, P. B., Lepš, J., E-Vojtkó, A., et al. (2020). Synchrony matters more than species richness in plant community stability at a global scale. Proc. Natl. Acad. Sci. 117, 24345–24351.

Google Scholar

Valone, T. J., and Balaban-Feld, J. (2018). Impact of exotic invasion on the temporal stability of natural annual plant communities. Oikos 127, 56–62. doi: 10.1111/oik.04591

CrossRef Full Text | Google Scholar

Van Holt, T., Moreno, C. A., Binford, M. W., Portier, K. M., Mulsow, S., and Frazer, T. K. (2012). Influence of landscape change on nearshore fisheries in southern Chile. Global Change Biol. 18, 2147–2160. doi: 10.1111/j.1365-2486.2012.02674.x

CrossRef Full Text | Google Scholar

van Ruijven, J., and Berendse, F. (2010). Diversity enhances community recovery, but not resistance, after drought. J. Ecol. 98, 81–86. doi: 10.1111/j.1365-2745.2009.01603.x

CrossRef Full Text | Google Scholar

Wang, S. P., and Loreau, M. (2014). Ecosystem stability in space: alpha, beta and gamma variability. Ecol. Lett. 17, 891–901.

Google Scholar

Wang, S. P., and Loreau, M. (2016). Biodiversity and ecosystem stability across scales in metacommunities. Ecol. Lett. 19, 510–518. doi: 10.1111/ele.12582

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Loreau, M., Arnoldi, J.-F., Fang, J., Rahman, K. A., Tao, S., et al. (2017). An invariability-area relationship sheds new light on the spatial scaling of ecological stability. Nat. Commun. 8:15211.

Google Scholar

Whalen, M. A., Aquilino, K. M., and Stachowicz, J. J. (2016). Grazer diversity interacts with biogenic habitat heterogeneity to accelerate intertidal algal succession. Ecology 97, 2136–2146. doi: 10.1890/15-1633.1

CrossRef Full Text | Google Scholar

Wickham, H., Averick, M., Bryan, J., Chang, W., Mcgowan, L. D. A., François, R., et al. (2019). Welcome to the tidyverse. J. Open Source Software 4:1686. doi: 10.21105/joss.01686

CrossRef Full Text | Google Scholar

Wootton, J. T. (2010). Experimental species removal alters ecological dynamics in a natural ecosystem. Ecology 91, 42–48. doi: 10.1890/08-1868.1

CrossRef Full Text | Google Scholar

Yachi, S., and Loreau, M. (1999). Biodiversity and ecosystem productivity in a fluctuating environment: the insurance hypothesis. Proc. Natl. Acad. Sci. U S A. 96, 1463–1468. doi: 10.1073/pnas.96.4.1463

PubMed Abstract | CrossRef Full Text | Google Scholar

Zelnik, Y. R., Arnoldi, J.-F., and Loreau, M. (2018). The impact of spatial and temporal dimensions of disturbances on ecosystem stability. Front. Ecol. Evol. 6:224. doi: 10.3389/fevo.2018.00224

PubMed Abstract | CrossRef Full Text | Google Scholar

Zelnik, Y. R., Arnoldi, J.-F., and Loreau, M. (2019). The three regimes of spatial recovery. Ecology 100:e02586. doi: 10.1002/ecy.2586

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y. H., He, N. P., Loreau, M., Pan, Q. M., and Han, X. G. (2018). Scale dependence of the diversity-stability relationship in a temperate grassland. J. Ecol. 106, 1277–1285. doi: 10.1111/1365-2745.12903

PubMed Abstract | CrossRef Full Text | Google Scholar