The experimental field study was carried out in two sites separated by ca. 22 km on the wave-exposed rocky shores of northern-central Chile (“Limarí” and “Punta de Talca”, ∼30°43′ S and 30°55′ S, respectively). In both sites, we assessed the effect of the experimental extinction of the dominant canopy-forming alga Mazzaella on community invariability. This red corticated alga inhabits the mid and low intertidal rocky shore between Coquimbo and Punta Arenas (Hoffmann and Santelices, 1997; Montecinos et al., 2012). Mazzaella can be highly dominant, with canopies covering up to ∼65% of secondary space (Santelices, 1990; Hoffmann and Santelices, 1997; Valdivia et al., 2021). The thallus can reach near to 30 cm in maximum length and the thick fronds can ameliorate environmental conditions beneath them (Hoffmann and Santelices, 1997). Particularly in our study sites, common sizes range from 5 to 15 cm in length (Supplementary Figure 1).

Direct harvesting of Mazzella constitutes a strong anthropogenic disturbance that drives abundance change in intertidal communities of the southeast Pacific coast (Buschmann et al., 2001). This activity occurs for the extraction of carrageenans used in food and pharmaceutical industry (Buschmann et al., 2001). Landings data from 2002 to 2017 shows that near to 3 tons of dry alga are extracted annually in Chile (Lopez, 2019).

Underneath Mazzaella canopies, a diverse assemblage of algae is observed, including the red algae Nothogenia fastigiata, Corallina officinalis var. chilensis, Gelidium sp. and Pyropia sp., and the green Chaetomorpha firma and Ulva spp. Also, the canopy proves habitat for macro and mesograzers such as the chitons Chiton granosus and C. barnesii, keyhole limpets (Fissurella spp.), the gastropod Siphonaria lessoni, the littorinid Austrolittorina araucana, and at least six species of the genus Scurria (Moreno and Jaramillo, 1983; Jara and Moreno, 1984; Moreno et al., 1984; Rivadeneira et al., 2002).

The experiment was established and monitored from January 2015 to August 2017, covering a sampling period of 32 months and at least three full seasonal cycles. The observations were conducted on 20 experimental plots equally distributed across both sites. Both sites were open access areas, where local fishermen and visitors are allowed to harvest marine resources directly (Buschmann et al., 2001). Also, both sites are exposed to similar environmental and oceanographic conditions (Aguilera et al., 2013; Valdivia et al., 2013, 2015).

In each site, we deployed at the mid intertidal zone ten 30 cm × 30 cm experimental plots, which were permanently marked with stainless-steel bolts. The size of the plots was chosen to minimize the impact on the local intertidal community. Experimental plots were placed within areas of high secondary cover (ca. 20–65% percent cover) of the locally dominant alga Mazzaella, and on flat and horizontal rocky surfaces, avoiding crevices and tide pools to reduce variability associated to substrate heterogeneity. Half of the experimental plots were haphazardly selected for the experimental treatment in which Mazzaella was completely removed with scrappers and chisels, leaving all the understory species intact (i.e., disturbed plots). The remaining half of the plots were kept unmanipulated (i.e., controls plots). Therefore, considering both sites together, our experimental design included a total of 10 control plots and 10 dominant species-removed plots. All settlers of Mazzaella were removed every three months in the disturbed plots. Three months between manipulations has been shown earlier to be an appropriate timespan to functionally exclude dominant species from intertidal communities in these sites and elsewhere (Bulleri et al., 2012; Aguilera et al., 2013; Valdivia et al., 2015).

All plots were sampled immediately before the removal of dominant canopies of Mazzaella, then two months after and then every three months until the end of the experiment in August of 2017. During each monitoring we estimated species abundance using a 30 cm × 30 cm frame divided in 25 fields with a monofilament line, following protocols described elsewhere (e.g., Broitman et al., 2011; Valdivia et al., 2020). All algae and invertebrates (c.a. >5 mm) occurring on each plot were identified in situ and the organisms were classified at the species level whenever possible. Abundance of algal species was estimated as substrate percentage cover (1% resolution), while grazers were registered as density (number of individuals per m²). All manipulations and observations were conducted during diurnal low-tide hours (ca. 1.5 m tidal range).

Abundance data of algae and grazers (percent cover and density, respectively) were transformed to biomass (dry weight g m^–2). This allowed us to analyze together mobile and sessile species. For algae, biomass was estimated through calibration curves generated from destructive samples obtained in the field (10–20 replicates for each of the 19 species identified in experimental plots). For algae whose abundances in the field were too low to generate a calibration curve, we used published linear and nonlinear models of percentage cover-dry weight for morphologically similar taxa (Ko et al., 2008). In addition, the biomass of grazer species was estimated as the product between each species’ density (ind. m^–2) and the average individual dry weight (g) obtained from Camus et al. (2013) and in situ samples (5–7 replicates for each of the main group of grazers: Littorinid snails, Scurria spp., Chiton spp., and Fissurella spp.). For biomass estimation, we excluded crustose algae due to logistic constraints, sessile invertebrates because their cover was always less than 1%, and predators because their frequency of occurrence and densities were very low (e.g., starfish and crabs). Community biomass was calculated for each plot as the sum of all individual species biomasses.

Community invariability was estimated for each plot as the inverse coefficient of variation (μ/σ) in community biomass over the 32 months of monitoring. This included together the biomass of the dominant canopy, understory algae, and grazers. This measure of invariability standardizes the temporal mean of community biomass to its temporal standard deviation; therefore, higher values indicate that community biomass is more constant over time (Tilman, 1999).

Population invariability was estimated following Thibaut and Connolly (2013) as the inverse of the weighted mean coefficient of variation in the biomass of all populations (species) within a sampling plot:

where μ_i and σ_i denote the temporal mean and standard deviation of biomass of each population i, N is the number of populations, and μ_c is the temporal mean of community biomass. Larger values of population invariability indicate that, on average, population abundances are more constant over time.

Following Sasaki et al. (2019), population asynchrony (φ_P) was estimated as:

where φ_P is the variance of community biomass (SxT2) divided by the squared sum of the standard deviations of all N individual populations (S_xi). Values around 0 indicate that biomass fluctuations among populations are perfectly synchronized, while values close to 1 will indicate that are perfectly asynchronized.

Invariability of dominants was estimated for control plots as the inverse coefficient of variation (μ/σ) of Mazzaella’s biomass over time. For disturbed plots we estimated the invariability of the “alternative dominants” as the invariability of the summed biomass of the foliose algae Ulva spp. and Pyropia sp., which were the most abundant algae when Mazzaella was absent. The separate estimation of dominant’s invariability in control and disturbed communities was carried out to assess, in separate statistical models, the contribution of dominant and alternative dominants to the invariability of the whole community, as well in the association among drivers that support invariability. Finally, mean species richness was estimated as the temporal mean of the number of taxonomic identities for each plot.

Additionally, local environmental variability during the experiment was assessed by analyzing sea surface temperature (publicly available at http://www.ceazamet.cl/) using OnSet HOBO^® water temperature dataloggers, which were deployed in both study sites, Limarí and Punta de Talca, at ca. 1m below low tide sea level following housed into concrete blocked anchored to the rocks with chain links. Deployment methods are extensively described elsewhere (Valdivia et al., 2013, 2015). Daily minimum, maximum, and average temperature data were derived from the temperature time series recorded during the experiment. Seawater temperature is a useful proxy for local community-wide environmental conditions as it correlates with other environmental variables such as nutrient and Chlorophyll-a concentration (Nielsen and Navarrete, 2004; Witman et al., 2008), and also through its influence on metabolic rates for ectotherms, animal and plant phenology, and population growth rates of coastal marine organisms (Strathmann et al., 2002; Baldanzi et al., 2018; Suárez et al., 2020).

Before testing our working hypotheses, we verified whether it was possible to pool the data from both sites. This was done to increase replication and thus statistical power for hypothesis testing. We used one-way analysis of variance (ANOVA) to compare mean percentage cover of Mazzaella, mean species richness, population asynchrony, and community invariability between both sites. For cover of Mazzaella, we analyzed only the control plots at each site (10 plots in total). For the other variables, we used the entire dataset (20 plots in total). Since we did not detect statistical differences between both sites in any variable (see Supplementary Table 1), we decided to analyze and report the result of both sites together.

We used information theory to assess the strength of empirical evidence supporting each of our hypotheses (Burnham and Anderson, 2002). We developed a set of general linear models (LM) linking community invariability to different combinations of predictor variables, in accordance with each hypothesis (see models in Table 2). These models were constructed on the basis of empirical and theoretical evidence supporting the premises described in the introduction section (see hypotheses H1, H2, and H3 in the introduction and Table 1). Model parameters were estimated through maximum likelihood. All numerical variables were centered and standardized before fitting.

Models were ranked according to bias-corrected Akaike Criterion of Information (AICc; Burnham et al., 2011). Model selection was based on Δ AICc, expressed a AICc_i − AICc_min, the Akaike weight, which represents the probability of each model as wi=Prob{gi|data}=li∑j=1Rlj, where g_i is model i and l_i = L(g_i|data = e^−0.5Δ_i, and the evidence ratio of the top model, calculated as w_top/w_j (Burnham et al., 2011). Among the ranked models, we selected those with a Δ AICc < 6, as our “top model set” (Richards, 2005, 2008; Harrison et al., 2018).

Since AICc tends to favor too complex models (i.e., models including many variables; Burnham and Anderson, 2002), we used the “nesting rule” to reduce the number of models used for inference from the top model set (Harrison et al., 2018). The nesting rule consists of selecting the model with the strongest empirical support (lowest AICc) and that is nested within more complex models with delta AICc < 6 (Richards, 2005, 2008; Harrison et al., 2018).

In addition, we used LMs to assess the effect of seasons and treatments on community biomass and that of each functional group. Model diagnostics were checked by visual inspection of quantile-quantile plots and fitted vs. adjusted residual plots. All the graphs and analyses were developed and fitted in R programming environment 4.0.2, using the sjPlot, ggplot2 packages (R Core Team, 2020). Community stability and synchrony were computed with “codyn” R package (Hallett et al., 2016). Linear models were computed by the “stats” package (R Core Team, 2020) and the “MuMIn” package was used to estimate AICc and Akaike weight scores (Barton, 2020).

Within the top model set (delta AICc < 6; Table 3), we selected model 12 because it had the strongest empirical support and it was a more parsimonious version of the others (i.e., nesting rule). This model included the separate effects of Mazzaella removal and population asynchrony on community invariability. The empirical support (evidence ratio) for model 12 was 4.6 times that of model 13, and 17.6 times that of model 15 (both in the top model set). Model 12 accounted for 75 % of community invariability (R² = 0.75, Table 4). The removal of Mazzaella decreased community invariability (Figure 1; difference between means = 1.16 [SE = 0.24], t-value = −4.9). Independent of the experimental removal, community invariability increased lineally with increasing population asynchrony (Figure 1; slope = 0.67 [0.12], t-value = 5.6).

In control plots, mean biomass of dominant canopies ranged between 30 and 300 g m^–2 during the experiment (see purple symbols in Figure 2A). Despite the large range of variation, canopies dominated in biomass most of the time. An initial decrease in canopy biomass in controls was followed by a peak increase in biomass of foliose algae; however, foliose algae biomass remained near to zero afterwards (cyan symbols in Figure 2A). The biomass of filamentous algae ranged from 0 to ∼70 g m^–2 and remained as subdominants during all the experiment, except when canopies decreased almost to zero in the early autumn of 2016 (blue symbols in Figure 2A). Branched algae ranged from 0 to ∼100 g m^–2 and fluctuated almost synchronically with dominant canopies during 2017 (yellow symbols in Figure 2A; see species composition of each group in Supplementary Table 2).

When canopies of Mazzaella were experimentally removed, the community became dominated by foliose, branched (corticated), and filamentous algae. We observed a large initial increase in foliose algae, which fluctuated from 20 to 250 g m^–2 in biomass and were the group with the highest contribution to community biomass until mid 2016 for this treatment (Figure 2B). Filamentous and branched algae were subdominant in these plots and fluctuated from 10 to 70 g m^–2 in biomass, showing a low contribution to overall community biomass. After a decrease in mid 2016, filamentous and branched algae remained in the same range of biomass and did not compensate the decrease in foliose algae (Figure 2B).

In controls, the biomass of macrograzers ranged from 5 to 55 g m^–2, peaking in early autumn of 2015 and early winter of 2016 (light-blue symbols in Figure 2C). Mesograzers, on the other hand, were less variable and ranged from 15 to 25 g m^–2 (pink symbols in Figure 2C). When Mazzaella was experimentally removed, macrograzers almost disappeared (0–5 g m^–2). An opposite response was observed for mesograzers, such as Siphonaria lessoni and Scurria spp., which increased in density and biomass, ranging from 10 to 40 g m^–2 (Figure 2D) and sometimes even overshooting the controls (from mid 2015 to late 2016).

Mean sea surface temperature ranged from 12.5 to 16.5° C during the experiment (Figure 3). The lowest temperature was reached during austral Autumn months, while the maximum temperature was reached in Summer months, particularly in January 2016. Community biomass showed no differences between treatments in the cold seasons (Figure 4). Between-treatment differences were evident in summer, where control plots reached almost twice the biomass of disturbed plots (model R² = 0.203; Figure 4 and Supplementary Table 3). The difference in community biomass between treatments observed in Summer was highly influenced by the increase in canopy biomass in controls (R² = 0.490; Supplementary Figure 2 and Supplementary Table 4). Seasonal effects were also observed for filamentous algae, in which biomass increased in Winter regardless of the experimental treatment (R² = 0.104; Supplementary Figure 2 and Supplementary Table 4). On the other hand, foliose algae showed the highest biomass in Autumn. Also, differences between treatments were observed in the biomass of foliose algae, with disturbed plots reaching almost two to three times the biomass of foliose algae in controls (R² = 0.206; Supplementary Figure 2 and Supplementary Table 4). The biomass of macrograzers was not affected by seasons but only by the experimental treatment (R² = 0.155), while the biomass of mesograzers tended to increase in Spring in the removal treatment (Supplementary Figure 2 and Supplementary Table 5).

Following the calculation of invariability (mean divided by the standard deviation; Tilman, 1999), the removal of canopies reduced community invariability respect to controls by decreasing the mean and/or increasing the standard deviation of community biomass (Bulleri et al., 2012). The removal of canopies led to a community dominated by low-biomass opportunistic foliose algae and mesograzers (Supplementary Table 2). Indeed, the removal of canopies decreased almost to zero the biomass of macrograzers, which are active consumers of algae such as the encrusting Hildenbrandia sp., calcareous branched Corallinaceae, and foliose Ulva spp. (Camus et al., 2013). Thus, the absence of these macrograzers in the removal treatment could explain the increase of foliose, branched, and filamentous species.

The species that became dominant under the canopy-removal condition were algae characterized by high colonization and mortality rates, low resistance to grazers, and physiological adaptations to environmental stress (Littler and Littler, 1984; Santelices, 1990; Duffy and Hay, 1994). In agreement with their opportunistic lifestyle, these algae also display large surface areas (i.e., foliose algae), which allows high productivity and growth rates (Littler and Littler, 1984; Steneck and Dethier, 1994). Consequently, the demographic stochasticity of these species may increase biomass variability at the community level (i.e., reduce community invariability; Stachowicz et al., 2008), as predicted by theoretical models of invariability (Loreau and de Mazancourt, 2013). Particularly, species such as Ulva spp. and Pyropia orbicularis can show high growth rates under broad environmental fluctuations and abiotic stress, quickly colonizing available niches as opportunistic species (Rosenberg and Ramus, 1982; Pérez-Mayorga et al., 2011). Also, elevated temperatures can trigger reproduction and settlement of some Ulva spp. (Gao et al., 2017). These foliose algae usually characterize early successional stages of intertidal communities, before the establishment of larger competitors like canopies or invertebrates (Foster et al., 2003; Edwards and Connell, 2012). Probably the absence of Mazzaella allowed constant available space for recruitment of these foliose algae, but seasonally variable disturbances—such as Winter extreme waves or enhanced mesograzer activity in Spring—might have limited their biomass. In this way, the demographic dynamics of the newly dominant species may have increased the temporal variance in total community biomass, decreasing thus community invariability.

The positive asynchrony-invariability relationship has been previously predicted by theory and has been found in a wide range of ecosystem across marine and terrestrial habitats (de Mazancourt et al., 2013; Valencia et al., 2020; Kigel et al., 2021). The positive effect of population asynchrony on community invariability has been associated to competitive interactions, differential species’ environmental tolerances, and trophic dynamics (Doak et al., 1998; Klug et al., 2000; Lehman and Tilman, 2000; Loreau and de Mazancourt, 2013). For example, the feeding activity of grazers can limit the expansion of subdominant algae after experimental disturbances such as canopy removal, increased nutrient load, and ocean acidification (Connell and Russell, 2010; Ghedini et al., 2015). Similarly, we observed an increase in the abundance of mesograzers in the canopy-removal plots (Figure 2C), as a result of milder interspecific competition with macrograzers and the increased habitat complexity provided by branched algae (e.g., Nothogenia fastigiata, Gelidium sp., Ahnfeltiopsis sp. and Trematocarpus dichotomus; O’Brien and Scheibling, 2018; Bertolini, 2018). Mesograzers, in turn, could have “dampened” the overshoot of foliose algae, which could explain the pattern of temporal variation observed for these algae (Figure 2B). Indeed, intertidal mesograzers can exert an important control on the abundance and diversity of algae and even sessile invertebrates (Aguilera et al., 2015; Tejada-Martínez et al., 2016). Specifically, mesograzers like Siphonaria lessoni and Scurria spp. are generalist consumers that consume algae such as Ulva spp., Ulvella sp., Hildenbrandia sp., and Codium dimorphum (Camus et al., 2013). In general, the feeding activity of (macro and meso) grazers on propagules and macroscopic stages can alter the dynamics of algal populations (Iken, 2012). These effects can be particularly strong during early successional stages, and thus, on the composition of algal communities (Arrontes, 1999; Aguilera and Navarrete, 2007; Aguilera, 2011; Iken, 2012). Additionally, moderate levels of herbivory by mesograzers could be beneficial during succession of intertidal communities dominated by Mazzaella (Aguilera and Navarrete, 2007). Therefore, a proportional increase in grazing, a trophic compensation, could well have counteracted the expansion of opportunistic algae as a response to canopy removal.

Theoretical, laboratory and field-based studies demonstrate that species richness generates invariability (Romanuk et al., 2006; Tilman et al., 2006; Loreau and de Mazancourt, 2013; Hallett et al., 2014). In our study, however, the empirical evidence supporting the independent effects of canopy removal and asynchrony on community invariability were ca. 18 times stronger than that supporting the effects of species richness. In the same line, a recent study of plant communities at the global scale has indicated that synchrony could be more relevant than species richness per se (Valencia et al., 2020). Asynchrony is indeed one of the three main mechanisms underpinning the stabilizing effects of species richness (de Mazancourt et al., 2013). Our results suggest therefore that the other two mechanisms driving the diversity-stability relationship—overyielding and reduced observation error—may have been less relevant than asynchrony in the analyzed community.

In addition to dominants, rare species can have strong effects on community invariability. For example, recent theoretical work suggests that perturbations that target rare species can have disproportionate effects on community dynamics, due to indirect effects cascading through the assemblage (Säterberg et al., 2019). In addition, the abundances of rare species are usually highly variable over time, which can have a strong influence on invariability averaged across all populations and even lead to negative diversity-stability relationships at the community level (Arnoldi et al., 2019). In our study, crustose algae and low-biomass filamentous and branched algae represented rare species. Future research could focus on the potential effect of these species on population and community stability. Our results support the idea that multiple factors, in addition to species richness, should be investigated to understand community invariability and other aspects of stability (Ives and Carpenter, 2007).

Despite canopy-removal plots exhibited lower invariability than control plots, this community-level property was independent from the invariability generated by the canopy of Mazzaella. Mazzaella exerts indirect effects on community invariability via trophic links and competitive interactions with the understory species. Alternatively, it can be the case the contribution of Mazzaella to mean community biomass was much larger than the “dampening” effect of this species on biomass standard deviation. Dominant species have been shown to influence both, the mean and standard deviation of total community biomass (Grman et al., 2010; Valdivia et al., 2013).