Resistance of rocky intertidal communities to oceanic climate fluctuations

macroinvertebrates. Due to layering of macrophytes or the epibionts growing on mussels, cover could total >100% when summing all species surveyed. We identified biota to species whenever possible, but those that we were not able to identify reliably to species or were very rare were lumped at the genus level or higher for diversity metrics (see S2 Table in S1 Appen- dix for species list). Oceanic climate variables Climate indices used here oscillate on a 20–30 year scale for the Pacific Decadal Oscillation (PDO), 7–15 years for the North Pacific Gyre Oscillation (NPGO), 3–7 years for El Niño Southern Oscillation (ENSO), and inter- and intra-annually for upwelling. Please see S2 Appendix for a brief background on these oceanic climate variables and how each may be affected by global climate change. We downloaded data for each of these oscillations for the northern California Current System between 2006 and 2020, including the PDO from the National Centers for Environmental Information at NOAA [13,79,80], the NPGO index from the NPGO portal [19], the ENSO from NOAAs Physical Sciences laboratory using Multivariate ENSO Index Version 2 (MEI.v2) [22,81], and upwelling as the monthly means of daily indices of the Biologically Effective Upwelling Transport Index BEUTI at 1˚ resolution, [82,83]. PDO, NPGO and ENSO are ocean basin-wide measurements [13,19,22], so the values were the same for all sites each year. For upwelling, we used the latitude appropriate for each cape (44˚N for Capes Foulweather and Perpetua, 42˚N for Cape Blanco and 40˚N for Cape Mendocino). For

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