Resistance of rocky intertidal communities to oceanic climate fluctuations
Data files
Jan 10, 2024 version files 3.46 MB
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 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.
README: Resistance of rocky intertidal communities to oceanic climate fluctuations
Data for Gravem et al.
This readme file was generated on 2024-01-08 by Sarah Gravem
Description of the data and file structure
Date of data collection: 2006-2020
Geographic location of data collection: Oregon and Northern California Coasts
Information about funding sources that supported the collection of the data: The National Science Foundation, The Packard Foundation, The Kingfisher Foundation
There are four type of data associated with this manuscript.
1) EnvFocalNorm.xlsx includes the peak seasonal mean environmental data. These are assigned to each quadrat for use in PRIMER-e, but vary by year and site. These include 8 variables (PDO, NPGO, ENSO, Mean Upwelling, SD Upwelling, Max Daily Air Temperature, Mean Daily Sea Surface Temperature, and years since sea star wasting disease). Data are normalized. "Call" column matches these data to a given community survey quadrat found in FunctionalGroupAbundance.xlsx and Diversity.xlsxfor a given site and year
2) FunctionalGroupAbund.xlsx shows taxon abundances by quadrat (0.25m2) of taxonomic functional groups surveyed yearly at each site. In percent cover for all sessile species and in densities for all mobile species. "Call" column matches these data to a given community survey quadrat found in EnvFocalNorm.xlsx and Diversity.xlsx for a given site and year
3) Diversity.xlsx is derived from FunctionalGroupAbund.xlsx and measures various diversity indices calculated by PRIMER-e for each quadrat surveyed. "Call" column matches these data to a given community survey quadrat found in EnvFocalNorm.xlsx and FunctionalGroupAbund.xlsx for a given site and year
4) Three CommTrajectory.xlsx files. Separate files for three different intertidal zones (high, mid low). Trajectories of community change in nMDS space calculated in PRIMER-e and analyzed using the adehabitatLT package in R. dist is akin to change in community structure in a given year (yearly vector length in nMDS space), and abs.angle is akin to direction of community change in a given year (yearly vector angle in nMDS space)
DATA-SPECIFIC INFORMATION FOR: FunctionalGroupAbund.xlsx
Number of variables: 83
Number of cases/rows: 4587
n/a indicates metric not available
Variable List:
- Call: Unique Identifier for each Quadrat includes Site_Year_Exposure_Zone_Area_Quadrat
- Bare Space through Ulva: Percent Cover of a given species
- Predatory sea stars through Nudibranchs: Density of a given species per m2
- cape: Promontory/cape/point the site was located on
- cape number: Order of cape north to south
- site code: Code abbreviation for the site
- site number: Order of site north to south
- Site: Name of site
- lat: Latitude of site in decimal degree
- RndLat: Rounded latitude (to match BEUTI index)
- lon: Longitude of site in decimal degrees
- State: US State
- Date: Date of survey
- year: Year of survey
- Area: Area of survey within site
- Zone: Intertidal zone of survey as low, mid or high
- ExposureCode: Exposure of site as protected or exposed
- Exposure: Exposure of site as protected or expo
- Quadrat: Quadrat number. All were 0.25m2
- PDO_MeanSeas: Pacific Decadal Oscillation Index for Oct – Mar the year survey taken
- PDO_SDSeas: Std Deviation of the Pacific Decadal Oscillation Index for Oct – Mar the year survey taken
- NPGO_MeanSeas: North Pacific Gyre Oscillation Index for Dec – Mar the year survey taken
- NPGO_SDSeas: Std Deviation of the North Pacific Gyre Oscillation Index for Dec – Mar the year survey taken
- MEI_MeanSeas: Multivariate El Nino Souther Oscillation index for Dec – Apr the year survey taken
- MEI_SDSeas: Std Deviation of the Multivariate El Nino Souther Oscillation index for Dec – Apr the year survey taken
- BEUTI_MeanSeas: Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- BEUTI_SDSeas: Std Deviation of the Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- CUTI_MeanSeas: Coastal Upwelling Transport Index for March to survey month the year survey taken
- CUTI_SDSeas: Std Deviation of the Coastal Upwelling Transport Index for March to survey month the year survey taken
- MeanDailyWaterTemp_AllYear: Average daily water temperature for the year the survey was taken
- MeanMaxDailyPeakAirTemp_PeakSeason: Average daily maximum air temperature from April to the survey month for the year the survey was taken
- SDDailyWaterTemp_AllYear: Std Deviation of the Average daily water temperature for the year the survey was taken
- SDMaxDailyPeakAirTemp_PeakSeason: Std Deviation of the Average daily maximum air temperature from April to the survey month for the year the survey was taken
- YearSinceSSWD: Year since the outbreak of sea star wasting disease
- AreaUnique: Unique identifier for the area concatenated from SiteCodeZoneExposureCodeAreaYear
- Siteyear: Site and Year
- capeyear: Cape and Year
DATA-SPECIFIC INFORMATION FOR: EnvFocalNorm.xlsx
Number of variables: 26
Number of cases/rows: 5167
n/a indicates metric not available
Variable List:
- Call: Unique Identifier for each Quadrat includes Site_Year_Exposure_Zone_Area_Quadrat
- PDO_MeanSeas: Pacific Decadal Oscillation Index for Oct – Mar the year survey taken
- PDO_SDSeas: Std Deviation of the Pacific Decadal Oscillation Index for Oct – Mar the year survey taken
- NPGO_MeanSeas: North Pacific Gyre Oscillation Index for Dec – Mar the year survey taken
- NPGO_SDSeas: Std Deviation of the North Pacific Gyre Oscillation Index for Dec – Mar the year survey taken
- MEI_MeanSeas: Multivariate El Nino Souther Oscillation index for Dec – Apr the year survey taken
- MEI_SDSeas: Std Deviation of the Multivariate El Nino Souther Oscillation index for Dec – Apr the year survey taken
- BEUTI_MeanSeas: Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- BEUTI_SDSeas: Std Deviation of the Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- CUTI_MeanSeas: Coastal Upwelling Transport Index for March to survey month the year survey taken
- CUTI_SDSeas: Std Deviation of the Coastal Upwelling Transport Index for March to survey month the year survey taken
- MeanDailyWaterTemp_AllYear: Average daily water temperature for the year the survey was taken in Celsius
- MeanMaxDailyPeakAirTemp_PeakSeason: Average daily maximum air temperature from April to the survey month for the year the survey was taken
- SDDailyWaterTemp_AllYear: Std Deviation of the Average daily water temperature for the year the survey was taken
- SDMaxDailyPeakAirTemp_PeakSeason: Std Deviation of the Average daily maximum air temperature from April to the survey month for the year the survey was taken
- YearSinceSSWD: Year since the outbreak of sea star wasting disease
- cape: Promontory/cape/point the site was located on
- cape number: Order of cape north to south
- site code: Code abbreviation for the site
- site number: Order of site north to south
- Site: Name of site
- State: US State
- Date: Date of survey
- year: Year of survey
- Area: Area of survey within site
- Zone: Intertidal zone of survey as low, mid or high
- ExposureCode: Exposure of site as protected or exposed
- Exposure: Exposure of site as protected or expo
- Quadrat: Quadrat number. All were 0.25m2
- EpidemicPhase: Before or after the outbreak of sea star wasting disease at the site
- AreaUnique: Unique identifier for the area concatenated from SiteCodeZoneExposureCodeAreaYear
- Siteyear: Site and Year
- capeyear: Cape and Year
DATA-SPECIFIC INFORMATION FOR: CommTrajectory_High.xlsx, CommTrajectory_Mid.xlsx, CommTrajectory_Low.xlsx
Number of variables: 27, 27, 27, respectively
Number of cases/rows: 57, 55, 139, respectively
n/a indicates metric not applicable
Variable List: (same for all 3 files)
- x: average x coordinate of community structure for a given site and year in nDMS space
- y: average y coordinate of community structure for a given site and year in nDMS space
- date: Date and time of survey as posixct
- dx: Change of mean community coordinate in horizontal nMDS distance between years
- dy: Change of mean community coordinate in vertical nMDS distance between years
- dist: Change of mean community coordinate in euclidian distance in nMDS distance between years
- dt: Time elapsed as posixct
- R2n: Squared diet
- abs.angle: Absolute Ange of trajectory in radians
- rel.angle: Relative Ange of trajectory in radians
- id: site identifier
- burst: site identifier
- pkey: site and date identifier
- year: Year of survey
- site code: Code abbreviation for the site
- Cape: Promontory/cape/point the site was located on
- Site: Name of site
- Latitude: Latitude of site in decimal degrees
- Longitude: Longitude of site in decimal degrees
- PDO_MeanSeas: Pacific Decadal Oscillation Index for Oct – Mar the year survey taken
- NPGO_MeanSeas: North Pacific Gyre Oscillation Index for Dec – Mar the year survey taken
- MEI_MeanSeas: Multivariate El Nino Souther Oscillation index for Dec – Apr the year survey taken
- BEUTI_MeanSeas: Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- BEUTI_SDSeas: Std Deviation of the Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- MeanDailyWaterTemp_AllYear: Average daily water temperature for the year the survey was taken in Celsius
- MeanMaxDailyPeakAirTemp_PeakSeason: Average daily maximum air temperature from April to the survey month for the year the survey was taken in Celsius
- YearSinceSSWD: Year since the outbreak of sea star wasting disease
DATA-SPECIFIC INFORMATION FOR: Diversity.xlsx
Number of variables: 42
Number of cases/rows: 4587
n/a indicates metric not available
Variable List:
- Call: Unique Identifier for each Quadrat includes Site_Year_Exposure_Zone_Area_Quadrat
- S: Total species
- N: Total individuals
- d: Species richness (Margalef)
- J': Pielou’s evenness
- H'(loge): Shannon-Weiner Diversity Index
- 1-Lambda': Simpson Diversity Index
- cape: Promontory/cape/point the site was located on
- cape number: Order of cape north to south
- site code: Code abbreviation for the site
- site number: Order of site north to south
- Site: Name of site
- lat: Latitude of site in decimal degrees
- RndLat: Rounded latitude (to match BEUTI index)
- lon: Longitude of site in decimal degrees
- State: US State
- Date: Date of survey
- year: Year of survey
- Area: Area of survey within site
- Zone: Intertidal zone of survey as low, mid or high
- ExposureCode: Exposure of site as protected or exposed
- Exposure: Exposure of site as protected or expo
- Quadrat: Quadrat number. All were 0.25m2
- PDO_MeanSeas: Pacific Decadal Oscillation Index for Oct – Mar the year survey taken
- PDO_SDSeas: Std Deviation of the Pacific Decadal Oscillation Index for Oct – Mar the year survey taken
- NPGO_MeanSeas: North Pacific Gyre Oscillation Index for Dec – Mar the year survey taken
- NPGO_SDSeas: Std Deviation of the North Pacific Gyre Oscillation Index for Dec – Mar the year survey taken
- MEI_MeanSeas: Multivariate El Nino Souther Oscillation index for Dec – Apr the year survey taken
- MEI_SDSeas: Std Deviation of the Multivariate El Nino Souther Oscillation index for Dec – Apr the year survey taken
- BEUTI_MeanSeas: Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- BEUTI_SDSeas: Std Deviation of the Biologically Effective Upwelling Transport Index for March to survey month the year survey taken
- CUTI_MeanSeas: Coastal Upwelling Transport Index for March to survey month the year survey taken
- CUTI_SDSeas: Std Deviation of the Coastal Upwelling Transport Index for March to survey month the year survey taken
- MeanDailyWaterTemp_AllYear: Average daily water temperature for the year the survey was taken in Celsius
- MeanMaxDailyPeakAirTemp_PeakSeason: Average daily maximum air temperature from April to the survey month for the year the survey was taken in Celsius
- SDDailyWaterTemp_AllYear: Std Deviation of the Average daily water temperature for the year the survey was taken in Celsius
- SDMaxDailyPeakAirTemp_PeakSeason: Std Deviation of the Average daily maximum air temperature from April to the survey month for the year the survey was taken in Celsius
- YearSinceSSWD: Year since the outbreak of sea star wasting disease
- AreaUnique: Unique identifier for the area concatenated from SiteCodeZoneExposureCodeAreaYear
- Siteyear: concatenated Site and Year
- capeyear: concatenated Cape and Year
- SiteyearZone: concatenated Site and Year and Zone
Sharing/Access information
Licenses/restrictions placed on the data: None
Links to publications that cite or use the data: TBD
Recommended citation for this dataset: Resistance of rocky intertidal communities to oceanic climate fluctuations. Sarah A. Gravem, Brittany N. Poirson, Jonathan W. Robinson, and Bruce A. Menge
Code/Software
See methods in manuscript for analysis details and software used
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, 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.25m2 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 log10-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 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 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.
Usage notes
Microsoft Excel
R and R Studio
PRIMER-e with PERMANOVA+ Add-on