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.

<br> <em>群落数据</em> 我们于2016年5月至6月，分别在FC与SH的可抵达低潮带冲蚀水道沿线的每个样地中，布设了25个随机设置的0.25 m²样方，以此量化隐蔽群落结构。我们使用不锈钢拉固螺钉与编号标签标记所有样方。通过样点截线法（point intercept method）随时间监测每个样方的盖度百分比与物种丰富度。盖度百分比通过ImageJ（2018 V 1.5版本）对每张照片上15×15网格内的225个样点进行叠加计算得到。每个样点下的类群被鉴定至尽可能低的分类阶元，或归类为裸岩、贝壳碎屑或沙子。这导致部分类群鉴定至物种水平，部分仅鉴定至门水平，因此与同类研究一致，我们将类群归并为更宽泛的分类学类别以开展分析。
<em>实验性干扰</em> 为模拟干扰并研究隐蔽群落的恢复速率，我们在每个样地的25个已标记样方中随机选取10个，使用刮刀移除所有生物后喷涂烤箱清洁剂（如Freidenburg等，2007）。该操作可移除100%的无脊椎动物与冠层/垫状藻类，仅在难以抵达的小型裂缝与缝隙中保留了不足5%盖度的结壳珊瑚藻（encrusting coralline algae）（图2）。剩余的14个样方作为未受干扰的对照样方。当存在冠层藻类（如<em>海带属（Laminaria）</em>物种、<em>无柄海叶藻（Hedophyllum sessile）</em>）时，我们通过目视估算其盖度，随后将其移开以拍摄照片。此类数据会被纳入盖度估算，因此总盖度可超过100%（如Guerry等，2009；如表1所示）。开阔表面恢复速率数据取自前期研究。
中潮带贻贝移除与恢复实验于1996-2000年间开展（Guichard等，2003）。每个实验包含5个重复的0.25 m²样方。在FC和SH样地，我们手动移除样方内的贻贝及其伴生生物。通过1至5个月间隔的摄影监测，记录贻贝对干扰的响应（即盖度变化）。
在另一项低潮带实验中，我们于2011年启动了年度干扰与恢复实验，共分析了7个年度实验的数据（即2011-12至2017-18年度）。每年6月/7月，我们移除FC和SH样地内5个0.25 m²样方中的所有生物。随后以1至5个月的间隔开展为期12个月的摄影监测，记录贻贝、其他固着无脊椎动物以及大型藻类的盖度，之后再次开展移除操作（Menge等，2022）。
<em>环境数据</em> <em>出露时间</em> 我们量化了出露时间（即平均出水时长，下文记为dE）与太阳辐射量（下文记为SI），以此评估栖息地尺度非生物因子（如干燥胁迫、热应激、光照）与群落结构恢复速率之间的关联。出露时间（dE）为基于每个样地的岸线高度与样地所在区域潮汐模式计算得到的月均出水小时数。
为估算出露时间，我们首先使用长距离自调平激光传感器（RL-H4C，拓普康，日本东京）测定每个样地相对于5个基准点的岸线高度。随后使用Trimble GPS系统（美国加利福尼亚州森尼韦尔市，测量精度约±30 cm）测定每个潮间带基准点的海拔（以平均最低低潮面为基准的米数）。接下来，我们通过对取自NOAA潮汐与海流数据库的潮汐数据拟合样条曲线，估算每个样地的出水时长。最后，利用样条曲线函数中每个样地的高度得到其每日出水时长占比的估算值，并按月对这些数据取平均，得到月均出露时间。我们意识到该方法仅能提供<em>估算</em>的平均出露时间，且天气驱动的变异会影响实际的真实出露时长（Harley & Helmuth, 2003）。出露时间反映了干燥胁迫，这通常是潮间带生物面临的强烈生理应激（Helmuth等, 2006; Williams等, 2013; Flores等, 2015）。
<em>太阳辐射量</em> 太阳辐射量（SI）与干燥胁迫、热应激以及光合作用相关，因此是潜在的关键非生物因子（Harley & Helmuth, 2003）。为测定SI，我们使用了Solar Pathfinder（美国田纳西州林登市）配套Solar Pathfinder Assistant Base软件（https://www.solarpathfinder.com）。该设备为凸面塑料穹顶，可反射每个样地360°视野内所有遮挡阳光的障碍物的全景图像（如悬垂岩架、贻贝床等）。通过导入穹顶照片并勾勒障碍物轮廓，软件可基于指定样地位置与样地近似方位角（两个样地均为正北偏西约167°），估算出以kWh/m²/天为单位的辐射量，且该估算值为月份与当日时段的函数。该测量未考虑样地被淹没时的光衰减效应，因此仅能近似实际的太阳辐射量。
<em>其他指标</em> 为记录样地间的其他环境差异，我们从俄亥俄州立大学PISCO监测项目中获取或采样得到了样地尺度的叶绿素a、营养盐（NO3+NO2）、颗粒有机氮（PON）、pH、碱度以及气温与水温数据。
<em>数据分析</em> 数据分析采用了搭载PERMANOVA+插件的PRIMER-e（《多元生态研究普利茅斯流程》，V7.0.21版本）、JMP软件（SAS研究所，V16.0.0，2021）以及R语言（v4.0.0）。我们通过对数据进行平方根转换并构建Bray-Curtis相似性矩阵（Bray-Curtis similarity matrix），量化了各样方间的群落相似性（McCune等, 2002; Clarke & Gorley, 2015）。
<em>群落结构的大尺度效应</em> 为探究未受干扰的冲蚀水道群落是否反映了大尺度过程（问题1），我们采用置换多元方差分析（PERMANOVA，Anderson等, 2008），检验了样地与时间对对照样方群落结构的影响。实验采用重复测量设计，将嵌套于样地的样方作为随机效应，将样地、时间（以采样周期作为有序因子）及其交互作用作为固定效应。我们采用999次置换的多元离散置换检验（PERMDISP），检验PERMANOVA中样地与时间点间的显著效应是否由方差异质性驱动（Anderson等, 2008; Viejo, 2009; Menge等, 2015）。基于距离的相似性百分比分析（SIMPER）可识别驱动样地与时间分离的分类类群（Anderson等, 2008; Viejo, 2009）。非度量多维标度（nMDS，non-metric multidimensional scaling）可可视化样地与时间如何与群落相似性相关联。我们在nMDS图上叠加向量，以展示在SIMPER分析中对样地间相异性贡献超过5%的类群。
<em>微生境尺度环境因子对群落结构的影响</em> 为检验微生境是否会影响未受干扰的冲蚀水道群落结构（问题2），我们采用PERMANOVA与PERMDISP分析了出露时间与太阳辐射量的效应（样地为固定效应，出露时间与SI为协变量）。鉴于我们关注二者对群落结构的累积效应，且对照样方的群落随时间变化较小（详见下文），我们首先对每个样方的物种盖度随时间取平均，并对2016年5月至2017年8月期间每个样方的月均出露时间与SI取平均。
<em>样地对群落恢复的效应</em> 我们比较了每个样地中移除样方与对照样方随时间变化的群落结构，以此探究隐蔽群落的恢复速率与模式是否与样地间大尺度过程（浮游植物丰度、补充速率）相关（问题3），分析采用PERMANOVA。PERMDISP检验了样地与处理组间的方差同质性。为可视化每个群落随时间的恢复过程，我们基于nMDS第一轴绘制了样地×处理组的质心，并使用向量轨迹连接各时间点。鉴于SH与FC的对照群落在nMDS1轴上占据不同位置，我们将群落恢复定义为其到达对应样地对照群落的nMDS1轴位置。我们采用广义线性混合模型（GLMM，generalized linear mixed model），检验了样地、处理组、时间（以有序因子表示）及其交互作用是否影响群落恢复（即nMDS1轴位置），其中将嵌套于样地与处理组的样方作为随机效应。通过Tukey多重比较事后检验，我们确定了移除样方与对照样方群落统计学上相似（表明已恢复）或存在差异（表明未恢复）的时间节点。我们采用双向SIMPER分析，探究了与特定样地和时间相关的类群。单向SIMPER分析识别了第0天与第468天（SH样地）、第0天与第469天（FC样地）间群落结构的处理组差异，以此评估移除组与对照组间的类群恢复情况。
为探究隐蔽群落多样性的时空变化，我们采用GLMM。多样性恢复评估采用Tukey多重比较事后检验，比较了时间×样地×处理组组合下的平均多样性。开阔表面群落恢复速率采用协方差分析（ANCOVA，analysis of covariance）进行分析，其中样地为固定效应，移除后的天数（低潮带实验中额外加入年份作为协变量）为协变量。分析前对数据进行了转换：贻贝盖度与其他无脊椎动物盖度采用ln(x+1)转换，大型藻类采用反正弦转换。
<em>环境因子效应</em> 为检验恢复速率是否受出露时间或辐射量影响（问题4），我们采用R语言的‘adehabitat’包（通常用于量化运动轨迹，Calenge, 2006）计算了nMDS1与nMDS2坐标的轨迹向量长度。我们通过GLMM检验了样地、出露时间与辐射量是否影响恢复速率（即向量长度）。我们采用协方差分析检验了环境指标的样地间差异，以及这些指标如何随实验天数变化，其中实验天数作为协变量。