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 temperatures 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 variation 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 oceanic climate fluctuations that occurred during the study. However, given ongoing intensification of climate change and increasing frequencies of extreme events, future responses to climate change seem likely.

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, Appendix S1: Table S1). From north to south, Capes included Cape Foulweather (sites Fogarty Creek and Boiler Bay), Cape Perpetua (sites Yachats Beach, Strawberry Hill and Tokatee 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 collection permits for this work were granted by the Oregon Department of Fish and Wildlife (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 Klootchman, 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.25m² quadrats (total of 4,588 quadrats). Because we sampled the low zone more thoroughly than the mid and high zones, this may affect metrics 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 quadrat, 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 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 Appendix S1: Table S2 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 Appendix S3 for 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), and upwelling as the monthly means of daily indices of the Biologically Effective Upwelling Transport Index BEUTI at 1° resolution. 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). We measured intertidal sea surface and air temperatures every 15 minutes at each site using temperature-logging sensors (HOBO Tidbits or HOBO Pendants, Onset Coarp., Bourne, MA). We separated air from water temperatures using a de-tiding method coupled with tide tables (http://tidesandcurrents.noaa.gov). 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, negative, and positive years demarcating the onset of, prior to, and after SSWD in 2013 at California sites and 2014 at Oregon sites. 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. 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, and then test whether temporal variation was associated with interannual oceanic climate variation. We performed all analyses in Primer-e 7 with PERMANOVA+ add-on, JMP Pro v15 (SAS Institute), R v4.0.0 (R Core Team) and R Studio v1.2.5042 (R Studio Inc.). Below, we summarize 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 Appendix S1: Table S2 for list) then log₁₀-transformed abundances to increase or decrease the influence of rare or common taxa, respectively. We quantified community similarity using Bray-Curtis 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). 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. We used similarity percentages analysis (SIMPER) to identify the functional groups substantially contributing to dissimilarity among or similarity within levels of a factor. 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 (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 analyses separately by zone.

Community associations with space and time. To investigate the relative importance of temporal 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 prevented 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 initial percent contribution calculation for each term) of the spatial and temporal components of variation were estimated as: spatial % contribution = Σ(cape + site + 0.5*(year x cape) + 0.5*(year x site), and temporal % contribution = Σ(year + 0.5*(year x cape) + 0.5*(year x site].

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 surveys). Because environmental data were on different scales, we normalized them (ranging from -2 to +2) before analysis using Primer 7 with the PERMANOVA+ add-on. We inspected draftsman plots and found no substantial collinearity among the 8 environmental variables that might interfere with models (All R² values < 0.4, except R² = 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 functional groups by zone (highest average dissimilarity among years). We then analyzed their sensitivity to annual environmental fluctuations with univariate generalized linear mixed models (GLMMs) using the glmer function in the lme4 package in R. Site was random and climate indices were fixed. We used a Bonferroni correction threshold of p <0.00625 (8 terms). For 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() function in the stats package and the ggplot2 package.

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-Curtis 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) to calculate annual rates of community shifts. 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).

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 diversity 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.

Microsoft Excel

R and R Studio

PRIMER-e with PERMANOVA+ Add-on