Impact of multiple disturbances and stress on the temporal trajectories and resilience of benthic intertidal communities

Coastal ecosystems face severe environmental change and anthropogenic pressures that affect both the structure and functioning of communities. Understanding the response and resilience of communities that face multiple simultaneous disturbances and stresses becomes essential. We observed the recovery of a rocky intertidal subarctic macrobenthic community dominated by a macroalgal canopy (Fucus spp.), a habitat-forming species, over a period of 14 months. Using 0.25-m plots, we ran an in situ onepulse experiment (removal of all materials to bare rock and then burning of the surface) followed by a full orthogonal factorial design of three press-type disturbances or stresses: grazer reduction, canopy removal, and nutrient enrichment. We evaluated the single and interactive effects of the three disturbances and stresses on species diversity and abundance structure. Of all the main effects, canopy removal has the most severe impact, resulting in decreased biomass, richness, and diversity, as well as an altered community structure. Canopy-removed plots had fewer invertebrates and more ephemeral algae; beyond this, however, there was minimal effect from grazer reduction and nutrient enrichment acting individually. We categorized the interaction types of all significant interaction effects: Canopy removal had a dominant effect over grazer reduction on richness, and it also dominated over nutrient enrichment on diversity and evenness. Nutrient enrichment and canopy removal had a negative synergistic interaction effect on richness at the end of the experiment. Without stressors, 11 months were required to achieve full recovery. The three stressors affected recovery time differently, depending on the identity and the number of stressors. Three stressors generally increased the time of recovery or even prevented recovery from being fully attained. Moreover, community structure and composition of plots subjected to the triple-stressor treatment had not fully recovered by the end of the study. Our results suggest that multiple stressors may interact on community indices and structure and that their interaction cannot be predicted from the outcome of single stressor studies. The inclusion of multiple disturbances and stresses in field experiments provides a better understanding of the mechanisms that shape community structure and their functioning following various forms


INTRODUCTION
Coastal ecosystems are subjected to both natural and anthropogenic disturbances and stresses, for example, eutrophication, sedimentation, and habitat loss Beck 2007, Halpern et al. 2008). These various changes can have drastic negative effects on ecosystems, such as biodiversity loss, that impact community functioning and stability (Stachowicz et al. 2002, Lotze et al. 2006, Hooper et al. 2012. Climate change and biodiversity loss are driving ecologists to study the relationship between biodiversity and ecosystem functioning , Gamfeldt et al. 2015. Biodiversity plays an important role in ecosystem stability (Loreau et al. 2001, Hillebrand et al. 2008, ecosystem functioning (Yachi andLoreau 1999, Hooper et al. 2005), offering resistance and resilience to disturbance and stress (Allison 2004, Aquilino andStachowicz 2012). Resilience (or engineering resilience) is an essential ecosystem function that helps maintain biodiversity and functioning. It is defined as the ability of a system to mitigate an impact and recover following a disturbance. Recovery is the ability of a system to return to a reference state following a disturbance (Ingrisch and Bahn 2018). Recovery can be measured in various ways. This includes temporal trajectories of different community characteristics, such as its abundance structure and recovery rate (Lotze et al. 2011, Duarte et al. 2015). All measured variables should be compared to an undisturbed reference (Underwood 1989, Ingrisch andBahn 2018) to account for spatial and temporal variation. Recovery can be influenced by the intensity and timing of a disturbance, the presence of intact stands of habitat-forming species, and the availability of propagules (Oliveira et al. 2011, O'Leary et al. 2017). On rocky shores, the rapid recovery of communities (within 1-5 yr) following a disturbance may be due to the open nature of this system that provides high numbers of propagules and larvae from unaffected shores , Crowe et al. 2000.
The presence of habitat-forming species can increase community stability (Maggi et al. 2009, Bulleri et al. 2012). Habitat-forming species (or foundation species; ecosystem engineers sensu Jones et al. 1994) play an essential role in ecosystem functioning by sheltering and protecting numerous organisms (Bertness et al. 1999). The loss of habitat-forming species negatively affects the surrounding community by reducing species richness and abundance as well as through changes in community structure and composition (Herkul and Kotta 2009, Rueda et al. 2009, Watt and Scrosati 2013. Intertidal zones of rocky shore habitats are often dominated by habitat-forming macroalgae that are considered as key species in these systems (Hawkins and Hartnoll 1983). These macroalgae provide complex habitats, dampen stresses-such as desiccation and extreme temperatures-limit water movement, provide shelter, and control biomass production and species diversity (Bertness et al. 1999, Jenkins et al. 1999, Benedetti-Cecchi et al. 2001, Heck et al. 2003, Eriksson et al. 2006. Presently, habitat-forming macroalgae are declining worldwide as a consequence of climate change and more regional and local-scale human impacts Beck 2007, Hawkins et al. 2009). Their loss can increase the recruitment of ephemeral algae (Schiel and Lilley 2007). Consequently, the disappearance of habitat-forming macroalgae has a deleterious impact on the associated assemblages through a reduction in community richness and abundance (Watt and Scrosati 2013) and an altering of composition structures (Lemieux and Cusson 2014).
Similarly, grazers play an important role in structuring communities within rocky intertidal shore habitats Underwood 1997, Atalah andCrowe 2010). The influence of grazers is particularly important during the first stages of algal colonization (Korpinen et al. 2007, Guerry 2008, and this influence varies depending on the environment, nutrient levels, and season (Lotze et al. 2001, Freidenburg et al. 2007. Grazer presence may accelerate the succession of macroalgae through the preferential consumption of ephemeral algae and thus affect algal density and diversity (Lubchenco 1983, Aquilino and Stachowicz 2012, Alestra and Schiel 2014. Coastal ecosystems experience constant stress due to human impacts (Halpern et al. 2015, Cloern et al. 2016). Among the many forms of stress, eutrophication exerts a marked pressure on coastal communities. Nitrogen fixation by humans has increased the usable forms of nitrogen for organisms by almost twofold over the last century (Vitousek et al. 1997, Fowler et al. 2013. Eutrophication caused by an excess of nitrogen impacts the benthic community structure by changing species composition and reducing macroalgal diversity (Worm andLotze 2006, Arevalo et al. 2007). Depending on the nutrient limitations at a site, higher nutrient loading may promote ephemeral algae settlement and growth (Bracken and Nielsen 2004, Kraufvelin 2007, Bertocci et al. 2017. Community resistance to eutrophication may be enhanced by habitatforming macroalgae (Eriksson et al. 2006). Grazers may reduce the effects of nutrient enrichment by preferential consumption of ephemeral algae (Atalah and Crowe 2010).
Multiple interacting anthropogenic and natural disturbances and stresses often co-occur in ecological communities. Their cumulative effects are often considered as being additive, although they may accumulate in a multiplicative manner, or they may not even accumulate, or they may show dominance from one stressor (Halpern et al. 2007, Côt e et al. 2016. Interactions may be synergistic or antagonistic, and they are unpredictable (Darling and Côt e 2008, Lyons et al. 2015, Côt e et al. 2016. Moreover, the frequency of synergistic interactions in the marine environment may be higher when communities are exposed to three stressors instead of two (Crain et al. 2008). However, the mention of synergistic interactions in the literature is growing much faster than mention of antagonistic or additive effects, and these claims are often being made without proper testing (Côt e et al. 2016). It is rare to find studies that include three or more disturbances or stresses in intertidal rocky habitats dominated by macroalgae (Strain et al. 2014); therefore, more of such studies are necessary. Single effects of canopy or grazer removal and enrichment are well known, yet few studies have looked at the interaction of biodiversity loss and stressors. Multiple interactive effects on both macroalgal and macroinvertebrate assemblage structures remain poorly studied, especially when following a major destructive event.
Here, based on measurements of diversity, structure, and resilience, we evaluate the response of subarctic rocky intertidal benthic communities subjected to single and interactive effects of canopy removal, grazer reduction, and nutrient enrichment. This study will improve our understanding of the role of top-down controls, habitat-forming species, and bottom-up forcing in the shaping of community structure. These roles and their interactions are not yet understood for the subarctic ecosystem of the St. Lawrence Estuary. This area was selected as it is exposed to both natural and anthropogenic stressors, including ice-scouring Bourget 1983, Bergeron andBourget 1984), expected increases in water movement (Savard et al. 2008), and eutrophication (Thibodeau et al. 2006, Gilbert et al. 2007). These stressors may affect the abundance and makeup of macroalgae and grazers in benthic intertidal communities. We hypothesize that in addition to a significant individual impact from canopy removal and grazer reduction, univariate and multivariate community characteristics and resilience will be even more affected by synergistic effects when coupled with the stress of nutrient enrichment. Although canopy removal and grazer reduction treatments are disturbances (see Grime 1977 andSousa 1984, where disturbance is related to biomass removal), only the terms stress and stressors will be used hereafter to simplify the reading.

Study site
We conducted the experiment from June 2012 to August 2013 on the south shore of the St. Lawrence Estuary near the municipality of Sainte-Flavie, Quebec, Canada (48°37 0 42.5″ N, 68°11 0 55.7″ W). The site is a flat, mid-intertidal rocky shore having a flora and fauna characteristic of a moderately wave-disturbed environment Bourget 1983, Bergeron andBourget 1984). This subarctic location has a temperature and salinity range from 4-16°C and 24-29 PSU, respectively (Fradette andBourget 1980, Archambault andBourget 1983). The tidal regime is mixed, dominated by semi-diurnal tides having a 3.5-m tidal range on average (see www.tides. gc.ca). Fucus spp. (F. distichus edentatus and F. vesiculosus) structures the natural communities of this mid-intertidal zone, and gastropod grazers (Littorina obtusata and L. saxatilis) and filter feeders (composed of Mytilus edulis, M. trossolus, and hybrids, hereafter referred to as Mytilus spp., see Moreau et al. 2005) are the dominant invertebrates. Ice often covers the shores of the estuary during winter (from mid-December until the end of March); this ice cover provides protection for the underlying biological assemblages against strong variations in water levels, storm waves, and extreme temperatures. However, the ice may also act as an indiscriminate disturbance factor on the flat rock surfaces and exposed crevices through heavy ice-scouring Bourget 1984, McKindsey andBourget 2001).

Experimental design
We used an orthogonal factorial experimental design to evaluate the loss of key species (macroalgal canopy [Ca], two levels; grazer [Gr], two levels; press-type disturbances) and nutrient enrichment ([Nu], two levels; press-type stressor) on the resilience of intertidal benthic communities in terms of their diversity and structure following a complete biomass removal (pulse disturbance) as a pretreatment: In June 2012, all organisms, algae, and sediments were scraped off the plots, and the remaining rock was burned once with a bitumen torch (Fig. 1). We added natural reference plots having no pretreatment to our design as a ninth treatment to evaluate community resilience (Fig. 1). All treatments were replicated four times (n = 4) and assigned randomly to 36 experimental plots (50 9 50 cm) on emergent rocky substrates. We positioned these plots within mature macroalgal bed zones respecting the following criteria: minimum of 80% cover of Fucus spp., homogeneous flat substrate, lack of pools or crevices, similar mid-intertidal height (1.33 AE 0.21 m; no differences of height among treatments, F 8,35 = 0.4878, P = 0.8540), and a minimum of 3 m between plots. We applied the treatments within 50 9 50 cm plots and only sampled the center, a 30 9 30 cm area, so as to avoid edge effects.
The canopy factor had two levels: Natural recruits of fucoids were either left intact (canopy untouched, C+) or completely handpicked (canopy removed, CÀ). We applied such maintenance twice a month from June to early September 2012, once in October 2012, and twice a month from May to the end of August 2013.
The grazer reduction factor had two levels: natural density (G+) or reduced density (GÀ). In the latter, we manually removed invertebrate grazers (periwinkles: L. obtusata, L. saxatilis, L. littorea, Margarites sp., Lacuna vincta; limpets: Testudinalia testudinalis, and isopods: Jaera albifrons) every 9-11 d in 2012 and every 4-8 d in 2013 (the reasons for this change in frequency are presented below). We maintained a reduced abundance of grazers with two physical barriers: A small twisted wire brush (2 cm diameter) screwed on the rock along the plot contour and a thin layer of a natural sticky barrier (5 cm wide; Tree Tanglefoot Insect Barrier, The Tanglefoot Company, Grand Rapids, Michigan, USA) that was renewed twice a month at the outer limit of the brush. Before installing the barriers, we cleared the surfaces, added a small quantity of concrete (when required) on rough surfaces (Poly-Plug Bomix, Daubois, Saint-Leonard, Quebec, Canada), and covered the bare rock or the concrete with epoxy (West Systems, Bay City, Michigan, USA) to ensure the adhesion of the Tree Tanglefoot. We preferred this cageless system to avoid any cage effects on the community, especially as the mesh size would be small (1-2 mm). We tested the method in 2012 with incomplete exclusion procedural controls (n = 4), and we observed no difference with the reference plots for all response variables, with the sole showing the three stress treatments (canopy, grazer, nutrient enrichment; two levels each) following pretreatment (plots scraped to bare rock and then burned), and reference plots left untouched (see Methods for details). Bottom row shows letter codes for treatments with one, two, or three letters representing the quantity of stress applied. Ctrl and Ref represent Control and Reference, respectively. exception in September 2012 for richness, obtained by destructive sampling (cf. Joseph and Cusson 2015). We reduced Littorina spp. abundance in the GÀ treatment for up to four consecutive days using this method; we therefore maintained the plots more frequently in 2013 (i.e., about 40% of the time in 2012 and 66% of the time in 2013 with significantly less individuals in GÀ than in G+; Cimon, Joseph, and Cusson, unpublished data). Moreover, the GÀ treatment significantly reduced the abundance of T. testudinalis and small individuals of Littorina sp. throughout the entire experiment (RM ANOVA with the same design as below (see Table 1, for example) with abundance transformed by square root, between subjects (all periods pooled together) and Gr contrast, respectively F 1,27 = 33.37, P = 0.0001 and F 1,27 = 122.58, P = 0.0001).
Two levels of enrichment were used: ambient (NÀ) and enriched (N+) conditions. Plots were enriched with 200 g of slow-release fertilizer pellets (N-P-K = 14-14-14, Smartcote Plant Prod, Canada) split equally into two screen-mesh bags installed in opposite corners inside the plots but outside the sampling area. Enrichment using slow-release fertilizer pellets has been used in many habitats (e.g., Worm et al. 2000). We replaced fertilizer pellets with washed pebbles in ambient nutrient treatments to consider any effects from the bags (i.e., additional substrate). We changed all bags each month during the sampling periods; we washed the pebbles and replaced the pellets with new ones. We collected, dried, and weighed the pellets at each replacement to estimate the amount of nutrient diffused. We observed an average of 32.5 AE 0.4% weight loss, with a total estimated diffusion of 8.66 AE 0.13 g of total nitrogen per month into each plot. This level of enrichment is similar to a moderate eutrophication in the St. Lawrence Estuary (Gilbert et al. 2007). Pilot tests in the field showed a threefold to sixfold increase in total nitrogen concentrations in water samples from an enriched quadrat compared to the natural concentration of the St. Lawrence Estuary. The tissues of F. distichus edentatus were analyzed for total nitrogen in a concomitant study, and they confirmed that the additional nutrients were assimilated by the algae (see Joseph and Cusson 2015).
We took out the brushes and bags on 20 October 2012, and we reinstalled them between 8 May and 11 May 2013 (a 29-week gap), to prevent the potential loss of material during winter due to ice-scour. There was at least one ice-scouring event during that winter. Reference plots were not the same in 2013 as those in 2012, as the latter were collected in September 2012 and were used as controls in another study (cf. Joseph and Cusson 2015). We selected new reference plots (n = 4) that encompassed our abovementioned criteria, except for canopy percent cover that averaged lower (72%) than our criteria of 80%. In addition, the references did not show any significant differences in community abundance structure between corresponding months from the two growing seasons.

Sampling
We sampled the plots using non-destructive techniques by visually estimating the percentage cover for algae and mussels and by counts for other invertebrates within a 30 9 30 cm quadrat separated into 25 units each representing 4% of the surface area. Each individual (having a size >1 mm) was identified to the lowest taxonomic level possible, usually species, using local guides and taxonomic keys. These monthly inventories took place from June to October 2012 and from May to August 2013: Period 1 (June 2-9; pretreatment), Period 2 (July 1-8), Period 3 (July 31-August 6), Period 4 (August 28-September 4), Period 5 (October 17-20), Period 6 (May 23-24), Period 7 (June 21-23), Period 8 (July 20-23), and Period 9 (August 18-21). We counted invertebrates, other than Mytilus spp., then transformed the counts into percentage cover (e.g., regression to convert density into percent cover with n = 129, R 2 = 0.66 for L. saxatilis and L. obtusata; arbitrary value of 0.25% for each individual of L. littorea, T. testudinalis, and Nereis sp.; 0.1% for each individual of L. vincta and Margarites sp.; 0.01% for each J. albifrons) to have the same abundance unit for all taxa and permit the calculation of diversity estimates. Total percentage cover can exceed 100% as we estimated the cover by taxa. We maintained (e.g., Fucus spp. removal and grazer reduction) the plots after performing the inventories. At the last sampling date (Period 9), right after the usual inventory, we destructively collected all biomass (except crustose species) within the 30 9 30 cm plot. In the laboratory, we filtered the samples through a ❖ www.esajournals.org 1-mm sieve, then counted and weighed all remaining organisms by taxon to the nearest 0.01 g (0.00001 g for tiny organisms). Biomass was converted into energy (kJ) by applying conversion factors from Brey et al. (2010).

Data analysis
We performed all analyses on non-manipulated associated species (unless mentioned); that is, we excluded Fucus spp. and removed grazers, such as Littorina spp., prior to any analyses (see the list above and Appendix S1: Table S1 for further details). We excluded Period 1 from the analyses as it was the application of the pretreatment (June 2012).
To test for simple or interactive effects of the treatments on total abundance (N), species richness (S), Simpson's index of diversity (1Àk), and Pielou's evenness (J 0 ), we used mixed models of repeated measures analysis of variance (RM ANOVA) with one fixed factor (Treatment; nine levels including the reference treatment) and eight repetitions over time (Period, fixed) using the identity of the plots as a random factor. To test our three factors of interest (Ca, Gr, Nu; two levels each), we ran orthogonal contrasts (a total of seven contrasts). If the main factor, Treatment or Period, was significant, we ran Tukey's HSD test to assess the differences between treatment levels or periods. Where the factor Period had a significant interaction with a contrast of interest (e.g., Period 9 Ca 9 Nu), we used comparisons of means (t-test, one test for each period) to determine at which period the difference occurred, and we then ran multiple comparisons of means (t-test, only required for interactions, e.g., Ca 9 Nu) with Bonferroni sequential corrections on the contrast for that specific period. We performed a 1-way ANOVA with the same design (Treatment and seven contrasts) for total biomass, as this variable was only available for the last period. We also used the same RM ANOVA as above to assess recovery time in terms of N, S, 1Àk, and J 0 , and we used comparisons of means (t-test) for each treatment vs. references (eight contrasts, e.g., CG vs. Ref) at each period. We considered that each treatment was recovered for a given variable when no significant difference between that treatment and reference plots was observed until the end of the experiment. We then compared the recovery time of each treatment to the controls to assess whether stressors had any impact on recovery time. Finally, we evaluated the resilience time at our site by comparing controls to references, using all species present, by applying the same method as above (RM ANOVA; with Treatment having only two levels). Normality and homogeneity assumptions were checked by graphical examination of the residuals (Quinn andKeough 2002, Montgomery 2012); only total abundance in percentage cover was transformed by square root prior to analyses. We characterized the interaction type of each significant double or triple interaction effect as either antagonistic, synergistic, additive, or dominant (sensu Côt e et al. 2016) using a calculated 95% confidence limit of the expected additive effect (Darling and Côt e 2008).
To examine the effects on community structure (for both percent cover and biomass, based on Bray-Curtis similarities), we used repeated measures of permutational multivariate analysis of variance (RM PERMANOVA; Anderson et al. 2008) run with 9999 permutations based on the same design as described above (Treatment + seven contrasts, eight repetitions over time). We tested the dispersion of data around a centroid using PERMDISP for each stress at each period. We evaluated recovery as presented above (1-way PERMANOVA; Treatment; eight contrasts of interest). For the latter, only 35 permutations per contrast were possible; as such, we used a Monte Carlo procedure to obtain P values. Abundance data as percent cover and biomass were pretreated using dispersion weighting by Treatment 9 Period (Treatment only for biomass) (Clarke et al. 2006) and transformed by square root according to shade plot method (Clarke et al. 2013), while data were transformed into presence-absence for the effects on compositional community structure. Dummy variables were used to deal with small species numbers (zero-adjusted Bray-Curtis; value of 0.1 for percent cover and 0.02 for biomass; Clarke et al. 2014). We visualized the effects of the treatments using non-metric multidimensional scaling (nMDS) ordinations. We evaluated the contribution of each taxon to the observed similarities/ dissimilarities between treatments using a similarity percentage analysis (SIMPER).
At the end of the experiment, we characterized the stress level among treatments using combined ❖ www.esajournals.org k-dominance curves of species biomass and percent cover data (abundance/biomass comparison or ABC curves). Stable (i.e., undisturbed) communities are usually dominated by the biomass of a few individuals (K-selected), while disturbed communities are dominated by many small individuals (r-selected, opportunistic; Warwick 1986). We used the W statistic, which is a summary statistic of ABC curves ranging from À1 to +1, where +1 represents undisturbed communities, and À1 represents disturbed ones (Clarke 1990). To test the effect on W values, we performed a 1-way ANOVA with the fixed factor Treatment and our contrasts of interest.
Univariate analyses were run using SAS university edition v3.5 (PROC MIXED), while multivariate analyses and ordinations were run using PRIMER+PERMANOVA v.7 (Anderson et al. 2008, Clarke andGorley 2015). We used a significance level of a = 0.05 for all statistical analyses. We applied no corrections to P values to our planned contrasts.

RESULTS
We identified a total of 39 different taxa (see the list in Appendix S1: Table S1) for all inventories with averages (AE SE) of 6.8 AE 0.2 species in 2012 and 8.6 AE 0.2 species in 2013 within the 30 9 30 cm plots. We observed 14 algal taxa, the most abundant being Fucus distichus edentatus, Stragularia clavata, and Porphyra sp., along with 25 invertebrate taxa with the grazers Littorina saxatilis, L. obtusata, Testudinalia testudinalis, and the filter feeder Mytilus spp. being the most common. For both percent cover and biomass data, invertebrates were generally more abundant when the macroalgal canopy was untouched, and ephemeral algae were more abundant when the canopy was removed.

Effects of single and multiple stresses
Abundance and diversity indices.-Our experimental treatments (canopy removal, grazer reduction, and nutrient enrichment) affected the associated community (i.e., non-manipulated species only) abundance and diversity profile, and some of their effects varied in time. As expected, abundance in percent cover, species richness, and Simpson's diversity values were all higher at the end of the experiment (e.g., from periods 2 to 9 with respective averages of 7 AE 3 to 25 AE 3 in percent cover, from 1.4 AE 0.2 to 3.4 AE 0.2 in richness, and from 0.14 AE 0.04 to 0.38 AE 0.03 for Simpson's diversity). We observed the highest values of richness (3.6 AE 0.3) and diversity (0.5 AE 0.03) in periods 8 and 6, respectively. Evenness values, on the contrary, remained rather stable with a maximum and only significant difference between periods 3 and 6 (0.48 AE 0.09 and 0.70 AE 0.04, respectively; F 1, 148 = 9.54, P = 0.0024).
Canopy removal affected abundance in percent cover differently over time (see Period 9 Ca, in Table 1a); its effect was only significant in Period 8 with percent cover values of the associated community more than doubled in the canopyremoved plots (CÀ, Fig. 2a). Canopy removal enhanced species richness values (all periods pooled together; 3.1 AE 0.2 compared to 2.4 AE 0.1 when the canopy was untouched; Table 1b). However, the effect of the interaction between canopy removal and grazer reduction varied over time (see Period 9 Ca 9 Gr; Table 1b). In Period 8, richness values were higher in the canopyremoved and grazer-reduced plots (CÀ GÀ) than with all other combinations of these two stressors, while grazer-reduced plots (C+ GÀ) were lower than with all other combinations of these two stressors (Fig. 2b). The resulting interaction (CÀ GÀ) is a dominance-type interaction as the effect of canopy removal overshadows the effect of grazer reduction (Fig. 2b). Canopy removal significantly increased diversity values (all periods pooled together; Simpson diversity: 0.38 AE 0.02 compared to 0.29 AE 0.02 when the canopy was untouched; Table 1c). The interaction effect between canopy removal and nutrient enrichment varied over time (see Period 9 Ca 9 Nu; Table 1c). In Period 5, diversity values were lower in the presence of both the canopy and nutrient enrichment treatments (C+ N+) than with all other combinations of these two stressors (see Fig. 2c). This result may be due to a significant low evenness observed at the same period in the same plots (cf. Fig. 2c, d). As with grazer reduction, canopy removal dominates the effect of nutrient enrichment for the cumulative impact (CÀ N+) for both diversity and evenness.
At the end of the study (Period 9), we observed different effects when using biomass data compared to percent cover. Canopy removal reduced biomass (Table 2a; Fig. 2e), but it did not affect the abundance in percent cover in this final period. Furthermore, treatment, canopy removal (Ca), and the triple interaction (Ca 9 Gr 9 Nu) all affected richness; a pattern that we did not ❖ www.esajournals.org observe when using percent cover data. The addition of cryptic taxa-mainly amphipods due to sampling methods-can partly explain these differences (see Appendix S1: Table S1 for further details). Canopy removal reduced the total biomass of associated taxa by more than 50% (Table 2a; Fig. 2e). Richness values were affected by the triple interaction (Ca 9 Gr 9 Nu in Table 2b) showing a negative synergism interaction type as we measured lower values in canopyremoved nutrient-enriched (CN) plots than in single treatments of enrichment (N) or grazer reduction (G) (Fig. 2f). Canopy removal reduced diversity (Table 2c; Fig 2g) and induced lower W statistic values ( Fig. 2h; F 1,27 = 4.94, P = 0.0349). Community structure.-Canopy removal significantly affected the associated community structure (Table 3a), beginning at Period 5 (about 135 d; Appendix S1: Table S2) with differences remaining significant (or close to with a P = 0.052 at Period 6) until the end of the experiment (Appendix S1: Table S2, Fig. 3a showing the canopy effect over time and Fig. 3b showing Period 9). Differences were mainly due to more ephemeral algae (e.g., Values are from data in percent cover from Period 8 (a-b), percent cover from Period 5 (c-d), biomass from Period 9 (e-g), and counts and biomass from Period 9 (h). Dark and light gray bars are the respective treatments with C+: canopy untouched; CÀ: canopy removed; G+: grazers untouched; GÀ: grazers reduced; N+: nutrients added; NÀ no nutrients added. See Fig. 1 for details of the different treatments in (f). The number of replicates used to obtain the averages was n = 16 in (a), (e), (g), and (h); n = 8 in (b), (c), and (d); n = 4 in (f). Different letters above the bars indicate significant differences (P < 0.05).
S. clavata and Petalonia fascia) as well as fewer filter feeders (Mytilus spp.) and carnivorous anthozoans (Aulactinia stella) in the CÀ treatment. Grazer reduction or nutrient enrichment did not affect the associated community abundance structure.
Canopy removal significantly increased the dispersion of the associated community structure in both periods 6 and 7 (Appendix S1: Table S3) where average distance among plots was higher in canopy-removed plots (distance from the centroid for canopy-untouched and canopy-removed plots, respectively: Period 6: 31 AE 3 and 43 AE 3; Period 7: 32 AE 3 and 44 AE 2). Grazer reduction significantly increased the dispersion of the associated community structure in Period 8 (Appendix S1: Table S3; Fig. 3c). We observed the highest abundance and richness of ephemeral algae in Period 8. Canopy removal significantly affected the associated community in biomass structure and increased dispersion among plots (Table 3c and in Appendix S1: Table S3; Fig. 3d). We did not observe any significant effect of grazer reduction or enrichment on the associated community biomass structure (Period 9). We observed comparable results with both percent cover and biomass data in compositional structure (all abundances being transformed into presence-absence data; Tables 3d and 3b). Canopy removal altered the species composition within the community throughout the experiment (Table 3b).

Community trajectory over time and community resilience
The total abundance, in percent cover, increased in all treatments over the first summer, yet understandably being at much lower levels in plots where species were removed (i.e., CÀ and GÀ treatments; Fig. 4a). After the winter (Period 6), we observed lower values and a higher variability due to moderate ice-scouring followed by a rapid recovery (Fig. 4a).
Of all treatments, those that included canopy removal (CÀ) had greater average dissimilarities over time when compared to the reference plots (Fig. 4b). Only the canopy-removed treatment, combined with another stress (i.e., CG, CN, and CGN), maintained a significant difference with the references throughout the entire experiment (Fig. 4b). These differences were mainly due to ephemeral algae (e.g., Ectocarpus sp., Ulvaceae, Porphyra sp.) being more abundant in the canopy-removed treatments, and invertebrates (e.g., Mytilus spp., A. stella, Pectinaria gouldii) being more abundant in the references.
Throughout the consecutive periods of the experiment, the trajectories of the average community structure within each treatment were all positively correlated with the trajectories of the controls (see correlations and P values in Appendix S1: Fig. S1). Thus, the period-to-period recovery process was similar to that observed in the controls. We summarized these trajectories using a second-stage nMDS; the nMDS illustrates that treatments having two stressors had more different period-to-period trajectories for community structure (relative to the controls) than single stressor treatments, as expected. Surprisingly, these differences were even greater than for the triple-stressor treatment (CGN; Fig. 5). When all species are considered, controls showed complete recovery after 11 months for all measured variables, although richness recovered after two months, while requiring four months for diversity. At the end of the experiment, after 14 months, controls also recovered in terms of biomass and cryptic species. When using only non-manipulated species, all measured variables in the controls were similar to the reference plots after 11 months (cf. Table 4). Prior to that, controls had lower richness and higher  evenness values compared to the reference plots. In general, richness and species composition required the most time to recover. The treatments of N and GN recovered faster than the controls for total percent cover, richness (N only), evenness, and diversity (see Table 4). All other treatments (C, CG, CN, and CGN) had longer or similar recovery times to the controls.

DISCUSSION
Macrobenthic community resilience depends of the nature and the number of stressors to which the community is subjected. Macroalgal canopy removal, that is, removal of habitat-forming species, had the greatest single effect by modifying in several ways such as diversity indices and community structure, as well as delaying or even preventing resilience. In general, we saw no significant individual impact of grazer reduction or nutrient enrichment. Some non-additive interactions occurred between canopy removal and grazer reduction as well as canopy removal and nutrient enrichment.
The community structure of the triple stress treatment (CGN) differed from that of the reference plots throughout the entire experiment suggesting a slower or even non-occurring resilience when faced with multiple stresses. Canopy removal alone (C) or combined with another stress (CG, CN, CGN) increased the resilience time, as did grazer reduction (G) on its own. Therefore, we partly confirmed our initial hypothesis that canopy removal and grazer reduction would have an individual impact on diversity indices, community structure, and resilience, and that communities would be affected to a greater extent via synergistic effects when coupled with nutrient enrichment stress.

Structure, abundance, and diversity of communities
Of all the main effects, canopy removal most affected the community diversity, abundance, and structure. Canopy-removed plots had a higher abundance (percent cover), richness, and diversity, as well as a different structure in terms of abundance of associated communities at specific times. The differences were mainly due to ephemeral algae being more abundant in the canopyremoved plots, and invertebrates being more abundant in untouched-canopy plots, as expected (see Schiel and Lilley 2011 for similar results). Indeed, Scytosiphon sp., Ectocarpus sp., Petalonia sp., Chordaria sp., and Ulvaceae were absent where the macroalgal canopy was untouched (C+). Although our pretreatment allowed space and light for ephemeral algae to establish, we did not see any delay in fucoid development in C+ plots due to potential competition with ephemeral algae; this observation contrasts with a number of other studies (e.g., Sousa 1979, Korpinen andJormalainen 2008). Therefore, the presence of ephemeral algae is not essential for fucoids to establish and for succession to occur on these sites. Indeed, the recruitment of ephemeral algae occurred mainly during the second summer (2013); this period of the year was too late to result in any competition in the C+ plots. The late recruitment of ephemeral algae could be due to the timing of the experiment. Archambault and Bourget (1983) followed experimentally denuded and reference communities in the St. Lawrence Estuary, and they observed that ephemeral algae were more abundant in experimentally denuded plots only the year after the denuding of the rock and after suffering an ice-scouring event. The reduced presence of invertebrates, when the canopy was removed, can be explained by the reduced habitat availability and increased stress level (e.g., higher temperatures, desiccation, and exposition to ultraviolet radiation; Bertness et al. 1999). Differences in average surface temperature between the reference plots and the canopy-removed plots ranged from 5°to 16°C, temperatures measured using an infrared camera. When using biomass (Period 9), we saw the opposite trend than with percent cover data: lower biomass, richness, and diversity in the canopyremoved plots. The same data transformed into percentage cover no longer showed significant differences in abundance and diversity. This indicates a smaller individual mass of species in canopyremoved plots as suggested by the lower W statistic values that are characteristic of a more affected community (Clarke 1990). The sampling method explains the opposite effect of canopy removal on richness between biomass and percent cover data as cryptic species were added and were generally more present in untouched-canopy plots (see Appendix S1: Table S1 for the list of species). Indeed, when we removed those taxa from the analyses, differences in richness were eliminated, and diversity values increased where the canopy was removed (results not shown). Canopy removal strongly affected the structure in terms of biomass for both location and dispersion (a higher deviation from the centroid when the canopy is absent) indicating a sign of non-resistance to a stress (Warwick and Clarke 1993).
We expected to observe a positive effect of the grazer reduction treatment on ephemeral algae recruitment and growth, as reported by many studies (e.g., Aquilino and Stachowicz 2012). The grazer reduction treatment had no single effect on community diversity indices except for a higher community dispersal at Period 8 (based on percent cover). The late recruitment of ephemeral algae could explain the lack of a single effect by grazer reduction in the first year. Small grazers, which dominated the plots but have a limited grazing capacity, could also explain this lack of effect as also seen by Archambault and Bourget (1983). Yet, there was an interaction effect on richness with canopy removal in Period 8 where two ephemeral species (Ectocarpus sp. and Ulvaceae) were only present in the canopy-removed grazer-reduced plots (C+ GÀ), and anemones (Aulactinia stella) were significantly more abundant in canopy-untouched grazer-untouched plots (C+ G+). Low richness values in reducedgrazer plots having an untouched canopy (C+ GÀ) are not easily explained. The latter generally had fewer encrusting algae, anemones, and barnacles than where grazers were not removed (C+ G+). The canopy removal effect overshadows the effect of grazer reduction when they are combined. This highlights the importance of the macroalgal canopy. On the other hand, the reduced presence of grazers significantly increased the abundance of ephemeral algae over the course of the entire experiment (results not shown). These results indicate that crawling grazers have no or little implication on community succession and structure in this rocky coastal system. Notes: Recovery was achieved when no significant differences (P > 0.05) were observed between a given treatment and the reference plots before the end of the experiment. If not reached by the end of the study, then the recovery time is deemed to be >14 months. Numbers in italic represent cases where recovery happened faster than in the controls, while numbers in bold signify where recovery was slower than observed in the control plots. Note the sampling gap between months 4 and 11 due to winter (i.e., it was not possible to determine recovery time values from months 5 to 10). See Fig. 1 for details regarding the treatment codes.
Primary producers of coastal zones are primarily nitrogen limited (Howarth 1988). Many studies have shown an increase in ephemeral algaewhich are fast-growing species-with high nutrient loading, especially when coupled to low grazing or space availability (e.g., Korpinen et al. 2007). Surprisingly, the single effect of the nutrient enrichment treatment decreased community evenness only in periods 3 and 5; it had, however, an interactive effect with canopy removal that decreased diversity and evenness only when the canopy was present in Period 5 (C+ N+). Those interactions had a combined impact that was strongly dominated by canopy removal, indicating a low importance of nutrient enrichment when facing such high stress. Low evenness in Period 5 for the enriched canopy-present plots (C+ N+) induced a lower diversity as the encrusting algae Stragularia clavata dominated the enriched plots. Nutrient enrichment did not affect total Fucus spp. abundance even though we measured higher nitrogen content in the macroalgae of the enriched treatments (results not shown). Once again, the late recruitment of ephemeral algae in the first year could be responsible for a lack of response to nutrient enrichment, or that the enrichment level was lower than a specific threshold to have a significant effect on our assemblages, or that opportunistic algae were more limited by other abiotic factors on our site. Bertocci et al. (2015) have shown that nutrient enrichment combined with high disturbance intensity may enhance susceptibility to species invasion in rock pools. Our results indicate that our nutrient enrichment level may not have been sufficient to influence succession or community structure (see Bokn et al. 2003 for similar results).

Community resilience
The studied macrobenthic community showed a high capacity of resilience following a pulsetype disturbance. Richness in control plots recovered after only two months, even though the complete recovery of all tested variables required 11 months. Such recovery is relatively fast after the pretreatment removal of all biomass. In contrast, Schiel and Lilley (2011) detected community effects (lower fucoid abundance, lower richness, more ephemeral algae, and different community structure) eight years after a similar disturbance, while Aquilino and Stachowicz (2012) measured a recovery of 18 months for perennial algae percent cover. Subarctic communities recover quickly after ice-scouring events, especially when crevices can be used as refuges (Bergeron and Bourget 1986). Within all our plots, there were no large crevices, although there was some limited relief that increased substrate heterogeneity and promoted recruitment (personal observations). According to O'Leary et al. (2017), community resilience is promoted by the presence of habitatforming species, the supply of new recruits, and favorable physical settings. In our study, the recovery of the habitat-forming species Fucus occurred 11 months after the pretreatment; this coincided with the full recovery of the controls. Finally, gross primary production (see Joseph and Cusson (2015) for methods) did not differ between the controls and the reference plots at the end of the experiment, confirming a recovery in function (results not shown).
The addition of stressors generally slowed down resilience, in agreement with other studies (reviewed in O' Leary et al. 2017). Nutrient enrichment alone (N) or combined with grazer reduction (GN) sped up the recovery of abundance, richness, diversity, and evenness; however, these plots had a different composition and abundance structure than the references, indicating that recovery was incomplete. Indeed, the encrusting algae S. clavata dominated in nutrient-enriched plots, while Mytilus spp. dominated the reference plots. Moreover, nutrient-enriched plots had little cover of habitatforming macroalgae compared to the reference plots. This likely results in the lack of recovery in total biomass and in function (primary production). In contrast, removal of macroalgal canopy (C), grazer reduction (G), and all other combinations reduced resilience capacity. This indicates that even though some stressors did not affect the main tested general community metrics, they can affect stability in terms of recovery time. For this reason, studies need to address as many aspects as possible when examining communities, including stability (recovery and resistance) and function, to reveal a more complete set of responses of the assemblages to multiple stressors.

CONCLUDING REMARKS
Further manipulative studies, such as our study, are required to have a better understanding of the effects of multiple stressors on assemblages. Indeed, the lack of resilience of the triplestressor treatment demonstrates the necessity for studies including more than two stressors and, inevitably perhaps, longer-term studies that last until the combined treatments reach a complete recovery. Stressors affected species differently indicating that the nature of the stress has an influence on community characteristics and their functions. For instance, macroalgal canopy facilitates the recovery of the community. Single or multiple stressors differently affected richness, diversity, evenness, structure, and recovery. This confirms that complementary metrics are needed to document these effects on community stability and persistence. As such, we do not recommend using single diversity measures, such as species richness, especially when habitat-forming species are absent or cannot be included in the analyses, as those measures provide an incomplete picture and can be misleading in terms of the actual state of a community. Rather, we recommend prioritizing the addition of multivariate measures, such as community structure and composition, to include subtle changes in dominance and identity. As our study suggests, multiple stressors may interact differently (dominant, additive, synergistic) on community characteristics and resilience, and the effects cannot be predicted from studies on a single stressor or that assume an additive cumulated impact of multiple stressors.