Resistance of rocky intertidal communities to oceanic climate fluctuations

RESEARCH ARTICLE Resistance of rocky intertidal communities to oceanic climate fluctuations Sarah A. GravemID*, Brittany N. Poirson¤a, Jonathan W. Robinson¤b, Bruce A. Menge Department of Integrative Biology, Partnership for Interdisciplinary Studies of Coastal Oceans (PISCO), Oregon State University, Corvallis, Oregon, United States of America ¤a Current address: Washington Department of Natural Resources, Olympia, Washington, United States of America ¤b Current address: Washington State University Snohomish County Extension, Snohomish, Washington, United States of America * GravemS@oregonstate.edu Abstract A powerful way to predict how ecological communities will respond to future climate change is to test how they have responded to the climate of the past. We used climate oscillations including the Pacific Decadal Oscillation (PDO), North Pacific Gyre Oscillation, and El Niño Southern Oscillation (ENSO) and variation in upwelling, air temperature, and sea tempera- tures to test the sensitivity of nearshore rocky intertidal communities to climate variability. Prior research shows that multiple ecological processes of key taxa (growth, recruitment, and physiology) were sensitive to environmental variation during this time frame. We also investigated the effect of the concurrent sea star wasting disease outbreak in 2013–2014. We surveyed nearly 150 taxa from 11 rocky intertidal sites in Oregon and northern California annually for up to 14-years (2006–2020) to test if community structure (i.e., the abundance

of functional groups) and diversity were sensitive to past environmental variation. We found little to no evidence that these communities were sensitive to annual variation in any of the environmental measures, and that each metric was associated with < 8.6% of yearly varia- tion in community structure. Only the years elapsed since the outbreak of sea star wasting disease had a substantial effect on community structure, but in the mid-zone only where spatially dominant mussels are a main prey of the keystone predator sea star, Pisaster ochraceus. We conclude that the established sensitivity of multiple ecological processes to annual fluctuations in climate has not yet scaled up to influence community structure. Hence, the rocky intertidal system along this coastline appears resistant to the range of oce- anic climate fluctuations that occurred during the study. However, given ongoing intensifica- tion of climate change and increasing frequencies of extreme events, future responses to climate change seem likely. Introduction Many species on Earth are sensitive to climate change [1–5]. However, the response of entire ecological communities and ecosystems to the changing climate is underexplored. PLOS ONE PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 1 / 34 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Gravem SA, Poirson BN, Robinson JW, Menge BA (2024) Resistance of rocky intertidal communities to oceanic climate fluctuations. PLoS ONE 19(5): e0297697. https://doi.org/10.1371/ journal.pone.0297697

Editor: Vitor Hugo Rodrigues Paiva, MARE – Marine and Environmental Sciences Centre, PORTUGAL Received: June 28, 2023 Accepted: January 11, 2024 Published: May 29, 2024 Copyright: © 2024 Gravem et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All data are publicly available at https://datadryad.org/stash/share/ TqHULHzm5D-gNHicsFv7rZjegrQ7EcgV7fyL0T4k_ fc. Funding: This work was funded by the National Science Foundation Long-Term Research in Environmental Biology (LTREB) Program (DEB10- 50694, DEB15-54702, DEB20-50017) and partially supported by funds from the David and Lucile Packard Foundation, the Andrew W. Mellon Foundation, the Wayne and Gladys Valley

Determining the responses of complex ecological communities comprised of a network of interacting primary producers, herbivores, predators, omnivores, and detritivores is typically untenable. However, some complex communities are easier to study than others. For example, rocky intertidal community members are relatively small, typically have slow-moving or sessile habits, mature and reproduce relatively quickly, and have a compact habitat that enables research at a community scale [6,7]. One approach for assessing community responses to global climate change is to amass long-term datasets of whole communities and explore relationships with past environmental variation [4,8]. Here, we report results of a 14-year study of responses of rocky intertidal com- munities in the northern California Current Large Marine Ecosystem (CCLME, coastal ecosys- tem from Baja California, Mexico, to British Columbia, Canada) to fluctuations in large-scale ocean climate and environmental patterns. These include the Pacific Decadal Oscillation (PDO), North Pacific Gyre Oscillation (NPGO), El Niño-Southern Oscillation (ENSO), as well as upwelling regime and variation (intermittency) in upwelling, maximum daily intertidal air temperatures, and average daily intertidal sea surface temperatures. We also tested the effects of the sea star wasting disease epidemic, which caused severe coast-wide declines in the key- stone predator Pisaster ochraceus in 2013 and 2014 [9,10]. Recent evidence of increasing instability

Our recent work shows some worrisome ‘early warning signals’ of impending change in this rocky intertidal system [11,12]. For example, annually repeated disturbance experiments to track recovery rates in the low zone at several sites in Oregon revealed increasing community instability (more variance in disturbed and undisturbed plots over time) and diminishing resil- ience (slower recovery from disturbance over time [11]). Further, overall declines and/or increasing variability in several key ecological metrics have occurred, including ecological sub- sidies to the nearshore from offshore (i.e. phytoplankton and invertebrate recruitment), mussel performance (body condition and reproductive output), and sea star performance (reproduc- tive output and predation rates [12]. Climate fluctuations may provide insight into drivers of community instability Oscillating oceanic phenomena like the PDO (20–30 yr cycle), NPGO (7–15 yr cycle), and ENSO (3–7 year cycle), and environmental conditions like upwelling and air or sea tempera- tures have well-known impacts on some marine communities [13–18]. Higher PDO indices generally indicate warmer sea surface temperatures [13]. NPGO is orthogonal to PDO and higher indices generally indicate windier conditions and higher primary productivity [19,20]. In the CCLME, ENSO results in alteration between warmer, less productive El Niño years and cooler, more productive La Niña years [15,21,22]. Upwelling is caused by alongshore wind stress [23] that draws deep, nutrient-rich water up to the surface, and is usually strongest in

June and July northern California and Oregon [24]. Variability in upwelling is also important, where intermittently strong upwelling (typical of central Oregon) provides both enough nutri- ents and enough wind relaxation for phytoplankton blooms to develop nearshore, while persis- tently strong upwelling winds (typical of southern Oregon and northern California) transports the blooms offshore [16,25]. Further, these environmental patterns and atmospheric climate are intricately linked, so global climate change is expected to alter them [17,26–29]. For exam- ple, the intensity and relative ecological dominance of NPGO versus PDO may be increasing [30,31], extreme El Niño and La Niña events are expected to increase [32]. On the other hand, PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 2 / 34 Foundation, the Kingfisher Foundation, and the NSF Biological Oceanography Program (OCE07-26983, OCE10-61233, OCE14-48913, and OCE17-35911). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

expectations for upwelling are mixed [12,33,34]. Warming of both air and sea temperatures is ongoing and is expected to continue with climate change [12,35]. The oscillations in the metrics above provide a regime of past climactic variability that may help us understand marine community sensitivity or resistance to future climate change. For example, if communities were notably sensitive to environmental fluctuations over past decades, we may assume that communities will respond readily to climate change in the near future [36,37]. Conversely, if communities remain unchanged despite repeated climate fluctu- ations or perturbations, we may assume that these communities are more resistant [38]. While this line of inquiry may provide some insight for the near future, we admittedly do not know whether the relationships between increasing climate stress and community structure are lin- ear. For example, communities may exhibit tipping points [39–41] where even gradual envi- ronmental change can lead to sudden community responses. Also, climate-related, prominent drivers of community change such as extreme events like heat waves, disease outbreaks, and natural disasters are clearly increasing in frequency [42–45]. Though predicting the effects of such episodic events is challenging, we may be able to infer whether communities will be responsive to future climate changes by testing their sensitivity to past variation in their envi- ronment. These include responses to extreme events in the CCLME, like the exceptionally

strong El Niño event in 2015–2016 [46] and the marine heat wave dubbed ‘the blob’ that hit the Northeast Pacific in 2014–2015 [47]. Sensitivity of processes and communities to oceanic climate variation Along the West Coast of North America, spatial differences in average oceanic conditions are known to be important drivers of intertidal community structure [25,48–53]. For example, community structure clearly varies along this coastline, and this has been attributed to several environmental drivers including temperature (air and sea), nutrient concentrations, upwell- ing, upwelling variation, precipitation, and tidal range [48,50,53]. Clearly, the average oceanic climate can have a strong influence on intertidal communities in different places, but here we focus on whether interannual fluctuations in oceanic conditions relate to community structure at a given location over time. Overall, we are interested whether intertidal communities are sensitive to interannual oceanic climate fluctuations, or instead remain consistent despite such changes. Recent evidence suggests that multiple ecological processes at population and organismal levels are sensitive to oceanic variation along the Oregon California coasts. For example, monthly mussel recruitment tended to vary with NPGO and upwelling but not PDO or ENSO [36,54]. Barnacle recruitment also appears to be sensitive to upwelling [34,55,56], but not in all cases [57–59]. However, order-of-magnitude increases in mussel recruitment associated with

changing climate regime in the early 2000s [20] did not translate into changes in community structure [20,36]. Further, while mussel recruitment was sensitive to environmental condi- tions, mussel cover was not, suggesting a decoupling between recruitment and community structure dynamics [60]. Oceanic climate can also affect growth rates and physiology of intertidal species. Higher phytoplankton abundance and warmer water during warm-phase PDO and El Niño years cor- related with increased growth rates of mussels [61,62]. The kelp Egregia menziesii grew faster in a strongly upwelled region north of Point Conception, California than it did in the weakly upwelled region south of Point Conception [63]. Changes in seasonal and yearly upwelling, nutrients, water temperature, ENSO and PDO correlated to changes in growth rates and Car- bon:Nitrogen ratios in multiple dominant macrophyte species along the Oregon Coast [64,65]. The studies above also highlight how studies lasting just several years to one or two decades PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 3 / 34

can successfully detect the ecological consequences of climate oscillations operating on multi- decadal scales (e.g. PDO, NPGO). Increasing evidence suggests that rocky intertidal community structure is already respond- ing to climate change. Smith et al. [66] documented lower diversity in California rocky inter- tidal communities in 2002 compared to the 1960s and 1970s and suggested either a warm phase of PDO or warming global temperatures were responsible. Poleward range shifts of rocky shore species has occurred in northern Europe [67], and the relative dominance of northern versus southern species has shifted in California [68]. Globally, many coastal ecosys- tems have shown sometimes dramatic susceptibility to warming climate [e.g. 69–72]. Here, we examine if annual shifts in intertidal community structure in Oregon and northern California correlate to oceanic environmental variation to identify which, if any, environmental factors may drive community change. By identifying these drivers, we may be better equipped to pre- dict and perhaps mitigate the effects of changing climate on these iconic communities. Approach Apart from the simple hypotheses that rocky shore communities did not (null) or did (alterna- tive) respond to recent changes in ocean climate, our approach was primarily focused on dis- covering patterns. To determine if intertidal community structure changed during a period of intensifying climate change [35], we used a 14-year dataset (2006–2020) taken at 11 sites in Oregon and northern California spanning 294 km (183 miles) of coastline. We investigated

whether these communities were sensitive to yearly variation in PDO, NPGO, El Niño, upwell- ing intensity, upwelling variability, air temperature and water temperature during this time. We also included a time series of sea star wasting disease occurrence in our analysis because it caused coast-wide declines in the abundance of the keystone predator, the sea star Pisaster ochraceus [9], so may have caused strong ecological effects in this system [73,74]. Though wast- ing disease was more severe in warmer locations [75], there is either no or modest evidence that the outbreak was triggered by environmental factors [9,10,76–78]. Because the periodicity of some climate oscillations (e.g., PDO, NPGO) is typically longer (10–30 yr) than most eco- logical data sets, our ability to investigate the full extent of their importance is limited. How- ever, when multidecadal datasets are unavailable, examining the relationships between ecological patterns and yearly variation in these climate indices over several years can provide insight into their potential importance [20,54,62]. We quantified percent cover of sessile taxa and density of mobile taxa yearly in the high, mid and low intertidal. We asked 1) how have intertidal community structure, abundance of focal functional groups, and diversity varied over time? We then investigated sensitivity of these communities to climate variation by asking 2) have the abundance of focal functional groups, rate of community change, or diversity cor- related with yearly fluctuations in PDO, NPGO, El Niño, upwelling intensity, upwelling vari-

ability, air and water temperatures, and sea star wasting disease? Materials and methods Community surveys Starting in 2006, summer (June-August) community surveys were conducted annually at 11 sites on 4 capes (i.e., headlands or regions) from central Oregon to northern California (Fig 1, S1 Table in S1 Appendix). From north to south, Capes included Cape Foulweather (sites Fogarty Creek and Boiler Bay), Cape Perpetua (sites Yachats Beach, Strawberry Hill and Toka- tee Klootchman), Cape Blanco (sites Cape Blanco North and South, Port Orford Heads and Rocky Point), and Cape Mendocino (sites Cape Mendocino North and South). Scientific col- lection permits for this work were granted by the Oregon Department of Fish and Wildlife PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 4 / 34

(ODFW 19272, 21042, 21848) and the California Department of Fish and Wildlife (SC-4055, S-182930001-18295-001). With exceptions noted below, we surveyed communities separately in each of 3 intertidal zones, with the low zone defined as below mussel beds, the mid zone as mussel beds, and the high zone as above mussel beds. Because of time constraints or lack of habitat, we performed only low zone surveys but not mid or high surveys at Tokatee Klootch- man, Cape Blanco South, Port Orford Heads and both Cape Mendocino sites. In addition, we began most mid and high zone surveys in 2011, five years after the low zone surveys. Because the low zone is patchier and more diverse than the mid and high zones, we used three 30 m transect line(s) parallel to the shore in the low zone at least 20 m apart versus one 30 m transect line in mid and high zones. On each transect, we uniformly (every three meters along transect tapes in mid and high zones) or haphazardly (quadrat tosses in the low zone, avoiding sand, tide pools and large crevices) sampled 10 0.25m2 quadrats (total of 4,588 quadrats). Because we sampled the low zone more thoroughly than the mid and high zones, this may affect met- rics like biodiversity and community structure at the site level, so we analyzed zones separately in all except one model in this study (See Data Analyses). We performed surveys in the same general area each year, but the quadrat locations varied between years. In each gridded quad- rat, we visually estimated the percent cover to the nearest 1% of all visible species, including

sessile algae, sessile invertebrates and bare space and counted total numbers of all mobile Fig 1. Map of sites in Oregon and northern California. Capes included Cape Foulweather in greens (Boiler Bay and Fogarty Creek sites), Cape Perpetua in reds (Yachats Beach, Strawberry Hill and Tokatee Klootchman sites), Cape Blanco in blues (Cape Blanco North, Cape Blanco South, Port Orford Heads and Rocky Point sites), and Cape Mendocino in purples (Cape Mendocino North and South sites). See S1 Table in S1 Appendix for coordinates. https://doi.org/10.1371/journal.pone.0297697.g001 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 5 / 34

ENSO indices, we assigned the bimonthly designation (e.g. DecJan) to the latter month (e.g. January, [22]. We measured intertidal sea surface and air temperatures every 15 minutes at each site using temperature-logging sensors (HOBO Tidbits or HOBO Pendants, Onset Corp., Bourne, MA). We separated air from water temperatures using a de-tiding method coupled with tide tables (http://tidesandcurrents.noaa.gov) [62]. We then calculated the maximum daily air temperature and the mean daily water temperature at each site. For sea star wasting disease, we included the years before or since the SSWD outbreak at each site, with zero, nega- tive, and positive years demarcating the onset of, prior to, and after SSWD in 2013 at Califor- nia sites and 2014 at Oregon sites [9,84]. Since some of the indices have a seasonal component during which their effects are most prominent, we first determined the window of peak months as from Oct—Mar for PDO, Dec—Mar for NPGO, Dec—Apr for ENSO, March to the summertime survey month for upwelling, and April to the survey month for maximum air temperature [22,85–87]. We calculated the averages among months for each index during those peak months only, and then assigned this average to the subsequent annual community survey. Since upwelling intermittency can have strong effects on community processes [16,25], we estimated this metric using the standard deviation of the mean upwelling index from March to the survey month. Data analyses Our overarching goal was to investigate if intertidal community structure varied over time and

space, then test whether temporal variation was associated with interannual oceanic climate variation. We performed all analyses in Primer-e 7 with PERMANOVA+ add-on [88,89], JMP Pro v15 (SAS Institute), R v4.0.0 [90] and R Studio v1.2.5042 (R Studio Inc.). Below, we sum- marize our general approach for community analyses, then specific features of each model. Multivariate community structure analyses. General approach. For all statistical models, we first aggregated the species data to the functional group level (see S2 Table in S1 Appendix for list) then log10-transformed abundances to increase or decrease the influence of rare or common taxa, respectively. We quantified community similarity using Bray-Curtis PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 6 / 34

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

Table 1. Results of a) a PERMANOVA analysis testing the effects of zones, capes, sites, years, and their interactions on intertidal community dissimilarity in Oregon and Northern California from 2006–2021. Cape, site, and year were treated as random factors. b) Results of SIMPER analysis testing the contributions of taxa to within-zone similarity. Only the top ten taxa contributing most strongly to within-zone similarity are listed. Communities were analyzed at the transect level (the average of ~10 0.5 x 0.5m quadrats). a) PERMANOVA for Zone, Cape, Site, and Year Source df SS MS Pseudo-F P Estimate Sq.root % Explained Zone 2 281,830 140,920 18.43 0.001 1698.00 41.21 55.3 Cape 3 88,568 29,523 3.01 0.001 187.65 13.70 6.1 Year 14 53,610 3,829 4.10 0.001 97.61 9.88 3.2 Site[Cape] 7 53,955 7,708 11.77 0.001 190.69 13.81 6.2 Zone*Cape 4 24,899 6,225 2.28 0.001 146.05 12.09 4.8 Zone*Year 20 22,113 1,106 2.76 0.001 100.37 10.02 3.3 Cape*Year 38 31,580 831 1.26 0.008 23.00 4.80 0.7 Zone*Site[Cape] 6 14,146 2,358 6.31 0.001 164.28 12.82 5.3 Site[Cape]*Year 76 48,494 638 1.43 0.001 61.88 7.87 2.0 Zone*Cape*Year 35 13,968 399 1.07 0.299 10.66 3.27 0.3 Zone*Site[Cape]*Year 46 17,110 372 0.83 0.953 -53.33 -7.30 -1.7 Res 211 94,172 446 446.31 21.13 14.5 Total 462 790,780 3073.18 b) SIMPER For Within-Zone Similarity Taxon Av.Abund Av.Sim Sim/SD Contrib% Cum.% High Zone Bare Space 3.58 16.77 3.23 26.98 26.98 Balanus 2.57 10.11 2.08 16.27 43.26

Herbivorous snails 3.21 8.62 0.98 13.87 57.13 Limpets 2.62 8.49 1.44 13.66 70.78 Fucoid algae 1.4 4.09 1.05 6.59 77.37 Chthamalus 1.16 3.88 1.87 6.25 83.62 Mussels 1 2.85 1.22 4.58 88.2 Mastocarpus 0.67 1.93 1.12 3.1 91.3 Endocladia 0.52 1.27 0.85 2.04 93.35 Algal crusts 0.56 1.26 0.78 2.03 95.37 Mid Zone Mussels 4.05 26.13 2.74 42.71 42.71 Limpets 3.3 15.62 1.79 25.53 68.23 Predatory snails 1.35 4.59 1.07 7.5 75.73 Bare Space 1.13 4.21 1.33 6.88 82.61 Large barnacles 0.92 3.05 1.16 4.98 87.59 Gooseneck barnacles 0.79 2.31 0.79 3.78 91.37 Balanus 0.5 1.28 0.81 2.09 93.46 Chthamalus 0.31 0.8 0.76 1.3 94.76 Anemones 0.21 0.53 0.69 0.86 95.62 Endocladia 0.3 0.49 0.3 0.8 96.42 Low Zone Bare Space 2.83 9.87 3.29 18.14 18.14 Articulated coralline algae 2.27 7.38 2.48 13.56 31.7 Coralline algae crusts 1.56 4.21 1.59 7.73 39.43 Limpets 1.66 3.72 1.12 6.84 46.27 Hedophyllum 1.36 2.98 0.97 5.47 51.74 Fine branching reds 1.09 2.78 1.37 5.1 56.84 Mazzaella blades 1.08 2.68 1.28 4.92 61.76 Surf grass 1.05 2.08 0.8 3.82 65.58 Anemones 0.87 1.9 0.94 3.5 69.08 Algal crusts 0.83 1.9 1.14 3.49 72.57 https://doi.org/10.1371/journal.pone.0297697.t001 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 8 / 34

significant climate index terms, we modeled and graphed the bivariate correlation between the cover or density of the taxon and the climate variable at the transect level (using the lm() func- tion in the stats package and the ggplot2 package [95]). Temporal trajectories of community structure. Because strong effects of site and cape made temporal patterns difficult to distinguish in nMDS space, we calculated centroids of Bray-Cur- tis similarity for each site and year, then plotted year to year trajectories using nMDS and saved the coordinates. We used the vector lengths of these trajectories (quantified using as.ltraj function in the adehabitatLT package in R [96]) to calculate annual rates of community shifts [11,97]. For each zone, we then tested if the log-transformed rates of community shifts differed 1) among years, with cape and sites (nested within cape), as random variables, and 2) with the 8 climate variables, with site as a random variable (lmer function in lme4 package [94]). Taxon diversity. Diversity analyses used the Shannon diversity index (H’ calculated using natural log of species or taxon cover to the highest possible resolution) for each transect. We tested whether taxon diversity varied among sites, years and their interaction using 3 linear models (one per zone, lm() function in stats package [90]). Next, we analyzed whether diver- sity in each zone varied with site (as a random factor) and the 8 climate indices using linear mixed models (lmer function in the lme4 package [94]). For significant climate variable terms

(Bonferroni cutoff at P = 0.002), we regressed the focal climate indices on taxon diversity. Results Environmental trends Temporal patterns of climate indices varied considerably during the study period (Fig 2). For example, the PDO mostly was negative (cooler conditions) from 2006 to 2013, then positive (warmer) in 2015 to 2017 while gradually declining then moving to negative in 2018 to 2020 (Fig 2a). The NPGO shifted from negative (calmer conditions) in 2006 to 2007 to positive (windier) from 2008 to 2013, then switched to negative and generally declined from 2014 to 2020 (Fig 2b). Finally, the ENSO index oscillated between cooler La Niña (2008, 2009, 2011) and warmer El Niño (2010, 2016) conditions (Fig 2c). Upwelling mostly varied latitudinally, as expected, with mean upwelling being similar at the two northern capes (Capes Foulweather and Perpetua), stronger in southern Oregon (Cape Blanco), and strongest in northern Califor- nia (Cape Mendocino). It tended to decline from 2008 to 2012, spiked in 2013, then varied lit- tle annually after that (Fig 2d). Fig 2. Average yearly environmental indices (Mean ± SD) for a) Pacific Decadal Oscillation (PDO), b) North Pacific Gyre Oscillation (NPGO), c) Multivariate El Niño Southern Oscillation Index (ENSO) and d) Biologically Effective Upwelling Transport Index, colored by Cape (BEUTI). All yearly means were calculated from daily or monthly indices during seasonal peak months for each variable, with PDO index from October—March, NPGO index from December—March, ENSO from December—April, and

BEUTI from March through the survey month. These data were included in models analyzing community structure. https://doi.org/10.1371/journal.pone.0297697.g002 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 9 / 34

For temperature trends, both mean maximum air temperature (Fig 3a) and mean water temperature (Fig 3b) generally increased during the study period. The strongest increases occurred for maximum daily air temperature at Rocky Point, Port Orford Heads, Cape Blanco North, and Yachats Beach (P < 0.01, with R2 = 0.22. 0.12, 0.07, and 0.09, respectively). As reported earlier [12], a period of elevated water temperature occurred during the confluence of the 2014–16 marine heat wave and the 2015–16 El Niño (Fig 3b). As expected, warmer water temperature was correlated to stronger PDO (R2 = 0.523, P = < 0.001), but no other correla- tions of temperatures with climate oscillations were significant. Mean daily water temperatures were similar between capes and between sites (Fig 3b). Maximum daily air temperatures were generally warmer at Cape Blanco that the other capes, and air temperatures were fairly similar among sites within capes (Fig 3a). The exception was Cape Blanco, where maximum daily air temperatures were typically lower at Port Orford Heads and higher at Rocky Point and Cape Blanco South compared to the Cape Blanco North site. Spatial versus temporal community variation Community structure was more strongly associated with spatial (zone, cape and site) than temporal variables (year) (Table 1). As expected, zones were clearly separated visually (Fig 4), as were capes and often sites (Fig 5). However, visual separation among years was much less clear (Fig 6a–6c). Within each zone, cape and site were associated with more community vari-

ation than was year (S3a-S3c Table in S1 Appendix; High: 35.6% vs. 19.9%, Mid: 33.2% vs. 19.4%, and Low: 48.4% vs. 14.7%, for spatial vs. temporal variation, see calculation in methods). Community structure among intertidal zones As indicated above, community structure was strongly associated with intertidal zone, which contributed 55.3% of model fit when analyzing the contributions of zone, cape, site and year (Fig 4, Table 1; PERMANOVA: Zone: p = 0.001). In SIMPER analyses, the strongest drivers of within-zone community similarity (i.e., contributing >5%) for the high zone were bare space, barnacles, herbivorous gastropods, and fucoid algae (Fig 4; S3d Table in S1 Appendix; S1 Fig in S3 Appendix). In the mid zone, mussels, limpets, whelks, bare space, and large barnacles drove community similarity (Fig 4; S3e Table in S1 Appendix; S2 Fig in S3 Appendix). The low zone was typified by bare space, articulated coralline algae, crustose coralline algae, limpets, the kelp Hedophyllum sessile, and fine branching red algae (Fig 4; S3f Table in S1 Appendix; S3 Fig in S3 Appendix). Geographic variation in community structure Taxon abundance often differed by cape and site (Fig 5). Although community structure did not differ by cape in the high zone (Fig 5a; S3a Table in S1 Appendix) and functional group rel- ative abundance was similar among capes (S3d Table in S1 Appendix), community structure differed among sites within capes (Fig 5a; p < 0.007 for all site comparisons). These trends sug- gest that high zone variation was controlled at a localized scale.

In the mid zone, among-cape differences (Fig 5b; S3b Table in S1 Appendix) were driven by separation of Cape Blanco from Capes Foulweather and Perpetua (p = 0.023, p = 0.014, respectively), which did not differ from each other (p = 0.311). As in the high zone, sites dif- fered within capes (Fig 5b; S3b Table in S1 Appendix; p < 0.018 for site comparisons). While dominant mid-zone species were often similar among capes, the strongest contributors to community separation among capes were mussels and limpets at Cape Blanco (Fig 5b; S3e Table in S1 Appendix; 58.9% and 31.2% contribution, respectively), predatory snails and bare PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 10 / 34

Fig 3. Average yearly measures (Mean ± SD) among capes (panels) and sites (colors) for a) the maximum daily air temperature (Max Air Temp), and b) the mean daily water temperature (Mean SST). Cape Foulweather (greens) includes Boiler Bay (BB) and Fogarty Creek (FC), Cape Perpetua (reds) includes Yachats Beach (YB), Strawberry Hill (SH) and Tokatee Klootchman (TK), Cape Blanco (blues) includes Cape Blanco North(CBN), Cape Blanco South (CBS), Port Orford Heads (POH) and Rocky Point (RP), and Cape Mendocino (purples) includes Cape Mendocino North (CMEN) and South (CMES). Dashed lines are linear fits of the data. All daily values were calculated from loggers recording every 15 minutes, then the yearly means and standard deviations were calculated from daily values. Air temperature values were taken from warm seasons (from April to the survey month) and water temperature was taken for the year. Missing data are from years where we were unable to maintain or discontinued the instrument deployments. These data were included as covariates in models analyzing community structure. https://doi.org/10.1371/journal.pone.0297697.g003 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 11 / 34

space at Cape Perpetua (17.0% and 6.3% contribution, respectively) and gooseneck barnacles at Cape Foulweather (5.9% contribution). In the low zone, most capes differed from one another (Fig 5c; S3c Table in S1 Appendix) except Cape Blanco did not differ from Capes Foul- weather or Mendocino (p = 0.091 and 0.089 respectively). Sites again differed within capes (p < 0.032 for site comparisons). Like the other zones, in the low zone many of the same taxa dominated every cape, but their relative contributions to cape separation differed (Fig 5c; S3f Table in S1 Appendix). For example, Cape Foulweather was associated with relatively low bare space but with high coralline algae crusts, Hedophyllum sessile, and Cryptopleura ruprechtiana/ Hymenena flabelligera (Fig 5c; S3f Table in S1 Appendix, 11.4%, 10.9%, 7.7%, and 6.1% contri- bution, respectively) while Cape Perpetua had relatively more sea anemones, predatory sea stars, and gooseneck barnacles (10.0%, 5.1% and 5.0% contribution, respectively). Community structure over time Interannual differences occurred in all zones, with year accounting for 19.9%, 19.4%, and 14.7% of the community variation in high, mid and low zones, respectively (S3a-S3c Table in S1 Appendix; Year: p = 0.001 for each zone, see % calculation in methods). However, these Fig 4. Non-metric multidimensional scaling (nMDS) plots depicting Bray-Curtis dissimilarity in community structure among high (yellow), mid (green), and low (blue), and intertidal zones. Vector overlays depict species contributing >3% to within-zone similarity in SIMPER analyses (S2b

Table in S1 Appendix). Each data point represents a community at the transect level (the average of ~ten 0.5 x 0.5m quadrats). https://doi.org/10.1371/journal.pone.0297697.g004 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 12 / 34

Fig 5. Non-metric multidimensional scaling (nMDS) plots depicting Bray-Curtis dissimilarity among intertidal communities in a) high, b) mid, and c) low intertidal community structure in Oregon and northern California from 2006–2021. Colors depict Cape Foulweather sites in blues, Cape Perpetua sites in reds, Cape Blanco sites in greens, and Cape Mendocino sites in purples. Each data point represents a community at either the quadrat level (a & b, 0.5 x 0.5m) or transect level (c, the average of ~10 0.5 x 0.5m quadrats). Vector overlays (in black) depict species contributing most strongly to within-cape similarity in SIMPER analyses (S4d-S4f Table in S1 Appendix). https://doi.org/10.1371/journal.pone.0297697.g005 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 13 / 34

Fig 6. Non-metric multidimensional scaling (nMDS) plots depicting Bray-Curtis dissimilarity among intertidal communities in high (a & d), mid (b & e), and low (c & f) intertidal community structure in Oregon and northern California from 2006–2021. Left panels (a-c) depict dissimilarities among years, with each data point representing a community at either the quadrat level (a & b, 0.5 x 0.5m) or transect level (c, the average of ~10 0.5 x 0.5m quadrats). Left panel (a-c) vector overlays (in black) depict species contributing most strongly to within-cape similarity in SIMPER analyses (S4d-S4f Table in S1 Appendix). Right panels (d-e) depict community trajectories over time with each data point representing the centroid of community structure for each site and year (i.e. centroids from left panels). We calculated the length and absolute angle of these trajectories to estimate rate and direction of community shifts (respectively) for subsequent analyses. Right panel (d-e) vector overlays depict the association between the environmental variables and community dissimilarity. Right panel (d-e) colors depict Cape Foulweather sites in blues, Cape Perpetua sites in reds, Cape Blanco sites in greens, and Cape Mendocino sites in purples. https://doi.org/10.1371/journal.pone.0297697.g006 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 14 / 34

differences were not readily apparent in nMDS visualizations (Fig 6a–6c) and the temporal patterns varied with site in all zones and with cape in the low zone (S3a-S3c Table in S1 Appen- dix). In the high zone, only 2007 differed from each of the other years (p < 0.005), but this was because only Boiler Bay was sampled that year. In the mid zone, only 2012 differed from 2017 and 2018 (p = 0.038 and 0.044, respectively). On the other hand, in the low zone several years differed from one another, particularly 2006, 2008, and 2019 compared to other years. Community associations with oceanic climate Community and functional group patterns. In all zones, all oceanic climate fluctuations were related to community structure (p < 0.03 in most cases), but with the exception of years since sea star wasting disease in the mid zone (15.9% of variance explained), the amount of variance explained by each was low (ranging from 0.6 to 8.6%; Table 2). In the high zone, no single factor accounted for more than 3.8% of the variation in community structure, and together the 8 factors drove only 16.1% of the variation in community structure (Table 2a). In the mid zone, environment-community associations were stronger (Table 2b, p < 0.012). Col- lectively, the climate indices, upwelling and temperature accounted for 25.5% of the variance, with NPGO and PDO contributing the most strongly (8.6% and 6.1% of variation, respec- tively). As noted above, the years since sea star wasting accounted for a further 15.9% of the

variance explained, bringing the total variance explained by environment to 41.3% in the mid zone. In the low zone, all factors except upwelling variation and maximum daily air tempera- ture varied with community structure, but individual effects were weak with sea star wasting disease, NPGO and mean upwelling being the strongest drivers (Table 2c, 6.2%, 4.5%, and 4.2%, respectively). Collectively, environment explained only 20.1% of the low zone variance. Analysis of correlations between the most temporally variable functional groups (S4, S5 and S6 Figs in S3 Appendix) and the 8 oceanic environmental variables showed that several groups appeared sensitive to the environment in the high and mid zones, but fewer were sensitive in the low zone. In the high zone, both herbivorous snails and limpets were sensitive to several environmental variables (Fig 7; Table 3a), and the strongest trends were increases with windier phases of the NPGO (R2 = 0.157 and 0.122, respectively), decreases with years since SSWD (R2 = 0.239 and 0.151, respectively), and limpet abundance increased with higher upwelling vari- ability (R2 = 0.164). Finally, fucoid algal abundance in the high zone was weakly negatively cor- related with upwelling variability and mussel abundance was weakly positively correlated with PDO (Fig 7; Table 3a; R2 = 0.077 and < 0.001, respectively). Mid zone limpet abundance varied similarly to that in the high zone, increasing with windier NPGO and increased upwelling variation (R2 = 0.107 and 0.109, respectively),

decreasing with years since SSWD (R2 = 0.131), and additionally decreasing in warmer phases of PDO (R2 = 0.222, Fig 8, Table 3b). Mid-zone predatory snail abundance was lower when maximum air temperature was warmer (R2 = 0.117), mussel abundance increased with years since SSWD (R2 = 0.136), and the alga Endocladia muricata was more abundant with increased maximum air temperature (R2 = 0.255, Fig 8, Table 3b). In the low zone, individual functional groups were generally less sensitive to varying environments (Fig 9, Table 3c), with the only notable trends being decreased limpets with warmer sea surface temperatures (R2 = 0.123) and decreased surfgrass (Phyllospadix spp.) with increased upwelling variability (R2 = 0.111). Temporal trajectories of community structure. Trajectories of annual community shifts showed that communities did vary over time, but no consistent trend emerged within any zone or among years (Fig 6d–6f). In the high zone, the rate of community change was not cor- related with any of the 8 environmental variables (Table 4a). In the mid zone, community PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 15 / 34

change rates were weakly associated with mean ENSO index (R2 = -0.003; Fig 10a; Table 4b) and the years since sea star wasting disease (R2 = -0.019; Fig 10b; Table 4b). In the low zone, community change rate increased slightly with increasing maximum daily air temperature (R2 = 0.069; Fig 10c; Table 4c). Table 2. Results of 3 PERMANOVA analyses testing drivers of community dissimilarity in a) high, b) mid, and c) low intertidal zones in Oregon and northern Cali- fornia from 2006–2021. We tested the effects of site (nested within cape) as a random factor and the effects of climate variables for mean Pacific Decadal Oscillation index (mean PDO), mean North Pacific Gyre Oscillation index (mean NPGO), mean El Niño Southern Oscillation index (mean ENSO), mean and standard deviation of the Biologically Effective Upwelling Transport Index (mean and SD BEUTI, respectively), mean daily sea surface temperature (Mean SST), mean daily maximum air tem- perature (Max Air Temp), and years since the sea star wasting disease epidemic (Yr Since SSWD). Communities were analyzed at the quadrat level in the mid and high zones (0.5 x 0.5m) and the transect level in the low zone (the average of ~10 0.5 x 0.5m quadrats). Bold statistics indicate that >10% of variation was explained by a given factor. a) High Zone Source df SS MS Pseudo-F P Estimate Sq.root % Explained Mean PDO 1 7,842 7,842 8.92 0.001 55.0 7.41 3.8 Mean NPGO 1 7,583 7,583 8.62 0.001 55.1 7.42 3.8 Mean ENSO 1 6,526 6,526 7.42 0.001 32.5 5.70 2.2

Mean BEUTI 1 5,020 5,020 5.71 0.001 31.2 5.58 2.2 SD BEUTI 1 1,900 1,900 2.16 0.056 4.6 2.14 0.3 Mean SST 1 3,570 3,570 4.06 0.001 23.2 4.81 1.6 Max Air Temp 1 2,312 2,312 2.63 0.018 8.9 2.98 0.6 Yr Since SSWD 1 3,124 3,124 3.55 0.004 22.1 4.70 1.5 Site 5 121,360 24,272 27.60 0.001 333.6 18.27 23.1 Residual 560 492,460 879 879.4 29.65 60.8 Total 573 742,780 1445.4 b) Mid Zone Source df SS MS Pseudo-F P Estimate Sq.root % Explained Mean PDO 1 12,885 12,885 15.40 0.001 112.0 10.59 6.1 Mean NPGO 1 6,593 6,593 7.88 0.001 157.4 12.54 8.6 Mean ENSO 1 5,797 5,797 6.93 0.001 71.7 8.47 3.9 Mean BEUTI 1 5,090 5,090 6.08 0.001 43.9 6.63 2.4 SD BEUTI 1 6,355 6,355 7.60 0.001 35.1 5.93 1.9 Mean SST 1 2,742 2,742 3.28 0.012 28.6 5.35 1.6 Max Air Temp 1 3,296 3,296 3.94 0.002 16.3 4.04 0.9 Yr Since SSWD 1 8,274 8,274 9.89 0.001 292.2 17.09 15.9 Site 5 84,260 16,852 20.14 0.001 239.4 15.47 13.1 Residual 541 452,660 837 836.7 28.93 45.6 Total 554 663,550 1833.4 c) Low Zone Source df SS MS Pseudo-F P Estimate Sq.root % Explained Mean PDO 1 1,187 1,187 2.03 0.032 7.1 2.66 0.6 Mean NPGO 1 8,489 8,489 14.52 0.001 56.4 7.51 4.5 Mean ENSO 1 4,851 4,851 8.3 0.001 32.9 5.74 2.6 Mean BEUTI 1 3,065 3,065 5.24 0.001 53.2 7.3 4.2 SD BEUTI 1 587 587 1 0.43 0 0.13 0 Mean SST 1 2,579 2,579 4.41 0.001 23 4.79 1.8 Max Air Temp 1 730 730 1.25 0.252 1.5 1.21 0.1

Yr Since SSWD 1 8,977 8,977 15.35 0.001 78.6 8.87 6.2 Site 10 112,240 11,224 19.19 0.001 422.1 20.55 33.5 Residual 317 185,370 585 584.8 24.18 46.4 Total 335 365,900 1259.6 https://doi.org/10.1371/journal.pone.0297697.t002 PLOS ONE Resistance of rocky intertidal communities PLOS ONE | https://doi.org/10.1371/journal.pone.0297697 May 29, 2024 16 / 34

Taxon diversity Taxon diversity (H’) varied among sites, years, and their interaction in all three zones but most changes were small, accounting for little variance (S7 Fig in S3 Appendix; S4 Table in S1 Appendix). For example high zone diversity increased slightly with warmer PDO and cooler maximum air temperatures (R2 = 0.042 and 0.029, respectively; Fig 11a&11b; Table 5a). Mid zone diversity increased with windier NPGO and La Niña conditions and decreased with years Fig 7. Bivariate regressions between the percent cover of focal functional groups (significant functional groups in Table 2a with the 8 focal climate variables in the high zone. Each data point is the average abundance of a functional group for a site and year. Environmental drivers included Pacific Decadal Oscillation index (mean PDO), North Pacific Gyre Oscillation index (mean NPGO), Multivariate El Niño Southern Oscillation Index (mean ENSO), Biologically Effective Upwelling Transport Index (mean BEUTI), the standard deviation in BEUTI (SD BEUTI), the maximum daily air temperature in celcius by site (Max Air Temp), the mean daily water temperature in celcius by site (Mean SST) and the years since the initial sea star wasting disease outbreak (Yrs since SSWD). Colors correspond to functional groups, with fucoid algae in green, mussels (Mytilus spp.) in in blue, herbivorous snails in pink, and limpets in purple. Units for fucoid algae and mussels are in percent cover while limpets and herbivorous snails are in density m-2. Some high