Data Processing:
We conducted separate statistical analyses for surfgrass and mussel tide pools for all response variables due to the difference in effects that each foundation species has on the ecosystem. Change in foundation species loss between the before and after removal periods varied continuously across treatment groups from -27.5% to 100% change in cover in the surfgrass tide pools and -10.4 to 98.5% change in cover in the mussel pools (Fig. S1, Fields & Silbinger 2022), with negative values indicating there was an increase in foundation species cover between time points. Therefore, we investigated foundation species loss as a continuous gradient rather than using categorical ANOVA-style design typical for BACI studies (Smith 2002).

Multivariate plots and analyses were used to assess how community composition and biogeochemistry parameters changed between time periods. We used a principal coordinate analysis (PCoA) in the vegan package (Oksanen et al. 2007) to visualize how the gradient of percent foundation species loss altered the change (after – before removal periods) in the cover of sessile organisms and counts in mobile organisms by functional group. Mobile community data were square root transformed to reduce the effect of rare values. One surfgrass control tide pool (Pool ID 18) was removed from the mobile community analysis because of outlier values of limpets and littorine snails. A PERMANOVA was used to test if foundation species loss altered functional community composition, with tide height and tide pool size included as covariates. For the biogeochemistry data, principal component analyses (PCA) were used to reduce the dimensionality and visualize the multivariate data. Data was centered and standardized before visualizing the PCA. Specifically, shifts between the before and after-removal period in biogeochemistry (DO, nutrients, and pH) and temperature mean, maximum, and variance of n = 8 (four day and four night) time points were visualized to observe differences for local tide pool conditions and the adjacent ocean from water sampling events.

We used general linear models to test the effect of foundation species loss on change (after – before removal periods) in species richness, percent maximum light, and maximum temperature in tide pools, with tide height and tide pool volume as covariates. Model assumptions for general linear models were verified by investigating residual plots of the model for homogeneity of variance and normality and data that did not meet assumptions were transformed. Multi-collinearity was tested with Pearson correlations between each predictor. All data was processed and analyzed using R statistical program (v 4.0.2; R Core Team 2020).

We conducted a piecewise structural equation model (SEM) using PiecewiseSEM 2.1.0 (Lefcheck 2016) understand how mussel and surfgrass percent loss and tide pool physical characteristics affected community changes, physical environmental conditions, biogeochemistry, and ecosystem metabolism of tide pool communities. Piecewise SEMs combine multiple linear models into a single causal framework (Shipley 2009, Lefcheck 2016). Because the model is pieced together with local estimation, it is more flexible than traditional SEMs and allows for inclusion of non-normal distributions, random effects, nested models, and smaller sample size within the causal model network (Lefcheck & Duffy 2015). Specifically, we tested direct and indirect pathways of surfgrass and mussel percent loss on ecosystem metabolism (NEC and NEP) mediated by changes in micro/macroalgae cover, temperature, nutrients (dissolved nitrogen to phosphate ratio), pH, and physical parameters of the tide pool. We conducted a multigroup analysis to test if each path varied by day and night, as expected from previous tide pool studies (Kwiatkowski et al. 2016, Bracken et al. 2018, Silbiger & Sorte 2018, Wolfe et al. 2020). The fit of each model was determined with a Shipley’s test of d-separation, which verifies there are no missing inferences from the model using a Fisher’s C statistic, where p>0.05 indicates that the model is a good fit for the data (Shipley 2009). However, since our study’s ratio of total number of samples (16) to the number of variables (9) is 1.78 versus the recommended 5 (Grace et al. 2015, Lefcheck 2016), the smaller sample size could influence the goodness of fit of the model. Path coefficients were standardized allowing for comparisons of the magnitude of effect among groups.

For all model components, we averaged values over the low tide period (n=3 - 4 values per tide pool for pH, nutrients, and temperature; n= 2 - 3 values for NEC and NEP), and then calculated the difference between the before and after removal time period (n = 16 pools for the surfgrass model and n =15 pools for the mussel model). Positive values indicate an increase, zero indicates no change, and negative indicates a decrease in that parameter after the removal period. For the community composition components of the model, we used the sessile community metrics that changed the most between the before and after removal period, including, surfgrass, mussel, and non-surfgrass, non-calcifying producers (all fleshy micro and macroalgae). We excluded coralline algae from the analysis due their different functional roles within the tide pool (e.g., calcification and production and to reduce the number of pathways in our analysis). Tide pool size, represented as volume, and tide pool height was used as a covariate in all models because physical parameters can affect temperature and biogeochemistry (Legrand et al. 2018, Wolfe et al. 2020).

Due to the collinearity between nutrient species, the dissolved nitrogen to phosphate ratio ([nitrite + nitrate + ammonium]: phosphate) was used as the measurement for nutrients in the SEM models. Because temperature and light were highly correlated and change in light was zero in most pools at night, only maximum temperature was used with the SEM framework (Fig. S3, Fields & Silbinger 2022). We used existing knowledge of intertidal systems to create hypothesized paths within each of the models (Grace 2008, Grace et al. 2012):

Micro/macroalgal cover ~ Foundation species loss + Volume + Tide height [1]

Model 1 represents the hypothesis that loss of surfgrass and mussel cover would increase producer cover by increasing free space for fast colonizers, with the volume of the pool and tide height as a covariates (Dethier 1981, 1984).

Maximum Temperature ~ Foundation species loss + Volume + Tide height [2]

Model 2 tests the understanding that temperature will increase with surfgrass and mussel loss (Stephens & Bertness 1991, Shelton 2010), with size of the pool and tide height as covariates.

N:P ~ Foundation species loss + Volume + Tide height [3]

Model 3 represents the hypothesis that the N:P ratio will increase with surfgrass loss because the loss of dominant producers will lead to less nitrogen uptake (Terrados & Williams 1997, Ramirez-Garcia et al. 2002). We hypothesized that the N:P ratio will decrease with mussel loss due to the mussels’ ability to recycle nitrogen (Nielsen 2003, Bracken & Nielsen 2004, Pfister 2007, Pfister & Altabet 2019). N:P ratio may be affected by the size of the tide pool due to diffusive properties of tide pools (Hurd 2000) and tide height due to emersion time.

NEP ~ Maximum temperature + Micro/Macroalgal Cover + N:P + Tide Height;

NEP~ N:P + Micro/Macroalgal Cover + Tide Height [4]

Model 4 represents the different mechanisms that influence NEP. NEP is increased by producer cover, temperature, and nutrients (N:P) and affected by tide height (Pfister 2007, Kwiatkowski et al. 2016, Takeshita et al. 2016, Bracken et al. 2018, Duarte & Krause-Jensen 2018, Silbiger & Sorte 2018, Pfister & Altabet 2019, Wolfe et al. 2020). Surfgrass and mussel loss were both correlated with most of the components of the NEP model and were therefore excluded from the NEP model. Since volume was already accounted for in the calculation of NEP, it was removed from the model equation. Maximum temperature and N:P were correlated within mussel pools, so N:P ratio was used because there is more background knowledge on mussels’ effect on production via altering nutrients (Pfister 2007, Pfister & Altabet 2019).

pH ~ NEP + Foundation species loss + Volume + Tide height [5]

Model 5 represents our understanding that there is a positive relationship between NEP and pH (Silbiger & Sorte 2018). Production increases pH, while respiration decreases pH due to the uptake and release of CO2 respectively. Foundation species loss, size of tide pool, and tide height may also affect pH. Foundation species loss was also included within the pH model as a missing link determined by tests of directed separation.

NEC ~ pH + Maximum temperature + Tide Height [6]

Model 6 represents our understanding that NEC increases with pH and temperature and may be affected by the tide height (Kwiatkowski et al. 2016, Silbiger & Sorte 2018, Wolfe et al. 2020). Since volume was already accounted for in the calculation of NEC, it was removed from the model equation. Foundation species loss was correlated with both pH and maximum temperature and therefore was excluded from the model.
R statistical program (v 4.0.2; R Core Team 2020)

BCO-DMO Processing Description:
- Adjusted field/parameter names to comply with BCO-DMO naming conventions
- Added a conventional header with dataset name, PI names, version date