Community Structure and Recovery in Cryptic Intertidal Communities Reflects Dynamics on Open-Surface Communities

<br> <em>Community Data </em> We quantified cryptic community structure in 25 haphazardly installed plots (0.25 m2) at each site in May-June 2016 along the length of accessible low zone surge channels at FC and SH, respectively. We marked and labeled plots with stainless steel lag screws and a numbered tag. We determined percent cover and species richness in each plot over time using the point intercept method. Percent cover was estimated by overlaying 225 points in a 15 x 15 grid per photo using ImageJ (2018 V 1.5). Taxa under each point were identified to the lowest possible taxonomic resolution or identified as bare rock, shell debris, or sand. This resulted in some taxa identified to species and others to phyla, so as in similar studies we lumped taxa into broader taxonomic categories for analysis. <em>Experimental disturbances</em> To mimic disturbance and study recovery rates in cryptic communities, we cleared n = 10 randomly chosen plots out of the 25 marked plots at each site by removing all biota with scrapers then spraying the plot with oven cleaner (e.g., Freidenburg et al., 2007). This removed 100% of invertebrates and canopy/turf algae and removed all but &lt;5% cover of encrusting coralline algae in inaccessible small cracks and crevices (Fig. 2). Remaining plots (n = 14) were undisturbed controls. When canopy algae (e.g., <em>Laminaria </em>spp., <em>Hedophyllum sessile</em>) were present, we estimated its percent cover by eye, then moved it aside for the photo. These data were added to percent cover estimates, so total cover could sum to &gt; 100% (e.g., Guerry et al., 2009) (e.g., Table 1). Open-surface recovery rates were obtained from prior research. Mid zone mussel removal and recovery experiments were run from 1996-2000 (Guichard et al., 2003). Each experiment consisted of 5 replicate 0.25 m2 quadrats. Plots were manually cleared of mussels and associated biota at FC and SH. Mussel responses to disturbances (i.e., changes in % cover) were monitored photographically at 1-to-5-month intervals. In a separate low zone experiment, annual disturbance and recovery experiments were initiated in 2011. We analyzed results of 7 annual experiments (i.e., 2011-12 to 2017-18). Biota were cleared each June/July from five 0.25 m2 plots at FC and SH. Percent cover of mussels, other sessile invertebrates, and macrophytes were monitored photographically at 1-to-5-month intervals for 12 months, then removals were reinitiated (Menge et al., 2022). <em>Environmental data</em> Emersion. We quantified emersion time (average time out of water, hereafter dE) and solar irradiation (hereafter SI) to assess the association between habitat-level abiotic factors (e.g., desiccation, thermal stress, light) and community structure recovery rates. Emersion (dE) was the monthly average number of hours out of water based on the shore level of each plot and the tidal patterns at the site. To estimate emersion, we first measured the shore level of each plot in relationship to five benchmarks using a long-range self-leveling laser sensor (RL-H4C, Topcon, Tokyo, Japan). We then measured the elevation of each intertidal benchmark (m above mean lower low water) using a Trimble GPS system (Sunnyvale, California, USA), which is accurate to approximately ± 30 cm. Next, we estimated the time out of water for each plot at each site by fitting a spline curve to tide data taken from NOAA Tides and Currents. Finally, daily estimates of the percent time each plot was out of the water were obtained using the height of each plot in the spline curve function and then averaging these data by month. We recognize that this method provides <em>estimated</em> average emersion, and that weather-driven variation affects the actual amount of true emersion (Harley &amp; Helmuth, 2003). Emersion reflects desiccation, often a strong intertidal physiological stress (Helmuth et al., 2006; Williams et al., 2013; Flores et al., 2015). Solar Irradiation. Solar irradiation (SI) is related to desiccation, heat stress, and photosynthesis and thus can be a potentially important abiotic factor (Harley &amp; Helmuth, 2003). To measure SI, we used a Solar Pathfinder (Linden, Tennessee, USA) coupled with Solar Pathfinder Assistant Base Software (https://www.solarpathfinder.com). The Pathfinder is a convex plastic dome that reflects a panoramic image of any obstacle in the 360° horizon that will block the sun in each plot (e.g., overhung ledges, mussel beds, etc.). By importing a photograph of the dome and outlining the obstacles, the software estimates irradiance in kWh/m2/day as a function of month and time of day based on the specified location (i.e., plot) and approximate azimuth of the site (~167° from N for both sites). This measurement did not account for light attenuation when plots were submerged so only approximated actual SI. Other metrics. To document other environmental differences between sites, we sampled or obtained site-level data on chlorophyll-a, nutrients (NO3 + NO2), particulate organic nitrogen (PON), pH, alkalinity, and air and water temperature from the OSU PISCO monitoring program. <em>Data analyses</em> Data analysis employed PRIMER-e with PERMANOVA+ (Plymouth Routines in Multivariate Ecological Research, Version 7.0.21), JMP software (SAS Institute Inc., Versions 16.0.0, 2021) and R (v. 4.0.0). We quantified community similarity among plots using square root-transformed data and a Bray-Curtis similarity matrix (McCune et al., 2002, Clarke &amp; Gorley, 2015). Large-scale effects on community structure. To determine if large-scale processes were reflected in undisturbed surge channel communities (Question 1), we tested the effect of site and time on community structure in control plots using permutational multivariate analyses of variance (PERMANOVA, Anderson et al., 2008). We used a repeated measures design with plot nested within site as a random effect, and site, time (as an ordinal factor using sampling period), and their interaction as fixed effects. Permutational tests of multivariate dispersion (PERMDISP) with 999 permutations tested if significant effects among sites and time points in PERMANOVA were driven by variance heterogeneity (Anderson et al., 2008; Viejo, 2009; Menge et al., 2015). Distance-based similarity percentage analysis (SIMPER) indicated which taxonomic groups were driving site and time separations (Anderson et al., 2008; Viejo, 2009). Non-metric multidimensional scaling (nMDS) enabled visualization of how site and time were associated with community similarity. We used vector overlays on nMDS plots to show taxa contributing &gt; 5% to among-site dissimilarity in SIMPER analyses. Microhabitat-level environmental effects on community structure. To test if microhabitat influenced undisturbed surge channel community structure (Question 2), we analyzed the effect of emersion and Solar Irradiation using PERMANOVA and PERMDISP (site = fixed, emersion and SI were covariates). Since we were interested in their cumulative effects on community structure, and because communities varied little over time in control plots (see below), we first averaged the species cover in each plot over time and averaged the emersion and SI for each plot by month from May 2016 and August 2017. Site effects on community recovery. We compared community structure over time in cleared and control plots at each site to examine if rates and patterns of cryptic community recovery were related to large-scale between-site processes (phytoplankton abundance, recruitment rates) (Question 3) using PERMANOVA. PERMDISP tested for homogeneity of variances among sites and treatments. To visualize the recovery of each community over time, we plotted site x treatment centroids using nMDS axis 1 with vector trajectories connecting time points. Since SH and FC control communities occupied different positions along nMDS1, we considered communities recovered when they reached the nMDS1 position of their respective control plots. Using GLMM, we tested if site, treatment, time (as an ordinal factor), and their interactions affected community recovery (nMDS1 position), including plot nested with site and treatment as a random effect. Tukey's pairwise post-hoc tests determined times when removal and control plot communities were statistically similar (indicating recovery) or different (indicating non-recovery). We used two-way SIMPER analyses to investigate which taxa were associated with a given site and time. One-way SIMPER identified between-treatment differences in community structure between day 0 versus day 468 [SH] and 469 [FC) to assess taxon recovery in removals vs. controls. To investigate spatio-temporal changes in cryptic community diversity, we used GLMM. Assessment of diversity recovery used Tukey’s pairwise post hoc test to compare average diversity by time x site x treatment. Open surface community recovery rates were analyzed using ANCOVA (analysis of covariance) with site as fixed factors and days since clearance (plus year in low zone experiments) as covariates. Data were ln (x+1)-transformed (mussel cover, other invertebrate cover) or arcsin-transformed (macrophytes) for analysis. Environmental effects. To test if recovery rate was affected by emersion or irradiance (Question 4), we calculated trajectory vector lengths of nMDS1 and nMDS2 coordinates using the “adehabitat” package in R, typically used to quantify movement (Calenge, 2006). GLMM tested if site, emersion, and irradiance influenced recovery rate (vector length). We tested for between-site differences in environmental measures and how these varied with days of the experiment using analysis of covariance, with days as the covariate.