Ecological theory predicts that kelp forests structured by trophic cascades should experience a recovery and persistence of their foundation species when herbivores become rare. Yet, climate change may be altering the outcomes of top-down forcing in kelp forests, especially those located in regions that have rapidly warmed in recent decades, such as the Gulf of Maine. Here, using data collected annually from 30+ sites spanning >350 km of coastline, we explored the dynamics of Maine’s kelp forests in the ~20 years after a fishery-induced elimination of sea urchin herbivores. Although forests (dominated by Saccharina latissima and Laminaria digitata) had broadly returned to Maine in the late 20th century, we found that forests in northeast Maine have since experienced slow but significant declines in kelp, and forest persistence in the northeast was juxtaposed by a rapid, widespread collapse in the southwest. Forests collapsed in the southwest apparently because ocean warming has – directly and indirectly – made this area inhospitable to kelp. Indeed, when modeling drivers of change using causal techniques from econometrics, we discovered that unusually high summer water temperatures the year prior, unusually high spring water temperatures, and high sea urchin densities each negatively impacted kelp abundance. Furthermore, the relative power and absolute impact of these drivers varied geographically. Our findings reveal that ocean warming is redefining the outcomes of top-down forcing in this system, whereby herbivore removal no longer predictably leads to a sustained dominance of kelp – the ecosystem’s foundation – but instead has led to a waning dominance (northeast) or the rise of a novel phase state defined by “turf” algae (southwest). Such findings indicate that limiting climate change and managing for low herbivore abundances will be essential for preventing further loss of the vast forests that still exist in northeastern Maine. They also more broadly highlight that climate change is “rewriting the rules” of nature, and thus that ecological theory and practice must be revised to account for shifting species and processes.

Methods for: Quantifying kelp forest trajectories through time and space

Associated files:

DMR_benthicsurvey_algae_urchin_randomsites_alldepths_2001_2018.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2016.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2017.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2018.csv

To quantify changes in Maine’s kelp forests through time and space, we conducted in situ surveys (via SCUBA) of rocky reefs every year from 2001 to 2018 (via long-term monitoring by the Maine Department of Marine Resources; hereafter “ME-DMR”).  From May to June of each year, ME-DMR divers visually surveyed up to 140 randomly selected sites, some of which were sheltered from waves and/or in rivers or estuarine environments and were subsequently excluded from our analysis (see below and associated code).  At each site, the diver began their survey at 15 m depth (mean lower low water, hereafter “MLLW”) and swam perpendicular to shore, placing a 1 m² quadrat at predetermined intervals (associated with 15 m, 10 m, and 5 m depth isobaths). The diver recorded aggregate percent cover of subsurface canopy-forming macroalgae and understory macroalgae (respectively) within the quadrat (n = 10 quadrats per depth per site, except three sites where n = 3-10).  Subsurface canopy-forming macroalgae included the kelps Saccharina latissima, Laminaria digitata, Alaria esculenta, and Agarum clathratum (order Laminariales) as well as two Desmarestia species.  However, because Desmarestia accounted for a very small amount of macroalgae surveyed (% cover in all species-specific surveys from 2016-2018: mean: 2.6; SD: 4.2; range: 0.0 to 26.3; see associated code), we hereafter refer to the ME-DMR canopy-forming macroalgae data as “kelp”.  To specifically align ME-DMR survey data with the focal areas of past studies (i.e., wave-exposed reefs located on offshore islands and at the ends of peninsulas), we scored the wave exposure of each site and its distance from shore (using a scale of 1 to 5 for each metric), then removed all sheltered sites (exposure score of 2 or less) and inland sites (coastal score of 2 or less).  To classify the wave exposure of each site, we estimated the total arc (in degrees) where a site was exposed to the open ocean (i.e.,  a ‘fetch’ of greater than 10 km).  Sites with a score of 1 had no open ocean fetch whereas those with a score of 5 had close to 180º of exposure to open ocean. To classify distance from shore, we set the exposed mainland coast at a level of 3 and then assigned offshore islands with higher scores and sites within rivers or estuaries with lower scores. Site classifications including location, exposure score, and distance from shore score can be found in the datasets. While rare, any quadrats conducted in soft-bottom or unstable (i.e., pebble) habitat were excluded. 

For analysis, we divided the coast into six sub-regions based on hydrology, geographic breaks (e.g., embayment) and local convention.  From southwest to northeast these are York, Casco Bay, Midcoast, Penobscot Bay, Mount Desert Island (hereafter “MDI”) and Downeast.  Casco Bay sub-region extends from Cape Elizabeth to Cape Small and is a well-defined embayment.  Likewise, Penobscot Bay extends from Port Clyde to Isle au Haut and is another large, well-defined embayment.  Downeast Maine extends from the Schoodic Peninsula to the U.S.-Canadian border.  The Midcoast and MDI sub-regions are interspersed between these bays. The York sub-region extends from Cape Elizabeth to the Maine-New Hampshire border. These sub-regions also broadly follow existing state fishery management zones (e.g., Downeast and MDI correspond to zones A and B for the American lobster fishery, respectively).

Excluding sites for analysis – based on the aforementioned criteria – produced a filtered ME-DMR dataset containing n = 3 to 31 sites per sub-region per year (except 2015-2018, when York and Casco Bay were not surveyed).  We augmented this filtered ME-DMR dataset with other, species-specific algal community data collected by the authors in 2016, 2017, and 2018 from various sub-regions (including York and Casco Bay) that met the same criteria (for data collection methods, see next section), to bolster our analyses and inference (see associated code).  The resultant combined dataset contained between 31 and 67 wave-exposed, outer coastal study sites per year for analysis.  Within the dataset, we only analyzed data collected at  5 m depth, because (a) only this depth was consistently surveyed across years and programs and (b) it predictably harbors stable substrate (ledge, boulder) – a prerequisite to forest development.

Methods for: Quantifying changes in algal community composition and abundance

Associated files:

Rasher_Steneck_benthicsurvey_algae_allsites_alldepths_2004.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2016.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2017.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2018.csv

To quantify changes in algal community composition through time and space, we conducted high resolution in situ surveys (via SCUBA) of rocky reefs along the coast of Maine during the summers of 2004, 2016, 2017, and 2018.  Survey sites in 2004 (n = 6 to 31 sites per sub-region; 100 total) were selected using a random-number generator to determine the site’s longitude within a predetermined segment of coastline; the most wave-exposed point closest to this longitude became the dive site.  In 2016, 2017, and 2018, we randomly resampled a subset of the 2004 sites.  In 2016, we resurveyed five sites in three sub-regions (Midcoast, MDI, and Downeast) and in 2017, we resurveyed 22 sites in three sub-regions (Penobscot Bay, MDI, and Downeast).  In 2018, we resurveyed all sub-regions (n = 29 sites, 5 to 6 sites per sub-region) except for York; data for York in 2018 were acquired from author JEKB (https://catalogue-temperatereefbase.imas.utas.edu.au/geonetwork/srv/eng/catalog.search#/metadata/a869dfaa-59a7-4d2f-aaa6-47c2dae5741a).  Sites were sampled at 5 to 7 m depth (MLLW; hereafter “5 m depth”) as well as at 10 m depth MLLW each year.  However, when analyzing changes in community composition, we focused only on data from 5 m depth (for reasons described above).  The one deviation from this approach was for York, where surveys were conducted at 7 to 10 m depth.  As with ME-DMR data, prior to analysis we scored the wave exposure of each site and its distance from shore, then removed any surveys of inland sites, sheltered sites, and the few quadrats conducted in soft-bottom or unstable habitat, from each dataset (see associated code).

At each site and depth, a diver quantified the species identity and percent cover of all kelps and other canopy-forming macroalgae (i.e., Desmarestia spp.) found within a 1 m² quadrat.  This process was repeated at predetermined intervals along a 10 m transect (n ≥ 8 quadrats per depth per site).  In addition, within a 0.25 m² subsection of every quadrat, the diver assessed the identity and percent cover of all macroalgae that defined the understory.  Lastly, within each 1 m² quadrat, the diver categorized the substrate type [ledge, boulder, cobble, etc., using the scale in Wentworth 1922 (https://doi.org/10.1086/622910).  The only deviation from this protocol was in York in 2018, where we instead used a Universal Point Contact method (n = 80 points along a 40 m transect). 

Kelps were identified to species.  Understory algae were identified to genus or species, except for filamentous red algae, which were grouped.  For the kelp timeseries analysis and associated driver model (see associated code), percent cover estimates of each kelp species were aggregated as total % cover of kelp and combined with the ME-DMR data.  In models assessing changes in algal community composition (see associated code), percent cover estimates were left at the species level or grouped at the genus or higher level for those that could not be resolved in the field (see associated code).  Changes in algal community composition were evaluated only between 2004 and 2018 – the only years when all six sub-regions were surveyed at this higher taxonomic resolution. 

 Methods for: Quantifying sea urchin dynamics through time and space

Associated files:

DMR_benthicsurvey_algae_urchin_randomsites_alldepths_2001_2018.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2016.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2017.csv ;

Rasher_Steneck_benthicsurvey_algae_urchin_allsites_alldepths_2018.csv ;

Urchin_Diameters_2001-2018.xlsx 

We enumerated the density of adult green sea urchins (S. droebachiensis with a test diameter greater than 20 mm) found within each 1 m² quadrat sampled during the annual ME-DMR survey, as well as during the species-specific algal surveys conducted from 2016 to 2018 (but not 2004).  As with algal data, we filtered the combined sea urchin density dataset to include data only from 5 m depth and those sites that met our aforementioned criteria (see associated code).

When sea urchins were present, we also measured their size (test diameter, to nearest mm) in a subset of the ME-DMR quadrats (n = 2 per depth per site) or – in the case of the high-resolution algal surveys – all of the quadrats surveyed from 2016-2018 (n ≥ 8 per depth per site).  As with the density data, we filtered the combined sea urchin size dataset to include data only from 5 m depth and those sites that met our aforementioned criteria (see associated code).

Methods for: Quantifying seawater temperature changes through time and space

Associated files:

NOAA_temperature_allsites_2001_2018.csv

To assess the possible effects of changing seawater temperatures on kelp, we obtained publicly available seawater temperature data for 2001-2018 from NERACOOS and NOAA oceanographic buoys (National Buoy Data Center 1971) located in all six sub-regions.  Measurements were made at depths shallower than our kelp forest surveys (1 m vs. 5 m depth) but were nevertheless utilized here, because (a) in situ seawater temperature data at the reef level are lacking in this region, and (b) during summer, Maine’s coastal thermocline is typically below 10 m depth (Brown and Irish 1993).  Further, empirical buoy measurements are generally more accurate than satellite derived estimates of SST (Vinogradova et al. 2009).  We thus calculated daily mean ‘Near-Surface’ Seawater Temperatures (NSSTs) from hourly measurements taken at 1 m depth.  In Casco Bay, NSST data from NDBC 44007 were primarily used, whereas CASM1 provided data during several weeks when NDBC 44007 was out of service.  The Downeast buoy (NDBC 44027) had multi-day or month-long gaps during 8 of 18 study years.  We interpolated these missing data using nearby MDI data by way of linear regression (R² = 0.943).  Otherwise, data originated from a single buoy in each sub-region (see associated code).  While missing daily data occurred up to 8-12% for some sub-regions, we averaged data together in spring and summer, ultimately only losing 1 to 2 data points per sub-region.

ME-DMR kelp forest surveys occurred primarily in May of each year from 2001 to 2018.  We hypothesized two components of annual temperature change influenced the percent cover of kelp observed during each survey: (1) summer (June 1 to August 31) seawater temperatures experienced the year prior, which may have affected adult kelp biomass (via direct mortality) or impacted reproduction in the prior year, and (2) spring (March 1 to May 31) temperatures in the 2 to 3 months leading up to and during the survey, which may have impacted adult and juvenile kelp growth during this time – the peak of primary production – either via thermal stress or due to reduced nutrient availability.  We thus calculated mean and maximum spring NSSTs, as well as 1-year time lagged mean and maximum summer NSSTs, for each sub-region from 2002 to 2018 (see associated code).  For each, we also calculated sub-regional means across the timeseries, as well as annual deviations (anomalies) from the sub-regional mean, for our model (see associated code).  We considered other temperature calculations, including ‘heat degree days’ but these alternative metrics were too collinear with means to be used in our models.

We repeated this process using ¼º resolution temperature data from the National Oceanographic and Atmospheric Administration’s Optimum Interpolated Sea Surface Temperature (OISST) product (Reynolds et al. 2007, Banzon et al. 2016, Huang et al. 2021) in order to evaluate if more spatially extensive modeled seawater temperature data would provide clearer answers within our models (see associated code).  All results were qualitatively the same as using buoy temperature data and thus we chose to use the buoy data for simplicity.