Analysis
To quantify temporal variation in sessile community composition, we analyzed the percent cover in eleven photographs of each permanent quadrat taken between 29 March 2008 and 14 March 2011. The percent cover of sessile organisms was estimated visually from photographs, using a method developed by [29] and modified by [12]. Taxa were scored only if they were attached to rock or encrusting algae, i.e., epibiotic taxa do not occupy primary space and thus were not quantified. We defined available space as the substratum available for the recruitment and growth of macroalgae and sessile invertebrates [30], which included bare rock, calcified encrusting algae, and non-calcified encrusting algae. Encrusting algae are included in the definition of available space because there is very little bare rock in shallow hard-bottom subtidal habitats, and most invertebrates can overgrow coralline and non-calcified algal crusts [31]. In so doing we assumed that these algal crusts are functionally equivalent, in part for simplicity, but also because the extent to which various species of encrusting algae facilitate [32] or inhibit [33] the settlement of other sessile taxa is poorly understood in this community.

Linear mixed effects models and a model selection approach were used to address the effectiveness of experimental treatments and the primary hypotheses. First, we tested the effects of removal treatment on grazer densities (log-transformed) during the experimental period (March 2009 - March 2010). Second, we tested the effects of grazer removal on temporal variation in invertebrate and macroalgal cover during the experimental period. Third, we hypothesized that grazer removal would result in an increase in clonal ascidians and a concurrent decrease in the cover of available space.

We studied temporal variation in the sessile community before the experiment (28 March 2008-20 March 2009), during the experiment (21 March 2009 - 22 March 2010), and during one year of recovery after the experiment (23 March 2010 - 14 March 2011). Specifically, we quantified the percent cover of sessile invertebrates and macroalgae in quadrats. The comparison of sessile epifauna and macroalgae was of interest because urchins and chitons are studied most often in a macroalgal context [18,19], and because the abundance of these two sessile functional groups is strongly dependent on the orientation of the rock surface [36]. By targeting higher densities of grazers for removal, it could be argued that our approach led to a reduced chance of detecting a treatment effect in the absolute cover of invertebrates or macroalgae, and was thus a conservative design. This is because control quadrats and transects harbored fewer grazers to begin with, and thus the average difference in grazing pressure between control and removal treatments was lower than if replicates for removal had been selected at random. We compared the set of 19 nested, ecologically relevant models for the data collected during the experiment (March 2009-March 2010). All models included site, transect and quadrat as random effects; quadrat at a specific time point was the unit of replication in this analysis. The fixed effects varied between the models, and these details are listed for each model description (Tables 1, 2, S1 and S2). In all models, time was treated as a categorical factor because the dependent variables (e.g., grazer density, percent cover) were not expected to vary linearly through time (for example, due to seasonality). The percent cover of sessile invertebrates was not transformed, but we used a logit transformation for algal cover to correct the nonlinearity observed in a diagnostic plot of residuals against fitted values.

To address the hypothesis that the removal of grazers would change the cover of clonal ascidians and concurrently affect the amount of available space, we quantified annual changes in these functional groups during (March 2009 - March 2010; experiment) and after (March 2010 - March 2011; recovery) the experimental treatment. By targeting higher densities of grazers for removal, we would expect a larger relative change in the benthic community and in this case our experimental approach is biased towards the detection of an effect. All statistical models included site and transect as random effects; quadrat was the unit of replication in this analysis.

Model fit based on maximum likelihood scores was compared using the small sample unbiased Akaike information criteria (AICc), a metric that considers both model fit and complexity (i.e., number of parameters, K). The difference in AICc (Di) between each model and the best model (i.e., lowest AICc) was calculated to emphasize the most plausible models given the data (Di ,2). Finally, Akaike weight (wi), or the relative likelihood of each model, was obtained by normalizing the likelihood across the entire set of candidate models. Details on model selection are provided in Burnham and Anderson [34] and Johnson and Omland [35]. Residuals were inspected visually for normality and homoscedasticity. When necessary, dependent variables were log or logit transformed. Maximum likelihood was used to estimate parameters in all mixed effects models. Statistical analyses were conducted using the packages ‘lme4’ [37] and ‘vegan’ [38] in R 2.14 [39]. Data used in this manuscript are publicly available as Dataset S1.