Resistance of rocky intertidal communities to oceanic climate fluctuations

resemblance matrices. We tested the relative importance of temporal (year), spatial (e.g. zone, site, cape), or oceanic climate (climate metrics and years since SSWD) factors on community structure using multiple distance-based permutational analyses of variance (PERMANOVA) [91,92]. We used 999 permutations under a reduced model, conditional type 2 sums of squares, and employed multivariate pair-wise comparisons among levels of significant factors. Percent contribution of each model term to the overall fit was determined by dividing the term’s component of variation to total model variation [89]. We used similarity percentages analysis (SIMPER) to identify the functional groups substantially contributing to dissimilarity among or similarity within levels of a factor [88]. We visualized community separation for each term using non-metric multidimensional scaling (nMDS) plots with vector overlays for those functional groups having the strongest correlations to nMDS cloud separation. We used PERMDISP to determine if each significant term in PERMANOVAs met the assumptions of multivariate homogeneity [89,93] (see SI for details). Community structure among zones. Zone effects on community structure were analyzed by testing main and interacting effects of intertidal zones (high, mid or low), capes, sites (nested within capes) and years, with cape, site, and year as random factors using PERMANOVA. Since the dataset was large and thus computer processing was prohibitively time-consuming, we averaged each functional group by transect rather than including each quadrat separately.

Unsurprisingly, community structure was overwhelmingly driven by intertidal zone (Table 1, see Results). Since our primary focus was on temporal variation, we performed further analy- ses separately by zone. Community associations with space and time. To investigate the relative importance of tem- poral and spatial (among-site) variation on community structure, we tested the effects of cape, site (nested within cape), year, and their interactions (all as random effects) on community structure by zone using PERMANOVA. Here again the large size of the low zone dataset pre- vented quadrat-level model runs for PERMANOVAs and made visualizing trends in nMDS plots difficult. Thus, we averaged each functional group by transect for all low zone PRIMER analyses. Mid and high zone analyses were performed at the quadrat level since transect level averaging severely reduced model degrees of freedom. Relative importance (see above for ini- tial percent contribution calculation for each term) of the spatial and temporal components of variation were estimated as: spatial % contribution = S(cape + site + 0.5*(year x cape) + 0.5* (year x site), and temporal % contribution = S(year + 0.5*(year x cape) + 0.5*(year x site] [89]. Community associations with oceanic climate. We used PERMANOVA to investigate if changes in community structure over time in each zone were correlated with the focal climate indices that the organisms would have been subjected to over the prior year (i.e. between sur- veys). Because environmental data were on different scales, we normalized them (ranging

from -2 to +2) before analysis using Primer 7 [88] with the PERMANOVA+ add-on [89]. We inspected draftsman plots and found no substantial collinearity among the 8 environmental variables that might interfere with models (All R2 values < 0.4, except R2 = 0.523, 0.456, 0.410, 0.404 for PDO vs. SST, PDO vs. ENSO, NPGO vs. Year since SSWD and BEUTI vs. SD BEUTI, respectively). In each model, we included cape and site (nested within cape) as random factors and the eight climate indices as fixed factors (mean PDO, mean NPGO, mean ENSO, mean BEUTI, SD BEUTI, maximum air temperature, mean sea surface temperature, and years since SSWD). To identify candidate functional groups that may have been associated with environmental changes, we used SIMPER to identify the 10 most temporally variable func- tional groups by zone (highest average dissimilarity among years). We then analyzed their sen- sitivity to annual environmental fluctuations with univariate generalized linear mixed models (GLMMs) using the glmer function in the lme4 package in R [94]. Site was random and climate indices were fixed. We used a Bonferroni correction threshold of p <0.00625 (8 terms). For PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 7 / 34