Our study sites encompass three inhabited islands that are part of the Baja California Pacific Islands biosphere reserve in Mexico: Isla Todos Santos (ITS), Isla San Martin (ISM), and Isla San Jeronimo (ISJ) (Figure 1). Located in the CCS, these islands stretch along a 250 km latitudinal gradient (29.8°–31.8°N) near the southern biogeographic distribution limit of boreal populations of M. pyrifera. They are internationally recognized for their high diversity, the abundance of their flora and fauna, the relative natural integrity of their ecosystems, and their high levels of endemism (Samaniego-Herrera et al., 2007). The marine areas adjacent to these islands, support various fisheries, most of them targeting species dependent on M. pyrifera as a source of food and habitat (Montaño-Moctezuma et al., 2014).

We obtained satellite-derived data of sea surface temperature by using the GHRSST Level 4 AVHRR Optimally Interpolated Global Blended Sea Surface Temperature Analysis (GDS version 2), from the National Centers for Environmental Information (2016), available at the physical oceanography distributed active archive center (PODAAC) of the NASA Jet Propulsion Laboratory (Last access: 2019-01-19). Daily sea surface temperature data are available at a 0.25° spatial resolution for 37 years (from September 1981 to present) (Reynolds et al., 2007; Banzon et al., 2016). To have a historical perspective of the strength of the 2014–2016 warming event, we calculated sea surface temperature anomalies for each 0.25° grid cell by subtracting the 1982–2018 local mean. We used a Fourier low-pass filter to show the contribution of interannual (basin-scale) processes (e.g., El Niño and the Blob) by removing signals with periods shorter than 18 months (cut off frequency), retaining inter-annual variability.

We identified MHWs, defined by Hobday et al. (2016) as discrete prolonged anomalously warm-water events in a location, using the same daily sea surface temperature database for the 1982–2018 period. As with other studies of MHWs, we focused on upper percentiles of SSTAs, rather than on absolute SSTs. This approach accounts for the likelihood that marine populations are adapted to local temperature conditions, with anomalies indicating deviation from normal conditions that might elicit location-specific responses to thermal stress. We used the Matlab toolbox (Zhao and Marin, 2019) to identify periods of time where temperatures were above the 90th percentile threshold relative to the climatology (seasonally varying mean), and with a duration of at least 5 days. We identified the start and end dates and calculated the duration and intensity for each MHW event.

We collected ecological data from the three islands. At each island, two sampling sites were selected (i.e., northern and southern sites). Within each of these sites, we used underwater survey techniques (depth ranges of 10–15 m) to characterize the assemblages of invertebrates, fish, and algae at three transects. Data were collected by experienced, trained divers to ensure constancy in survey methods, and species lists compiled over the years. We estimated abundance of M. pyrifera and fish in situ by census transects and assessed abundance of invertebrates and algae by analyzing video transects. To corroborate the results of the in situ estimate of M. pyrifera abundance and better understand the effects of the MHWs, we estimated kelp biomass using remote-sensed Landsat imagery. We conducted the field surveys in late autumn 2013 and 2016, but due to bad weather in 2016, we had to extend the field work to early winter 2017 (henceforth “2016”).

We used multispectral Landsat imagery to estimate M. pyrifera canopy biomass at 30-m resolution for a 10-year period (from January 1, 2009 to December 31, 2018) following Cavanaugh et al. (2011). We analyzed each atmospherically corrected Landsat image using multiple endmember spectral unmixing (MESMA) (Roberts et al., 1998). This analysis enables estimates of kelp canopy area and biomass using observed relationships between MESMA-derived canopy fractions and diver-measured canopy biomass (Cavanaugh et al., 2011). These relationships were developed using data from the Santa Barbara Channel. There are 2 Landsat “scenes” that cover our study sites, and each scene has about 60 images between 1984-present, allowing us to determine seasonal canopy biomass for each 30 m × 30 m pixel.

For all analyses comparing data collected before and after the MHWs, we used the parallel approach of modeling the response as a function of the crossed fixed factors Island (three levels: ITS, ISM, and ISJ) and Time (two levels: Before and After). In each instance, we assumed the six replicates taken at each time at each island (three replicates taken from two sites at each island) to comprise an independent, representative sample for that island. For univariate counts, we fit generalized linear models (GLMs) with Poisson errors, and simplified these based on standard log-likelihood ratio tests. For these models, we assessed approximate significance of pair-wise differences based on 95% confidence intervals. For multivariate responses, we used permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2001) in the R package vegan (Oksanen et al., 2013) on the basis of Bray-Curtis distance matrices derived from square-root transformed data, and used non-Metric Multidimensional Scaling (nMDS) to visualize multivariate patterns in community structure. Where post hoc tests were required, we employed pairwise PERMANOVAs, with false-discovery-rate adjustment for multiple tests (Benjamini and Hochberg, 1995), using functions in the R package RVAideMemoire (Hervé, 2018). We used Dufrene-Legendre Indicator Species Analysis (Dufrêne and Legendre, 1997) to distinguish which species were contributing to significant clusters identified by these post hoc tests, using functions from the R Package labdsv (Roberts, 2018).

The time series of Landsat estimates of kelp biomass were the only long-term data available to explore whether observations from the period after the MHWs were beyond the bounds of natural variation for the period before the MHWs. We could not use standard linear modeling approaches in this instance, because observations would be expected to contain significant temporal autocorrelation, violating assumptions of independence. We therefore used automatic algorithms from the R package tsoutliers (López-de-Lacalle, 2019) to detect outliers from ARIMA (0,1,1) time series fitted to the log₁₀(x + 1)-transformed data from each island. This approach has the advantage over a random walk model of smoothing the fit using an exponentially weighted moving average of preceding values, rather than placing excessive emphasis on the most recent observation. Outliers detected in this way can be assumed to arise from a process not inherent to that which generated the other observations in the time series (i.e., they are outside the natural variation associated with the time series).

We programmed all analyses in R version 3.5.2 (R Core Team, 2013), setting α = 0.05.

Along the Baja California coast, 2013 was in a slightly cold stage, with negative (cold) anomalies and no evidence of the “Blob” (Figure 2A), until 2014, when the warming event became evident (Figure 2B). By 2015, the additional El Niño effect boosted the warm anomalies and extended their distribution offshore and along the coastal waters of Baja California (Figure 2C). Toward the end of 2016, this effect had weakened (Figure 2D), but warm water anomalies remained throughout 2017 (Figure 2E and Supplementary Figure S1). After removing the annual cycle, the average low-pass anomalies were more than 1°C warmer than the extreme El Niño 1982–1984 and 1997–1998 events (Supplementary Figure S1). The 2014–2016 warm-event anomaly registered on the three islands was the largest in extent and most prolonged in duration observed to date, with SSTAs slightly lower at the northernmost site (ITS) than the southernmost sites (ISM and ISJ), with anomalous temperatures peaking on summer 2015 at 2.15, 2.32, and 2.39°C above average, respectively (Supplementary Figure S1).

Around 50% of the identified MHW days between 1982–2017 registered at the three study sites occurred during 2014–2017 (Figure 2F and Supplementary Figure S3A). The period 2014–2015 accounted for 545 MHW days at ITS (58%), 466 at ISM (51%), and 476 at ISJ (53%) (Figures 3A,C). The year with the highest number of MHW days registered was in 2015 for ITS (261 days), ISM (242 days) and ISJ (256 days) (Figures 3A,C). For ITS, the longest continuous MHW comprised 152 days, lasting from November 2014 until late April 2015, while ISM and ISJ registered events lasting 123 and 124 consecutive MHW days, respectively, between late August and late December 2015 (Supplementary Tables S1–S3). In autumn 2015, the highest maximum temperature during a MHW were identified for the three islands, peaking at 5.6, 5.9, and 5.7°C, respectively (Figures 3A,B and Supplementary Tables S1–S3). These values are higher than the maximum registered during previous strong 1982–1984 (4.1, 3.3, and 3.3°C) and 1997–1998 El Niños (3, 3.4, and 3.7°C) (Supplementary Figures S2, S3B).

We detected significant interactions between Island and Time in our Poisson GLMs models of the numbers of M. pyrifera individuals per transect as well as the numbers of fronds per individual before and after the MHWs (Table 1). For the northernmost island, ITS, the number of individuals increased significantly after the MHWs. In contrast, we found significant decreases at the southern islands, ISM and ISJ, with only a few individuals per transect observed after the MHWs (Figure 4A). For all the islands we found a strong, consistent and significant decline in the numbers of fronds per M. pyrifera individual (Figure 4B).

Remote-sensed kelp biomass estimates were highly variable, but we clearly detected the MHWs signal at both ISM and ISJ (Figure 4C). As was the case for number of M. pyrifera individuals (Figure 4A), kelp biomass at ITS contrasted from the pattern at ISM and ISJ, although in this case, the time series analysis suggested no long-term trend (Figure 4C). Irrespective, it was clear that kelp biomass during and after the MHWs at the southernmost islands, ISM and ISJ, was well beyond the bounds of natural variability observed prior to the MHWs. Qualitatively, satellite-derived kelp biomass estimates allowed us to document the loss of the kelp forests during the MHWs at all islands. Immediately after summer 2015, the satellite images did not detect kelp at the surface of the three islands, and the biomass remained negligible until the end of 2016 (Figure 4C). By winter and spring 2017, kelp biomass at ITS had recovered to previous seasonal values, but at ISM and ISJ there was no recovery, and values stayed very low until the end of the time series (end of 2018). Note that very low values of biomass do not necessarily mean that the kelp forest completely disappeared, since satellite images do not detect sub-surface smaller individuals.

The 2014–2016 warming was the most intense and sustained event, with the highest percentage of MHW days registered, in the 37-year historical series for the CCS. Contrary to past El Niños, which appeared after negative temperature anomalies, the strong 2015–2016 El Niño began with an already warmer ocean (“Blob”) (Jacox et al., 2016). The unprecedented combination of these two events preceded the loss of M. pyrifera after summer of 2015, corresponding with the maximum temperature anomalies, which peaked at 5.9°C warmer than the climatology (this heatwave lasted for around 4 months). These anomalies were ∼2°C warmer than any previous registered (El Niños), most likely with temperatures exceeding physiological thresholds for growth, survival, and reproduction of M. pyrifera (Zimmerman and Robertson, 1985; Hernández Carmona, 1988; Tegner et al., 1997; Hernandez-Carmona et al., 2001). The absence of significant wave disturbance during this event (Reed et al., 2016) suggests a combination of thermal stress and nutrient-poor waters (Whitney, 2015) as the primary drivers behind kelp canopy disappearance in Baja California. Our results are consistent with previous work that reports the sensitivity of M. pyrifera in southern California and Baja California to warm, nutrient-poor waters (Zimmerman and Kremer, 1984; Hernandez-Carmona et al., 2001) during El Niño events (Dayton and Tegner, 1984; Dayton, 1985; Zimmerman and Robertson, 1985; Tegner and Dayton, 1987, 1991; Hernández Carmona, 1988; Tegner et al., 1996; Ladah et al., 1999; Hernandez-Carmona et al., 2001; Edwards, 2004; Edwards and Hernandez-Carmona, 2005; Edwards and Estes, 2006).

Reed et al. (2016) found no ecosystem effects in M. pyrifera kelp forests at Santa Barbara, California, to the same extreme warming event. They argued that the variability of giant kelp biomass was high during the 15 studied years and that the observed decline was within normal variability. However, fewer impacts to kelp forests located at Santa Barbara are expected since they are located ∼1000 km north of the southern distribution limit (i.e., the upper limit of thermal tolerance of M. pyrifera). There, average SSTs are ∼3.3°C cooler than at our study sites, providing a greater population-level thermal safety margin (Pinsky et al., 2019). We found similar results in our northernmost site, ITS, located ∼400 km south of Santa Barbara (nevertheless satellite derived biomass estimates remained very low during autumn 2015 and all 2016). However, in the southernmost sites, ISM and ISJ, located 550 and ∼650 km south of Santa Barbara, respectively, we clearly detected the MHWs signal on M. pyrifera forests. Thus our results support previous studies in the CCS and other regions (Edwards, 2004; Edwards and Estes, 2006; Johnson et al., 2011; Wernberg et al., 2013, 2016; Horta e Costa et al., 2014), which reported kelp forest located near temperate-subtropical transition zones to be more vulnerable to historical warming events than those at higher latitudes. For example, following the MHWs of 2011 in Western Australia, Wernberg et al. (2016) reported a regime shift from a kelp to a seaweed turf dominated ecosystem in sites located near temperate-subtropical transition zones, but little or no change in communities at higher latitudes. Similarly, after 2017/2018 MHWs in New Zealand, Thomsen et al. (2019), found that Durvillaea spp. species located at warmer sites were more affected than those found in cooler ones. Although Reed et al. (2016), challenged the status of M. pyrfiera as a sentinel species for climate change, our work suggests that sentinel status should be evaluated in a site-by-site or regional basis, and reinforces the status of Baja California as a sentinel area for climate change (Hobday and Pecl, 2014).

Macrocystis pyrifera forests in the CCS have proven to be dynamic and resilient ecosystems (Dayton et al., 1992; Filbee-Dexter and Scheibling, 2014), recovering quickly after previous El Niño events. La Niña conditions that usually follow El Niño events bring cold, nutrient-rich waters, representing favorable conditions for M. pyrifera (Dayton et al., 1992; Tegner et al., 1997; Edwards and Estes, 2006). The lack of La Niña conditions after the 2014–2016 warming event and the persistence of anomalously high SST during 2017 (around 50 days registered as MHWs) is probably challenging the capacity of M. pyrifera to recover.

Nevertheless, after the warming event M. pyrifera recovered at ITS, the northernmost site. We found higher densities of individuals than before the MHWs, but with significantly lower number of fronds per individual, indicating the dominance of young M. pyrifera individuals. The presence of some productive individuals (≥16 fronds/individual) implies that M. pyrifera did not lose its reproductive capacity, since sporophylls begin to differentiate and have fertile sori when they have around 16 fronds/indivdual (Dayton et al., 1992). We found a strong effect of the MHWs on M. pyrifera densities and number of fronds per individual at the southernmost sites, ISM and ISJ, suggesting a very low density of fertile individuals. Biomass satellite time-series corroborate that kelp forest at ITS recovered after the MHWs, while at ISM and ISJ biomass values remained below average seasonal values throughout 2018. The three islands have similar mean SSTs, suggesting similar thermal safety margins, however, SSTAs and the maximum intensities registered during the MHWs were slightly higher at the southernmost sites than at ITS, which could explain the different latitudinal responses of the kelp forest to the extreme warming. However, it is also true that the differences in the anomalies were very small (∼0.2°C) between islands and other factors such as local seascape complexity, which might create “microclimates” (Woodson et al., 2018; Pinsky et al., 2019), could also contribute to the observed differences. Previous studies found that deeper sites represented local refugia and individuals living there could find and survive strong warming events such as El Niño, possibly serving as a source of spores and recruits to colonized shallower rocky areas (Ladah et al., 1999). Further research could test whether the observed differences are a result of local seascape variability, latitudinal gradients in SSTAs, or perhaps a combination of both.

The extreme conditions resulted in dramatic changes in the invertebrate, fish, and algal assemblages, with the loss of over half of invertebrate species. The significant reduction or disappearance of M. pyrifera can have profound effects on other members of the community (Graham, 2004) due to a release of resources (space and light) (Reed and Foster, 1984; Dayton, 1985) and a decrease in food availability (Ebeling et al., 1985; Harrold and Reed, 1985). Although we found evident changes, we cannot separate between the effects of high SSTs and loss of canopy on the community, since both are synergistically affecting all food web levels. Importantly, we found that M. pyrifera dynamics at the two sothernmost islands in our study immediately after the MHWs were beyond the limits of natural variation in the previous 5 years, and remained that way until the end of our study in 2018. Although all other metrics we analyzed represented only two periods, a year before and a year after the MHW, this finding that kelp biomass was significantly and persistently impacted by the MHW gives us some confidence in attributing community-wide impacts to the effects of the MHWs.

After the warming event, almost all sessile invertebrates disappeared, and all cnidarians were absent. Sessile invertebrates are highly sensitive to extreme temperature conditions due to their inability to escape from unfavorable conditions (Przeslawski et al., 2008) and their dependence on phytoplankton and suspended organic matter as food sources that might decrease during warming events (e.g., during El Niño events, when upwelling is suppressed). Reports of mass mortality of sessile species and bleaching of corals in temperate and tropical ecosystems (Berkelmans and Willis, 1999; Hoegh-Guldberg, 1999; Cerrano et al., 2000; Marshall and Baird, 2000; Perez et al., 2000; Garrabou et al., 2001, 2009; Caputi et al., 2014; Couch et al., 2017; Hughes T. P. et al., 2017) suggest that SSTs at our sites probably peaked at values in excess of their physiological tolerances (Lejeusne et al., 2010). Moreover, almost all sea star species disappeared or dramatically decreased their abundance. In 2014, an epidemic outbreak of wasting disease from Alaska to Mexico led to mass mortalities of many sea star species (Hewson et al., 2014). Anomalously warm water conditions in the CCS have been suggested to favor the spread of the densovirus linked to this disease and thus increase the mortality rate in sea stars (Eckert et al., 2000; Hewson et al., 2014).

For the entire community, we report the loss or decrease in abundance of many species with cold-water affinities and the appearance or increase in abundance of some species with warm-water affinities. The tropicalization of fish species near temperate transition zones has already been attributed as the primary driver of composition change in the fish assemblage (Wernberg et al., 2013, 2016; Caputi et al., 2014; Horta e Costa et al., 2014; Vergés et al., 2014; Smale et al., 2019). The loss of the ten fish species having the most northerly distribution centroids suggests SSTs peaked above their thermal tolerances and favored species with stronger affinity to warmer waters, such as H. semicinctus (Fishbase). We also reported the presence of the sea-urchin C. coronatus at ITS and ISM, which has the southernmost distribution range of invertebrate species described in this study. Similar patterns favoring C. coronatus have been reported during El Niño 2009–2010, where populations in one site at Oaxaca (Pacific southern Mexico) increased (López-Pérez et al., 2016), and during El Niño 2015 in southern California, with reports of poleward range expansion (Freiwald et al., 2016). Another example is a species from the same genus, Centrostephanus tenuispinus, which was the main winner after an extreme MHW in 2011 on the west coast of Australia (Smale et al., 2017).

The effect of warm water temperatures was also evident in kelp forest subcanopy species L. farlowii and E. arborea. L. farlowii was the characteristic subcanopy species before the MHWs, probably due to the absence (or very low densities) of the competitively superior Pterygophora californica (Tegner et al., 1997). In 2016, abundance of L. farlowii dramatically decreased, despite favorable niche conditions (light and space) after the loss of M. pyrifera forest. The sensitivity of Laminaria species to warm water conditions (Schiel et al., 2004) might explain its decrease after the extreme warming. On the contrary, abundances of E. arborea, known to be favored during El Niño conditions due to its tolerance to warm, nutrient-poor conditions (Dayton et al., 1984; Hernandez-Carmona et al., 2000; Matson and Edwards, 2007), increased at the two southernmost islands.

The influence of the warming event on the understory algae was less evident and more variable than for the rest of the community. For the southernmost islands (ISM and ISJ), we found a decrease in cover of coralline crustose algae, which have been reported to be sensitive to MHWs (Martin and Gattuso, 2009; Smale and Wernberg, 2013; Voerman et al., 2013; Wernberg et al., 2013, 2016; Short et al., 2015; Thomsen and South, 2019). However, the role of heatwaves in regulating abundance of coralline algae is still not clearly understood, since other factors, such as the loss of canopy cover (Irving et al., 2004, 2005; Flukes et al., 2014) or ocean acidification (McCoy and Kamenos, 2015), might also be influential for this group. The presence of coralline crustose algae can be facilitated by canopy-forming kelp species as they exclude other macroalgal competitors (Thomsen and South, 2019). This could explain the decrease of this taxon at the southernmost islands after the MHWs, while at the northernmost site, their cover and the abundance of M. pyrifera remained high. Little attention has been paid to coralline crustose algae, despite representing important habitats, and food sources for diverse epifaunal invertebrates (Graham et al., 2007; Thomsen and South, 2019), as well as providing chemical cues for abalone settlement (Morse et al., 1979; Miner et al., 2006). A generalized increase in the coverage of non-crustose coralline understory algae at all sites suggests that the combined loss of canopy kelp, sessile invertebrates, and coralline crustose algae might have modified light conditions and liberated free space for recruits (Reed and Foster, 1984), which probably favored the rest of the understory assemblage (Reed and Foster, 1984; Dayton et al., 1992; Tegner et al., 1997; Clark et al., 2004).

Finally, we report the appearance of two invasive macroalgal species, U. pinnatifida at ITS and S. hornerii at ISM and ISJ. In Baja California, U. pinnatifida was first detected in 2003 (Aguilar-Rosas et al., 2004) and S. hornerii in 2005 (Aguilar-Rosas et al., 2007). In 2013, both species appeared during monitoring surveys at the protected side of the three islands, but were absent from the exposed sides (our study sites) (Montaño-Moctezuma et al., 2014). U. pinnatifida is an opportunistic cosmopolitan invasive species (Lowe et al., 2000) able to tolerate a wide range of temperatures (Morita et al., 2003), allowing this species to colonize empty spaces left by canopy-forming kelps. Nevertheless, without disturbance, U. pinnatifida rarely establishes under a dense canopy (Valentine and Johnson, 2003; Thomsen and South, 2019; Thomsen et al., 2019). Similar findings have been reported in other parts of the world as a result of both MHWs (Thomsen et al., 2019) and manipulated removal experiments (Thomsen and South, 2019). For example, in summer 2017–2018, strong MHWs in New Zealand caused the decline of the canopy-forming kelp species Durvillaea spp., and a successive subsequent invasion by U. pinnatifida (Thomsen and South, 2019; Thomsen et al., 2019). Our study supports previous research that opportunistic and invasive species are usually outcompeted by native species, but adverse conditions can leave empty niche space and trigger their expansion (Valentine and Johnson, 2003, 2004; Occhipinti-Ambrogi, 2007; Epstein and Smale, 2017; South et al., 2017). Although it is difficult to predict the long-term effect of invasive species, more frequent MHWs might increase the probability that invasive species colonize and perhaps displace native species.

With one-time observation before and after, the ecological changes we describe cannot be confidently attributed to the MHWs. However, we chose our study sites on inhabited islands, since they are less impacted by cumulative human activities (Arafeh-Dalmau et al., 2017) than coastal areas, but are vulnerable ecosystems due to their isolation (Johnson and Black, 2006; Bell, 2008) making them ideal sites to understand kelp forests response to MHWs. Our study sites span a 250-km latitudinal gradient (almost half of the distribution of kelp in Baja California), allowing us to capture a good portion of the biogeographic variability of the region. There is also strong circumstantial evidence that the MHWs were, indeed, the driver of the observed changes: not only was there a corroborating signal in remote-sensed time series of kelp biomass, but we are unaware of any other plausible human-induced changes during the 2013–2018 period that could have altered the community structure in this way. Nevertheless, considering that the response of communities to the same disturbance might vary within the same region (Wernberg et al., 2013, 2016; Thomsen et al., 2019) and even within sites (Woodson et al., 2018), we acknowledge the need to understand the dynamics of these ecosystems at broader spatial and temporal scales and emphasize that our results should not be extrapolated without due care.

Kelp forest coverage has been reduced worldwide (Krumhansl et al., 2016; Wernberg et al., 2019), in part, due to climate change. Moreover, in some areas, their historical distribution limit has contracted poleward (Johnson et al., 2011; Wernberg et al., 2013, 2016), with clear signs of tropicalization of temperate ecosystems (Vergés et al., 2014). If positive temperature anomalies persist (2018 was still an anomalously warm year), the chronic effects of these warming events might shift Baja California’s ecosystems to a new alternative state, potentially precipitating range contractions for M. pyrifera and other species. Although future ecological states are difficult to predict in such dynamic ecosystems, our findings might be an early indication of a possible tropicalization of Baja California’s kelp forest (Cheung et al., 2013) under climate change, emphasizing the need to monitor and manage these ecosystems.

Southern kelp forests of the CCS are a shared marine resource between the south of California and the north of Baja California (Ramírez-Valdez et al., 2017). Advancing the understanding of the ecological impacts of MHWs (Hobday and Pecl, 2014) in this transboundary region to promote adequate management strategies will require a binational effort (Aburto-Oropeza et al., 2018). Thus, the protection of these valuable ecosystems needs immediate actions and close collaboration between southern California and Baja California (Arafeh-Dalmau et al., 2017). If we fail to act, once-abundant, extensive M. pyrfiera forests could experience important local declines, as they have in other regions (Johnson et al., 2011), leading to detrimental ecological and socio-political outcomes (Frölicher and Laufkötter, 2018).