Changes in Rocky Intertidal Community Structure During a Marine Heatwave in the Northern Gulf of Alaska

Introduction

The ocean’s climate varies naturally over a range of temporal and spatial scales, from seasonal cycles to interannual or interdecadal patterns, such as the El Niño-Southern Oscillation (ENSO) or the Pacific Decadal Oscillation (PDO) to long-term climate and ecosystem transformations termed “regime shifts” (Anderson and Piatt, 1999; Beaugrand, 2004; Litzow, 2006; Litzow and Mueter, 2014). These naturally occurring events overlay the anthropogenic trend of increasing temperature (and associated physical and chemical changes) resulting from climatic forcing mediated primarily by greenhouse gas emissions (Gleckler et al., 2012). Marine biological systems respond to these changes in a multitude of ways with the outcomes ultimately culminating in either an increase or decrease in population abundance or distribution (Fields et al., 1993; Kordas et al., 2011). Expansions of population distribution or increases in population numbers may be some of the positive outcomes and can include many fish (Murawski, 1993) and coral species (Baird et al., 2012). Negative outcomes can include reductions in phytoplankton production (Barber and Chavez, 1983) and associated decreases in forage fish, seabirds, and piscivorous fishes (Piatt et al., 1999, 2020). In nearshore systems, a reduction in canopy forming kelps as well as the grazers of these kelps is a common result of warming (Tegner and Dayton, 1987; Edwards, 2004; Filbee-Dexter et al., 2020).

Determining how various marine communities respond to a warmer ocean is further complicated when extreme events occur. Marine heatwaves are one type of extreme event, defined as prolonged anomalously warm water events that are distinct in their duration, intensity, and spatial extent (Hobday et al., 2016). Marine heatwaves are occurring on a global scale (Di Lorenzo and Mantua, 2016; Couch et al., 2017; Oliver et al., 2017; Benthuysen et al., 2020) and are impacting many different types of habitats (Le Nohaïc et al., 2017; Arias-Ortiz et al., 2018) and organisms (Short et al., 2015; Jones et al., 2018). The severity, coupled with the temporal and spatial extent of marine heatwaves, likely dictates a species’ response. Over the past century, there has been an increase in marine heatwave frequency and intensity (Frölicher et al., 2018). From 1925 to 2016, global average marine heatwave frequency and duration increased by 34 and 17%, respectively, resulting in a 54% increase in annual marine heatwave days globally (Oliver et al., 2018; Holbrook et al., 2020).

In the Northeast Pacific, an extreme marine heatwave occurred from 2014 through 2016. This Pacific Marine Heatwave (PMH) or Sea Surface Temperature Anomaly (SSTA, Arafeh-Dalmau et al., 2019), initially identified as the North Pacific “Blob,” was first observed in offshore waters in January 2014 and had encroached on coastal regions in the spring of that year (Bond et al., 2015) with ocean water temperatures peaking in 2016, exacerbated by an El Niño (Walsh et al., 2018). During 2016, sea surface temperatures (SST) in the Gulf of Alaska (GOA) were among the highest on record, and coastal areas of Alaska had their warmest winter-spring on record (Walsh et al., 2018). The Blob turned into a prolonged SSTA, also dubbed “Blob 2.0” (Amaya et al., 2020), with higher than normal SST resurging in 2019. Reports of impacts in the pelagic zone from the 2014–2017 PMH indicate similar biological responses to previous warming events at lower trophic levels, including a shift in phytoplankton species composition (Batten et al., 2018; Peña et al., 2019), a reduction in phytoplankton production (Gomez-Ocampo et al., 2018), and increases in harmful algal species (Vandersea et al., 2018). These effects also were evident in higher trophic levels with observations of poleward shifts in the distribution for migratory fish species, reductions in forage fish abundance and nutritional composition (von Biela et al., 2019; Barbeaux et al., 2020), along with redistribution and die-offs of sea birds [especially common murres (Uria aalge)] and California sea lions (Zalophus californianus) (Cavole et al., 2016; Walsh et al., 2018; Piatt et al., 2020) and a reduction in baleen whale reproductive output, especially humpback whales (Megaptera novaeangliae) (Cartwright et al., 2019).

Rocky intertidal habitats naturally experience a wide range of physical conditions and their communities are notably resistant to these varying conditions. At low tide, intertidal habitats are exposed to the air, which can desiccate organisms, while they also may be experiencing extreme heat or freezing temperatures (Carroll and Highsmith, 1996). While immersed, intertidal organisms may experience fluctuations in salinity, temperature, and wave forcing. The specific physical attributes that structure a particular intertidal community will vary depending on its location and its associated local environmental conditions (Cruz-Motta et al., 2010; Konar et al., 2016). On smaller spatial scales, zonation patterns along the tidal elevation gradient are common (Underwood, 1985; Hawkins et al., 1992; Bertness et al., 2006) and may be more obvious than regional differences (Konar et al., 2009). While understanding drivers of intertidal community dynamics is a continuing effort (Hacker et al., 2019; Lara et al., 2019), extreme events, such as marine heatwaves, are unique opportunities to study how intertidal communities respond to large-scale perturbations.

Assessing the effects of extreme warming events like marine heatwaves is especially compelling in rocky intertidal systems. The mostly sessile nature of the dominant organisms in rocky intertidal habitats does not allow them to move to thermal refugia as may be possible for open water organisms. Furthermore, although rocky intertidal organisms are likely pre-adapted to extreme environmental fluctuation (Wethey et al., 2011; Vinagre et al., 2016), documenting responses to warming is not new (e.g., Sanford, 1999; Harley, 2008) and some extreme responses in rocky intertidal taxa have been specifically attributed to recent marine heatwave events. For example, a marine heatwave impacting Australia and New Zealand coasts resulted in the loss of the local low intertidal bull kelp, Durvillaea spp. and replacement by the invasive kelp Undaria pinnatifida (Thomsen et al., 2019). During the recent Mediterranean heatwave, the presence of the invasive mussel, Xenostrobus secures, increased the survival of the native mussel, Mytilus galloprovincialis, indicating complex physical and ecological interactions (Olabarria et al., 2016). Strong increases in the recruitment of limpets and barnacles as well as unprecedented numbers of species moving northwards have been observed along the northern California coast during the PMH (Sanford et al., 2019). The recent sea star wasting syndrome outbreak from Mexico to Alaska coincided with the PMH, although causal relationships are uncertain (Harvell et al., 2019; Konar et al., 2019). Here, we examined information from an existing program monitoring intertidal habitats and their associated communities in the northern GOA to assess rocky intertidal community structure prior to, and response during and after the recent PMH. We tested the progression of percent cover of intertidal sessile taxa over time as well as changes in percent cover observed before and after the onset of the PMH. We also determined which species groups were particularly responsive to such warming events. In addition, we assess whether responses were similar across four different regions of the northern GOA.

Materials and Methods

Study Design

We collected species composition and relative abundance of sessile invertebrates and macroalgae (i.e., sessile community data) in rocky intertidal habitats from four regions in the northern GOA, spanning approximately 1,200 km of coastline, annually in the summer from 2012 to 2019 (Figure 1). Starting in 2012, two independent monitoring programs in the northern GOA (Rigby et al., 2007; Dean et al., 2014) were combined and expanded into the Gulf Watch Alaska (GWA) monitoring program (¹ Coletti et al., 2018). The four monitoring regions include Western Prince William Sound, Kenai Fjords National Park, Kachemak Bay, and Katmai National Park and Preserve (Figure 1). Within each region, there were five rocky intertidal sites except for Kachemak Bay, which had six sites. All sites were initially chosen based on the presence of hard substrate, exposure (relatively protected from high wave exposure), slope (not vertical), and extent (at least 100 m continuous rocky habitat), which range in geographic location from small embayments and fjords along the outer coast to more estuarine influenced coves in protected inland waters. Further descriptions of these regions and sites can be found in Konar et al. (2016).

During each sampling event, we estimated percent cover of sessile invertebrates and macroalgae, identified in the field to the lowest feasible taxonomic level (as low as species but sometimes phylum), in quadrats along permanently fixed 50 m transects at the mid and low intertidal strata (approximately at + 1.5 and + 0.5 m, respectively, relative to mean lower-low water). Percent cover was quantified using one of two methods to maintain consistency through time with the prior long-term monitoring programs in each region. Protocols at all sites except for those in Kachemak Bay were originally established as part of the National Park Service Southwest Alaska Inventory and Monitoring Program, then adopted in 2012 by GWA to support continuation of monitoring across the three regions. At these sites, percent cover was determined within twelve 0.25 m² quadrats in each stratum. Quadrats were systematically placed along the transects based on a randomly chosen starting point uniquely selected each year. Within quadrats, the occurrence, by layer (from surface to substrate), of macroalgae and sessile invertebrates was determined at 25 systematically placed points. Percent cover was calculated based on the proportion of layers and points occupied by each taxon, including a minimum percent (1%) for all species present within a quadrat, but not encountered in the 25 points (Dean and Bodkin, 2011). In Kachemak Bay, to preserve consistency with other historical data collected as part of the Census of Marine Life (Rigby et al., 2007), percent cover of the overstory kelp layer (i.e., primarily species of Alaria, Hedophyllum, Laminaria, and Saccharina), if present, and all other sessile invertebrates and algae were visually estimated from ten randomly placed 1 m² quadrats. Within each percent cover dataset resulting from either method, all layers were then combined, and the proportion of each taxon was determined for each quadrat to create a standardized percent cover for each replicate quadrat sample. Consistent observers were used when feasible across regions to minimize observer bias in sampling. In addition, a methods comparison was completed by Konar et al. (2016) to ensure data comparability across all study sites despite differences in the details of data collection methods.

At all sites, water and air (when exposed at low tide) temperatures were monitored with HOBO V2 temperature loggers (accuracy of ± 0.2°C; Onset Computer Corporation, Bourne, MA, United States) during the study period. They were placed inside a short length of perforated (so that water would flow around the sensor) 2.54 cm diameter gray PVC pipe that was securely bolted to a boulder or bedrock. We attempted to place the temperature loggers as close as possible to the 0.5 m mean lower low water tidal elevation, (the low tidal stratum in this study), by observing the water level at the site when the tide was predicted to be at 0.5 m, based on the nearest tide station’s predicted tide chart (Tides and Currents Software, NOBELTEC, Beaverton, OR, United States). Each logger was deployed for 1 year at their initial deployment. Loggers were retrieved and replaced annually thereafter, maintaining the original mounting location. Loggers were programmed to record temperature at 20-, 30- or 60-min intervals.

Prior to 2014, loggers were deployed without testing for possible temperature recording errors. While not required based on manufacturer’s specifications (ONSET Hobo Pro V2), we began testing loggers pre- and post- deployment for temperature reading errors. Any post-deployment errors that were found were applied to the entire temperature record from that deployment. After testing began, any temperature logger that did not meet the ± 0.2°C recording accuracy threshold at any pre-deployment temperature soaks were subsequently removed from the instrument pool. All temperature data presented here, hence, met the quality-control standards.

Since the HOBO temperature loggers were within the intertidal zone, they alternated between being submerged and being exposed to air under the influence of sea level changes (typically twice per day with the dominant semi-diurnal tide, although storm surges, and river freshets also impact sea level). Data are assumed to represent submerged water temperatures when the predicted tide level from the nearest tide station was ≥1.5 m. We excluded temperatures when the predicted tidal elevation was <1.5 m, as these data were likely air temperatures. From the resulting data, we plotted monthly water temperature averages and anomalies (relative to mean across all years of study) for each of the four regions. The length of each temperature time-series varied by site and most sites had some data gaps due to sensor malfunction, loss, or break in deployment.

Statistical Analyses

Percent cover data were analyzed using PRIMER v7 and PERMANOVA + (PRIMER-e ltd, Quest Research Limited; Anderson et al., 2008; Clarke et al., 2014; Clarke and Gorley, 2015). To examine changes in intertidal community structure among regions before and during/after the PMH, community data from mid and low tidal strata were analyzed independently from 2012 to 2019, based on prior work that showed significant community differences among tidal strata (Konar et al., 2009). As we were specifically interested in investigating changes with the onset of the PMH, categorical designations for time periods before (2012–2014) and during/after (2015–2019) the PMH were defined by Northeast Pacific and GOA-wide analyses (Di Lorenzo and Mantua, 2016; Hobday et al., 2018; Cornwall, 2019) and intertidal water temperatures recorded at our sites (Danielson et al., 2020). We designated 2012–2014 as pre-PMH since our summer sampling coincided with the start of the PMH in 2014 and we did not expect there to be a detectable response from the intertidal community in that same year. The standardized percent cover data from all regions were used so that the total percent cover from each sample fell between 0 and 1. Data were fourth-root transformed (√(√(x + 0.01)) and a dummy variable of 0.01 added to construct a zero-adjusted Bray-Curtis similarity matrix for both tidal strata. The dummy variable removes undefined dissimilarities in the matrix but does not influence the overall resemblance structure (Clarke et al., 2006). A similarity percentage analysis (SIMPER) was performed to determine which taxa contributed to at least 70% of observed variation in each region-PMH grouping for both elevations. Similarity in community structure before and during/after the PMH among regions was calculated with the SIMPER routine. A 3-factor permutational multivariate ANOVA (PERMANOVA; McArdle and Anderson, 2001) design was constructed for each stratum with region and year as fixed factors. For sites, all quadrats of a year and stratum were averaged and the site term included in the PERMANOVA as a random factor, nested within region to account for this random interaction, non-zero variance component when testing the terms region and year. The categorical factor PMH was added as a contrast to the factor “year.” The PERMANOVA was run for 9,999 permutations to determine differences in community structure. Pairwise tests using the categorical factor PMH by region were performed to test for significant differences in community structure before and during/after the PMH for each region and among regions. Visualizations of community structure were made for tidal strata separately by calculating average values for each region-year group and ordinating them in a non-metric multi-dimensional scaling (nMDS) plot. To complement comparisons in community composition with univariate diversity measures, we calculated the Shannon Wiener Diversity index and Pielou’s Evenness index for each region before and during/after the PMH and compared data for each region using t-tests. Means plots were used to visually determine PMH-related trends in taxa deemed important by SIMPER analysis.

Results

The PMH was evident in intertidal water temperature records by May 2014 with temperatures returning to pre-PMH levels in early 2017, but then warming again by late 2018 through summer 2019 (Figure 2). The mean monthly changes in temperature were observed synchronously across all regions. Rocky intertidal communities differed by region for both the mid and low intertidal stratum (PERMANOVA, p = 0.0001; largest pseudo-F values for the factor region in both strata, Table 1). The factor year also was significant in both strata, but the larger pseudo-F value associated with the contrast PMH indicated a larger effect of the PMH than interannual variation (Table 1). Similarly, the interaction effect of region x PMH was larger than the region x year interaction for both strata.

Increases in similarity were observed after the onset of the PMH, but rocky intertidal community structure still varied significantly before and during/after the PMH for each of the tidal strata for all regions (Figure 3). Community structure before and during/after the PMH remained significantly different in all regions and strata, except for Western Prince William Sound mid stratum (Table 2). Similarly, all between-region pairings in the before PMH and the during/after-PMH periods were significantly different (Table 3). Univariate diversity indices (Pielou’s Evenness and Shannon Wiener Diversity) changed in some regions and strata before and during/after the PMH (Table 4). Evenness increased in the Katmai mid and Western Prince William Sound low stratum but decreased in the Kenai Fjord mid stratum. Shannon diversity decreased in the Kenai Fjord mid but increased in the Western Prince William Sound low stratum (Table 4).

Before the PMH, Kachemak Bay, and Western Prince William Sound were the most dissimilar among the regions for both the mid and the low stratum (Figure 4). Regional trajectories began to converge more at both tidal strata and through time in the during/after-PMH period. At the mid stratum, Kachemak Bay showed the greatest variability with a relatively consistent trajectory through time, while, at the low stratum, variability through time was similar among regions (Figure 4).

Not all taxa responded consistently coincident with the PMH but there were some clear positive and negative responses (Table 5). A SIMPER analysis revealed that 17 of 81 taxa (or taxonomic groupings), plus bare substrate (open space), could explain most of the observed differences in community structure between before and during/after PMH periods, at both tidal strata across all regions (Table 5). Changes in the percent cover within each tidal stratum for each of the important taxa identified by the SIMPER analysis demonstrated that strata often showed similar trends. For example, mussels (Mytilus trossulus), barnacles (i.e., Balanus glandula, Semibalanus spp., Chthamalus dali), bare substrate, and some perennial species of algae such as the red Odonthalia/Neorhodomela spp. and Savoiea/Polysiphonia spp., increased between the two time periods but most other macroalgae, including the brown alga, Fucus distichus, and red algae Devaleraea/Palmaria spp. and Halosaccion glandiforme decreased between before and during/after the PMH (Table 5).

Although the percent cover of some taxa differed between the before and during/after PMH time periods, how those changes occurred over time were variable (Figure 5). Some changes were seen as gradual increases (Savoiea/Polysiphonia) or decreases (Devalerea/Palmaria, encrusting coralline algae, and Halosaccion). Others had a time lag associated with the change (bare substrate, barnacles, and Ulva/Monostroma). While some taxa decreased initially after the PMH, they appeared to be recovering or starting to recover in the later study years (Cladophora/Chaetomorpha, Fucus, and Odonthalia/Neorhodomela). Lastly, some taxa were associated with high interannual variability (Boreophyllum/Pyropia/Wildemania and non-coralline algal crust). In general, percent cover trends were similar in both strata, except for when a taxon did not have high cover in one of the strata (Alaria, Cryptosiphonia, and Gloiopeltis). The predominant sequence of changes that occurred across all regions included a loss of macroalgae with the onset of the PMH, which opened up bare space in 2016, leading to a peak in barnacles in 2017 and then in mussels in 2018 and 2019 (Figure 5).

Discussion

Intertidal organisms are adapted to and tolerant of extreme temperatures due to the cyclical exposure to air and water; however, sometimes when acting in concert with other processes, water temperature can drive the structure and function of these systems (McQuaid and Branch, 1984; Sanford, 2002; Wolfe et al., 2020). Extreme warm water conditions, such as experienced during marine heatwaves, can exceed the thermal tolerance of shallow-water marine taxa and lead to the loss of habitat-forming taxa and the homogenization of marine communities (Colossi Brustolin et al., 2019). For example, foundational seagrasses (Thomson et al., 2015) and kelp (Rogers-Bennett and Catton, 2019; Arafeh-Dalmau et al., 2020; Filbee-Dexter et al., 2020) collapsed during marine heatwaves, leading to persistent losses of ecosystem function and services. The present study adds to our understanding of such marine ecosystem responses, where intertidal foundation species such as rockweed (Fucus distichus) and other fleshy macroalgae decreased during the PMH. These losses coincided with a shift from an autotrophic to a more heterotrophic-dominated community, as filter feeders such as barnacles and mussels increased after the onset of the PMH. This shift was characterized by a greater homogenization (greater community similarity) of the intertidal community, as also observed for other marine systems elsewhere (Wernberg et al., 2013; Colossi Brustolin et al., 2019). In our study, however, community homogenization remained evident over a 5-year period during/after the PMH.

The increase in seawater temperature from the PMH across the northern GOA was detected as positive warm-water anomalies at our rocky intertidal sites starting in May 2014 that persisted through the end of 2017 (Figure 2; Danielson et al., 2020). Changes from generally colder to warmer water anomalies were coincident with major changes in community structure at our intertidal sites. Most effects of the PMH on intertidal taxa lagged the onset of the PMH event by at least a year; however, this observation may have been an artifact of the timing of our sampling in early to mid-summer, causing any effects of the PMH that may have occurred later in the summer or fall of 2014 to remain undetected until the following year. Despite the timing artifacts, there was an obvious change in rocky intertidal community structure. While we present only correlative evidence for a relationship between changes in community structure and the PMH, the scale of detectable responses suggests that changes in structure were attributable (either directly or indirectly) to the PMH, a major disturbance that caused synchronous responses across a large geographic scale.

Severe habitat degradation due to the loss of complex, three-dimensional foundation species in the wake of marine heatwaves has been reported for multiple subtidal systems across the world (Ainsworth et al., 2020). Specifically, kelps (Arafeh-Dalmau et al., 2020; Filbee-Dexter et al., 2020), seagrass (Arias-Ortiz et al., 2018; Strydom et al., 2020), and corals (Le Nohaïc et al., 2017; Fordyce et al., 2019) are biogenic habitat formers that have been recorded to suffer from marine heatwaves. In the rocky intertidal, the rockweeds are important habitat formers (Wikström and Kautsky, 2007), and Fucus distichus was one of the taxa that declined concurrent with the onset of the PMH. Fucus cover varies naturally based on life-history cycles over longer time frames in the GOA (Driskell et al., 2001; Klinger and Fukuyama, 2011), similar to other high-latitude systems (Vadas et al., 2004). However, against this background of natural variation, strong disturbances have elicited synchronous responses in Fucus spp. before, both to experimental manipulation (Klinger and Fukuyama, 2011), climate warming (Álvarez-Canali et al., 2019), and to anthropogenic stressors such as oil contamination (Highsmith et al., 1996; van Tamelen et al., 1997; Driskell et al., 2001). The synchrony of response in these previous studies scaled to the severity of the disturbance, with the most synchronous responses in direction (decrease or increase in cover) and strength coinciding with more severe, catastrophic events. We observed a similar response in Fucus cover, which declined across our study regions simultaneously, indicating that the PMH was likely severe enough to elicit such synchrony. While intertidal species naturally experience and are often adapted to temperature variation (e.g., Somero, 2002), many intertidal species already live at the upper limits of their temperature tolerance, which might make them particularly susceptible to severe warming events such as during the PMH (Lindstrom et al., 1999; Helmuth et al., 2006; Tomanek, 2010). Studies of Fucus have reported significant responses to changes in temperature ranging from physiological responses such as an increase in heat shock proteins (Smolina et al., 2016) and decrease in photosynthesis and growth (Graiff et al., 2015) to population-level responses through decreased post-settlement survivorship (Alestra and Schiel, 2015). However, Fucus populations are often genotypically distinct and have high phenotypic plasticity in their stress responses, representative of low dispersal and local-scale adaptations (Wahl et al., 2011; Smolina et al., 2016). The consistent response (decline; Figure 5) seen in this species across a large spatial scale emphasizes the overwhelming strength of an environmental stressor at play. Aside from recovery of this foundation species, another long-term concern of the PHW influence on Fucus is a possible genetic homogenization (Nicastro et al., 2013) that erodes the species’ typical resilience to a broad range of temperatures (Jueterbock et al., 2018; Coleman et al., 2020).

The decline of other foliose macroalgae (e.g., Devaleraea/Palmaria spp., Halosaccion glandiforme, Mastocarpus/Mazzaella; Figure 5) coupled with the declines in Fucus co-occurred with an increase in relative and absolute abundance of primary space-occupying invertebrate filter feeder species, specifically first to barnacles and then to Mytilus. Essentially, the community shifted from a mostly autotroph-dominated system to a mostly heterotroph/invertebrate-dominated system. Invertebrates, especially barnacles and mussels, also are important rocky intertidal foundation species and play an essential role in structuring the rocky intertidal through both biological and physical processes (Connell, 1961; Lubchenco and Menge, 1978; Petraitis and Dudgeon, 2020). These invertebrates increase habitat complexity in the intertidal and subsequently increase biodiversity (Archambault and Bourget, 1996; Mazzuco et al., 2020). However, while the function of foundation species per se may not be lost by a perturbation such as the PMH, community shifts from one to another type of foundation species typically are associated with shifts in overall biodiversity and ecosystem functioning that will need to be evaluated over time (Sorte et al., 2017). Therefore, while the PMH reduced fleshy algal foundation species, some invertebrate foundation species increased in abundance, either through reduced competition for space, reduced predation from the loss of top predators such as sea stars (Konar et al., 2019), or physiological benefits from warmer waters.

Barnacles are adapted to the temperature fluctuations in rocky intertidal systems and severe or irreversible physiological responses only occur at extremely high temperatures (Berger and Emlet, 2007). This tolerance of long-lived, adult barnacles to warm temperatures may explain their persistence during the PMH. Barnacles can be the first settlers after a perturbation (Highsmith et al., 2001; Traiger and Konar, 2018), and post-settlement barnacles can be tolerant of elevated temperatures (Findlay et al., 2010) but decreases in survival at high temperatures during the larval phase have also been observed (Kendall et al., 1985). Although barnacle cover increased after the PMH, if ocean warming trends continue, then barnacle recruitment may become limited over time, and overall population declines may become evident after a lag time or extended warm periods. This could explain the drop in barnacle abundance we observed after 2017, several years after the onset of the PMH.

The ecological role of mussels (Mytilus) in rocky intertidal systems was first assessed as a dominant space competitor that would reduce overall community diversity when occurring as monocultures (Paine, 1966). Later, however, they were recognized as important foundation species (Menge, 1976) that provide habitat and can increase diversity. The role and success of mussels in intertidal systems often are a result of complex ecological interactions and physiological responses (Olabarria et al., 2016; Moyen et al., 2019; Seuront et al., 2019). In our study, the observed increase in Mytilus cover may have been a result of a major recruitment event that was observed at the GOA-wide scale following the PMH (Bodkin et al., 2018). This was hypothesized to be driven by favorable conditions during larval life stages and high post-settlement survival. Thermal tolerance of newly recruited mussels is higher in cohorts that experience higher temperatures during the dispersal phase, indicating that selection occurs during dispersal (Sorte et al., 2018). Once settled, Mytilus can be more sensitive to desiccation than to heat (Jenewein and Gosselin, 2013). Temperature-tolerant phenotypes may have been able to disperse across GOA regions during the PMH and occupy the open space left by the decline of Fucus and other macroalgal species. While juvenile Mytilus survival is higher when associated with other species that influence microhabitat and climate such as macroalgae or other invertebrates (de Nesnera, 2016; Barbier et al., 2017), Mytilus abundance can also be dependent on predator abundance as they are important prey for sea otters, shorebirds, sea ducks, and sea stars (Marsh, 1986; Miller and Dowd, 2019). Mytilus abundance is currently high in the GOA, possibly because predator numbers of smaller mussels (e.g., sea stars) are low after recent sea star wasting impacts (Konar et al., 2019). Without these predators to control settling Mytilus, recruits might have persisted at much higher abundances.

Gulf of Alaska intertidal community changes not only represent a shift in the physical and taxonomic structure of rocky intertidal habitats but also in the functional structure of the community (Bremner et al., 2006; Scrosati et al., 2011). The switch from an autotroph- to heterotroph-dominated system discussed above will result in associated shifts in ecosystem function (Sorte et al., 2018). For example, a decline in species or functional richness or shift in relative abundances can lead to an imbalance of resource uses (under- or overutilization) in the system (Naeem et al., 1994). Both barnacles and mussels consume a substantial portion of macroalgal production in the form of detritus in addition to phytoplankton in many rocky intertidal systems worldwide (Bustamante and Branch, 1996; Tallis, 2009), including Alaska (Duggins et al., 1989). Typically, these consumers have sufficient trophic plasticity to compensate for the decline in one of the food sources; however, multiple food subsidies are considered essential elements of long-term food web stability (Huxel et al., 2002). Hence, if our observed loss of significant amounts of macroalgae were to persist over long time periods, destabilization of trophic links and community function could be expected.

The dispersal of heat-tolerant species such as discussed here for barnacles and Mytilus can lead to the homogenization (increased similarity) of communities across regions after a heatwave (Eggers et al., 2012; de Boer et al., 2014). Homogenization is indicative of strong metapopulation connectivity, where dispersal across regions supports increased biomass of select species with favorable traits (e.g., thermal tolerance) (Norberg, 2004; Webb et al., 2010). In GOA rocky intertidal habitats, significant differences in intertidal community structure existed across sites prior to the PMH. These differences in structure were likely shaped primarily by local-scale biotic and abiotic drivers (Konar et al., 2016). While some regional differences continued to persist after the PMH, overall similarity in community composition increased, indicating that the relative strength of local drivers of community structure was lessened by the PMH perturbation. Following the principles of community assembly theory, rocky intertidal systems can be hierarchically shaped by processes from system-wide, regional, to local scales (Martins et al., 2008; Weiher et al., 2011; Menge et al., 2015). Environmental conditions often separate rocky intertidal communities on regional or system-wide scales, if those conditions are sufficiently different (Blanchette et al., 2008; Martins et al., 2008; Merder et al., 2018). Site level drivers such as wave exposure, local disturbance, predation and competition can be important local drivers of rocky intertidal community structure (McQuaid and Branch, 1984; Harley, 2006; Cruz-Motta et al., 2010). The pre-PMH dissimilarity among our study regions has been reported previously, but static environmental conditions such as slope, substrate composition, exposure, distance to glaciers, and freshwater, only played a minor role in driving differences in community composition (Konar et al., 2016). This suggests that dynamic environmental conditions such as temperature, salinity, nutrient availability, and dissolved oxygen are important local-scale drivers of community structure. In addition to dynamic physical drivers, biological drivers such as predation influence local rocky intertidal community composition. For example, there was a marked decline in sea star abundance in nearshore waters that was observed coincident with the onset of the PMH, likely due to sea star wasting syndrome (Harvell et al., 2019; Konar et al., 2019). The decline in sea stars, especially species known to be effective mussel predators (Pisaster ochraceus and Evasterias troschelii), was coincident with a subsequent increase in mussel abundance (Coletti et al., 2018). During/after the PMH, increased similarity of the rocky intertidal community structure among the study regions suggests that strong, large-scale oceanographic events like the PMH may override local drivers to homogenize community structure across regions (also see Holbrook et al., 2019). Whether the changes were the direct result of high temperatures or were due to other associated physical or biological changes remains to be determined.

While increases in similarity occurred across GOA regions with the onset of the PMH, those increases in similarity were not homogeneous across tidal elevation strata. The mid elevation stratum community structure became more similar across regions and time than the low elevation stratum, potentially due to the smaller and somewhat different species pool at the mid stratum than at the low, where kelps and other fleshy macroalgae were more abundant. Species differences between these two tidal elevation strata and variability in how these strata respond to biotic and abiotic interactions has been reported elsewhere (Broitman et al., 2001). The changes among particular mid stratum species may be a result of these species being adapted to the high environmental variability typical of mid elevation intertidal communities. Whereas, the community at the low elevation stratum experiences less variability and is not as resilient to temperature disturbances as communities higher in the tidal zone (Somero, 2002). However, species that occur in both strata and general trends in community structure at each strata displayed similar trajectories in response to the PMH, i.e., a decrease in autotrophs and an increase in heterotrophs.

In summary, we observed the transition of rocky intertidal foundational species from autotrophs to heterotrophs coupled with the homogenization of community structure across four distinct regions, spanning over 1,200 km of coastline in the GOA after the onset of the PMH. Increases in similarity were driven by the loss of macroalgae leading to increases in open space (bare substrate), which allowed for the settlement of pioneer foundational heterotrophic species dominated first by barnacles and then mussels, restarting a well-documented processes of succession in rocky intertidal communities (Sousa, 1979; Farrell, 1991) that continued after the PMH had dissipated (Figure 6). While many taxa did not show any obvious or directional responses to the PMH, a few species highlighted in this paper showed strong responses. Strongest responses were detected in those taxa that dominate overall percent cover in the system, their high abundance emphasizing the effect. Other, less-abundant taxa certainly also responded to environmental perturbations but while their relative change may have been high, their absolute impact on community composition was limited. This study documents an important shift in intertidal communities associated with a marine heatwave.

Recent changes in water temperature, first with the slight cooling in 2017 but now again warming conditions (Amaya et al., 2020), may impede recovery of the rocky intertidal communities to a pre-PMH state if a similar disturbance response presented here is triggered to this new warming. Continued monitoring of intertidal community structure may elucidate further progression and key environmental and biological drivers. Understanding how natural variation intertwines with patterns of community response to a disturbance is critical when assessing the causes and drivers of trends in an ecosystem. Our findings indicate that warming water temperatures can trigger intertidal communities to undergo a predictable progression pattern across large spatial scales, but that local biotic and abiotic drivers will further influence how a community responds to and recovers from disturbances.

Data Availability Statement

Publicly available datasets were analyzed in this study. All data used for this study can be found through the Gulf Watch Alaska data portal (https://portal.aoos.org/gulf-of-alaska#) under the Nearshore Component. Alternatively, source data from western Prince William Sound, Kenai Fjords, and Katmai regions can also be found for rocky intertidal community percent cover 2006–2019 at (https://doi.org/10.5066/F7513WCB), and for temperature, 2006–2014 at (https://doi.org/10.5066/F77S7KXH), and 2014–2016 at (https://doi.org/10.5066/F7WH2N3T).

Author Contributions

BW, BK, KI, HC, DM, RS, and TD contributed to the conception of the study. BW, KI, and BK completed all analyses. BW, BK, KI, HC, DM, RS, DH, and ML contributed sections of the manuscript and provided critical review. All authors contributed to field research in remote wilderness areas.

Funding

This research described in this paper was in part, supported by the Exxon Valdez Oil Spill (EVOS) Trustee Council. However, the findings and conclusions presented by the authors are their own and do not necessarily reflect the views or position of the Trustee Council. This synthesis was supported by the EVOS Trustee Council for the Gulf Watch Alaska Program Nearshore Component through the following grants: 16120114-R & 16120114-X (2012–2016), 17120114-H (2017), 18120114-H (2018), 19120114-H (2019), and 20120114-H (2020).

Conflict of Interest

TD was employed by the company, Coastal Resources Associates.

The remaining authors declare that the research was conducted in the absence o any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank the many students, researchers, and volunteers for assisting with intertidal surveys over the years. We thank Distinguished Professor Dr. Marti Anderson, Massey University NZ, for her invaluable advice on statistical analysis. We also thank the scientific team of the Gulf Watch Alaska long-term monitoring program for their expertise and resources in support of this publication. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U. S. Government. This manuscript has been approved for publication consistent with USGS Fundamental Science Practices (http://pubs.usgs.gov/circ/1367/).

Footnotes

1. ^ https://gulfwatchalaska.org/monitoring/nearshore-ecosystems

References