Resistance of rocky intertidal communities to oceanic climate fluctuations

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