Marine heatwaves drive recurrent mass mortalities in the Mediterranean Sea.

Climate change is causing an increase in the frequency and intensity of marine heatwaves (MHWs) and mass mortality events (MMEs) of marine organisms are one of their main ecological impacts. Here, we show that during the 2015-2019 period, the Mediterranean Sea has experienced exceptional thermal conditions resulting in the onset of five consecutive years of widespread MMEs across the basin. These MMEs affected thousands of kilometers of coastline from the surface to 45 m, across a range of marine habitats and taxa (50 taxa across 8 phyla). Significant relationships were found between the incidence of MMEs and the heat exposure associated with MHWs observed both at the surface and across depths. Our findings reveal that the Mediterranean Sea is experiencing an acceleration of the ecological impacts of MHWs which poses an unprecedented threat to its ecosystems' health and functioning. Overall, we show that increasing the resolution of empirical observation is critical to enhancing our ability to more effectively understand and manage the consequences of climate change.

INTRODUCTION

Anthropogenic climate change has become a major environmental driver, strongly affecting marine biodiversity and ecosystem functioning, often interacting with other human‐induced stressors (Gissi et al., 2021; Smale et al., 2019). Besides gradual ocean warming, acidification, deoxygenation, and sea‐level rise (IPCC, 2019), human activities have caused a rapid increase in the frequency and intensity of extreme climatic events, such as marine heatwaves (MHWs; Collins et al., 2019; Darmaraki, Somot, Sevault, & Nabat, 2019; Oliver et al., 2018). MHWs are discrete periods of prolonged anomalously warm water at a particular location (Hobday et al., 2016). In the last 20 years, MHWs have globally doubled in frequency and have become longer‐lasting, more intense and more extensive (Collins et al., 2019). This increase in the frequency of MHWs has been associated with the rapid increase in mass mortality events (MMEs) worldwide (Garrabou et al., 2009; Hughes et al., 2017; Thomson et al., 2015; Wernberg et al., 2016). Understanding the emergent properties and dynamics of these recurrent climatic disturbances is critical for predicting spatial refuges and cumulative ecological responses, and for managing long‐term impacts on ecosystems and the provision of services they provide (Hughes et al., 2021; Smith et al., 2021).

The Mediterranean Sea is generally considered a relevant model to assess the ecological effects of climate change on marine biodiversity and to test potential adaptation and mitigation strategies that may be scaled up (Cramer et al., 2018). Although the Mediterranean Sea represents only 0.32% of the total volume of the oceans, its unique geomorphological history led to a markedly high level of biodiversity with 7%–10% of all known marine species and a large proportion of endemic species (Bianchi & Morri, 2000; Coll et al., 2010). While barely recognized by society at large, marine MMEs are one of the most predominant biological impacts of climate change in the Mediterranean Sea (Cramer et al., 2018; Marbà et al., 2015; Rivetti et al., 2014). Indeed, Mediterranean MHWs have triggered unprecedented climate‐driven MMEs during the last decades and their occurrence is expected to increase in the coming decades (Darmaraki, Somot, Sevault, Nabat, Cabos‐Navaez, et al., 2019; Garrabou et al., 2019). The most dramatic reported events in terms of geographic extent and number of affected species occurred in 1999 and 2003 along the Northwestern Mediterranean Sea (Cerrano et al., 2000; Garrabou et al., 2009; Perez et al., 2000). These two events affected more than 40 species from various taxa (e.g., Porifera, Cnidaria, Bivalvia, Bryozoa, Ascidiacea) across thousands of kilometers of coastline.

The vast majority of studies on MMEs in the Mediterranean have focused on foundation species, such as habitat‐forming gorgonians and seaweeds (Chimienti et al., 2021; Garrabou et al., 2019; Verdura et al., 2021). Some species associated with MMEs have exhibited mortality rates of up to 80% of the studied populations (Cerrano et al., 2000; Garrabou et al., 2009) with major consequences for the structure and functioning of benthic ecosystems (Gómez‐Gras et al., 2021; Verdura et al., 2019). Yet, MMEs occurring within other Mediterranean populations and communities remains poorly understood, offering an incomplete picture of ecosystem‐wide impacts. Likewise, MMEs were mainly documented in the western Mediterranean Sea (Garrabou et al., 2019). However, given that the eastern basin is warming faster (Pisano et al., 2020), our knowledge may be biased from the extensive and long‐term sampling efforts in the western compared to the eastern basin. Indeed, the south‐eastern basin is the trailing edge for many native Mediterranean species, including molluscs and sea urchins, which have collapsed in the past decades (Albano et al., 2021; Rilov, 2016; Yeruham et al., 2015). Linking the timing and dynamics of these collapses to gradual ocean warming or to past, undocumented MHWs is hampered by a lack of historical baselines.

Since the observation of the 1999 and 2003 MMEs, different research teams have been monitoring climate‐driven mortalities across the Mediterranean basin. Given the apparent increase of mortality events and concomitant extreme heat conditions during the 2015–2019 period, we sought to investigate patterns and trends in MHWs and MMEs over temporal, spatial, and depth scales. Specifically, we characterize the exceptional Mediterranean wide thermal conditions recorded between 2015 and 2019 using high‐resolution, satellite‐derived sea surface temperature (SST) and subsurface coastal seawater temperature data (from 5 to 40 m depth). Then, we quantify the occurrence of MMEs during this period using observational data across the Mediterranean basin and explore potential relationship with MHW activity. By assessing and integrating sea surface and in situ temperature data with mass mortality records across the basin, our study provides the most up‐to‐date account of the impacts of extreme warming events on Mediterranean marine organisms and ecosystems.

MATERIAL AND METHODS

Characterizing 2015–2019 temperature conditions

Temperature data sets

Satellite‐derived SST across the Mediterranean Sea was obtained from CMEMS (https://resources.marine.copernicus.eu/?option=com_csw&view=details&product_id=SST_MED_SST_L4_REP_OBSERVATIONS_010_021). The data consist of daily (night‐time), gap free, optimally interpolated foundation SST at ~4 km resolution from AVHRR (Advanced Very‐High‐Resolution Radiometer) with improved accuracy and stability over the 1982–2019 period (Pisano et al., 2016, 2020). The global ecoregion analyses were performed with the NOAA OISST daily 25 km resolution product (Banzon et al., 2016, 2020; Reynolds et al., 2007). This optimally interpolated remotely sensed product is a gap free level 4 operational product with records from late 1981 to a 1‐ or 2‐day lag from the current date.

Subsurface coastal seawater temperature data were obtained from the T‐MEDNet temperature database. T‐MEDNet is a long‐term collaborative initiative devoted to building a pan‐Mediterranean observation network on climate change effects in marine coastal ecosystems (www.t‐mednet.org). The T‐MEDNet temperature database consists of high‐frequency (hourly) time series obtained using HOBO data loggers (accuracy ±0.21°C) set‐up at standard depths along rocky walls by divers, generally every 5 m from the surface to 40 m depth. To characterize the 2015–2019 thermal conditions across depths, from the 75 sites currently in the network we retained seven sites that had time series over the full depth range with a minimum length of 10 years of available data including the 2015–2019 period. The selected sites were located in the Western Mediterranean ecoregion (sensu Spalding et al., 2007), and within this area the sites where located in the European coasts, namely in the Columbretes islands Natural Park (Spain), El Montgrí, Les Illes Medes I El Baix Ter Natural Park (Spain), northern and southern areas of the Cap De Creus Natural Park (Spain), Calanques National Park (France), Port‐Cros National Park (France) and Scandola Natural Reserve (France). From the hourly data series, daily average seawater temperature was calculated. Missing data surrounded by shallower and deeper observations were interpolated vertically using linear interpolation. The obtained daily data series presented high temporal coverage for both the 2015–2019 period (97.0 ± 5.3%, mean ± SD) and 2010–2019 period (95.4 ± 6.6%, mean ± SD). Overall, the in situ seawater temperature observations span a range of hydroclimatic conditions and seasonal temperature stratification patterns (Figure S1; Bensoussan et al., 2010, 2019).

Temperature data analysis

The temperature analysis for each year focused on the warm hydrological period when thermal maxima were observed, defined hereafter as the period from June 1 to November 30 [JJASON]. This includes the period of maximum heat for the near surface (usually in July or August) and for the subsurface coastal ecosystems (after thermocline deepening in Fall, Figure S1).

Individual SST analyses were conducted per pixel. Statistics at the ecoregion and Mediterranean basin levels were obtained by averaging the per pixel statistics within the given region. Spatial SST analysis was conducted both at Mediterranean basin and ecoregion levels. For the ecoregion level, to cope with north–south heterogenous temperature conditions, the vast Western Mediterranean was subdivided in two at the 39°N parallel limit. Hence, spatial analysis of SST was conducted on a total of eight marine ecoregions as follows: Alboran, Northwestern Mediterranean, Adriatic, Ionian, Aegean, Levantine, Tunisian, and southwestern Mediterranean (Figure 1).

FIGURE 1

Patterns of warming and marine heatwaves (MHWs) across the Mediterranean Sea. (a) The difference in temperature between the mean of 2015 to 2019 minus the mean of 1982 to 1986. Green lines and blue dots have been added to visualize, respectively, the Mediterranean ecoregions and the monitoring areas with in‐situ (along depth) temperature monitoring (also in panel b). (b) The highest category of MHW experienced from 2015 to 2019 at locations where only a category I MHW was experienced in 1982 to 1986. Light pink areas show when a category II or greater MHW occurred in the earlier period. (c) Annual sea surface temperature (SST) anomalies from 1982 to 2019 (climatology period is 1982 to 2019). Horizontal black bars show 5‐year averages over the time series. (d) The annual surface area (%) of the Mediterranean affected by category II or greater MHWs from 1982 to 2019. The color of the bar shows what proportion relates to the MHWs category. Horizontal black bars show roughly 5‐year averages in % of surface affected by category II or greater MHWs; finally the black dots show the average global ocean cover (%) per year for category II or greater MHWs. Category I moderate events are excluded here to allow for a better understanding of the increasing severity of MHWs. (e) The average per depth strata decadal warming trend for 7 nearshore sampling areas from the Northwestern Mediterranean sea. Error bars indicate the standard deviation. The asterisk (*) indicate that satellite data were used for the surface whereas subsurface data originate from in situ monitoring every 5 m from 5 to 40 m depth. (f) The average along depth number of MHW days for the months from June to November over the study period (2015–2019) for the seven Northwestern Mediterranean sampling areas. Error bars indicate the standard error. The asterisk (*) indicate that satellite data were used for the surface, whereas subsurface data originate from in situ monitoring every 5 m from 5 to 40 m depth. The map lines in panels (a) and (b) delineate study areas and do not necessarily depict accepted national boundaries.

The seasonally varying climatology for each SST pixel was determined by first finding the average daily temperatures from 1982 to 2011 and applying a moving average with an 11‐day window to them. A second pass of smoothing was then given using a 31‐day window. The same method was used to create the 90th percentile thresholds that were used to detect MHWs.

MHWs, and their categories, were detected in this study using the standard methodology outlined in Hobday et al. (2016, 2018) and implemented in R (R Core Team, 2021; Schlegel & Smit, 2018). Whenever SST values exceeded the 90th percentile threshold determined for each pixel for five or more continuous days over the JJASON period per year this was flagged as an MHW. The duration of a MHW is the count of the number of days from the beginning to the end of the event. The annual duration values were determined by adding up the duration of all MHWs over the JJASON period. The annual chronic heat stress of organisms during the JJASON period was quantified as the number of MHW days a site was exposed to. To provide a general quantitative understanding of the acute heat stress experienced by many different taxa the categories of MHWs were used. These are determined by the number of times the temperature anomaly is in excess of the distance between the seasonally expected temperature and the 90th percentile threshold. A weak positive anomaly will be flagged as a category I “Moderate” event. Whereas an anomaly twice as high as the aforementioned distance will be II “Strong,” three times is III “Severe,” and four times or greater is IV “Extreme.” The MHW duration and category values used for the global comparisons were extracted from the MHW Tracker (http://www.marineheatwaves.org/tracker.html; Schlegel 2020; Hobday et al., 2016, 2018).

To determine how rapidly the surface of the Mediterranean Sea may be warming, the average SST values from the first 5 years of the satellite period (1982–1986) were subtracted from the average values over the study period (2015–2019). These differences were calculated for the Mediterranean Basin, per ecoregion, and per pixel. Similar 5‐year average comparisons were performed for the annual average MHW days detected. Because the presence of moderate MHWs are relatively common, both in the Mediterranean basin and globally, a metric for the change in the category of MHWs over time required more nuance. It was decided that because category I “Moderate” MHWs are consistently present throughout the satellite record, they would not be included in comparisons of the first and last 5‐year periods. To determine how much the occurrence of MHWs increased over the surface of the Mediterranean Sea, the highest category event per pixel in the first 5‐year period (1982–1986) was checked against the highest category event in the current study period (2015–2019). If the highest category of MHW in the study period was not in excess of the highest category of MHW in the first 5‐year period, then it was considered that the intensity of MHWs on that pixel had not increased. There were no pixels in which the category in the first 5‐year period was greater than the current period. The MHW days and categories experienced per year (not only the JJASON period) in the Mediterranean were compared against the global averages by extracting these statistics from the MHW Tracker.

Considering the strong vertical temperature gradients across a few tens of meters during the seasonal stratification in the nearshore monitoring areas (Figure S1), we further investigated how the warming conditions at the sea surface are being transmitted to the subsurface layers taking advantage of high‐resolution in situ temperature series for decadal trends and MHW analysis. From the seven sites retained for analysis, the co‐located satellite SST data were retrieved to allow their comparison with the in situ data.

The decadal [2010-2019] warming trend was calculated from the linear regression of monthly anomalies versus time. As for the SST, MHWs were detected using the standard methodology outlined in Hobday et al. (2016) and implemented in Matlab (Zhao & Marin, 2019) using the 2010–2019 period for climatology and thresholds calculation. As above, but focusing on the years 2015–2019, the annual duration values were determined by summing the duration of all MHW days over the JJASON period. To determine the warming and MHW patterns across the vertical dimension, the statistics for the seven sites were pooled for the calculation of depth average and standard error values.

Assessment of MMEs during 2015–2019 across the Mediterranean

Field surveys

Data were obtained through benthic surveys conducted by 33 research teams from 11 Mediterranean countries. Surveys covered thousands of kilometers of coastline, spanning 13° of latitude (32°S to 45°N) and 40° of longitude (−5° W to 35° E) in the Mediterranean Sea. The data set represents the most comprehensive inventory of MME records for benthic species between 2015 and 2019 in the region. This database, as with the depth temperature monitoring, were undertaken as part of the T‐MEDNet initiative (www.t‐mednet.org). The surveys were conducted in 142 monitoring areas. Monitoring areas were considered as geographic areas (10–25 km coastline, for example, a marine protected area and the nearby coast) sharing common environmental features. The data set provides quantitative information to assess the patterns of MMEs regarding the following: (i) species and habitats; (ii) geographic and temporal patterns; and (iii) depth distribution (0–45 m depth; Garrabou et al., 2022).

Quantification of MMEs

Mass mortality records were obtained from quantitative benthic surveys assessing the mortality impact of local populations at specific times. Here, we consider a local population as a group of colonies/individuals or cover of the same species (ranging from tens to thousands colonies/individuals depending on the species) dwelling in a specific geographic location defined by spatial coordinates and depth range. Note that from the same geographic location we can obtain several mass mortality records obtained from surveys on different species and/or for the same species at different depth ranges and years.

The surveys quantified the percentage of affected colonies/individuals for each species either directly in the field or from photographic/video sampling. In this study, we considered affected colonies/individuals displaying signs of recent mortality such as denuded skeletons in gorgonians, sponges, and scleractinian corals and empty valves in bivalves attached to the substratum. The observations on mortality across taxonomic groups may depend on the timing of the surveys and the duration of the specimens' structures (e.g., soft vs. hard skeletons) attached to the substratum after their death. In the case of clonal species (e.g., gorgonians, bryozoans, sponges; Figure 2; Video S1), we considered an individual/colony to be affected by mortality when it showed recent tissue necrosis over 10% of its surface (Garrabou et al., 2009, 2019). In the case of Posidonia oceanica and habitat forming macroalgal species, increase in shoot mortality or a sharp decline on thallus densities, respectively, served as indicators (Marbà & Duarte, 2010; Verdura et al., 2021). In the case of mussel beds, decline in cover was considered. Depending on the abundance of the surveyed species/taxa in the monitored locations, the total number of surveyed colonies/individuals ranged from 10 for the less abundant (e.g., sponges) to hundreds of colonies/individuals for the most abundant (e.g., gorgonians), or in the case of mussels to 100% cover of the rock.

FIGURE 2

Ecological impacts of the 2015–2019 mass mortality events (MMEs). (a) A red gorgonian Paramuricea clavata colony exhibiting partial mortality (brown areas overgrown by epibionts) contrasting with live tissue (red areas; 35 m depth in Tavolara Island, Sardinia, Italy; photo by: M. Munaretto). (b) A colony of the reef‐building coral Cladocora caespitosa with signs of partial mortality (bare skeleton in white areas) contrasting with live tissue (brown areas; 3–4 m depth in Tabarca Island, Alicante, Spain; photo by A. Izquierdo). (c) A colony of the bryozoan Myriapora truncata showing mortality (white parts overgrown by epibionts; 15 m depth in Ibiza, Balearic Islands, Spain; photo by J. Garrabou), (d) An individual of the mollusk Spondylus gaederopus affected by mortality (white; 25 m depth in Ullastres, Catalan coast; photo by P. López‐Sendino). (e) The massive sponge Agelas oroides exhibiting mortality (brown spongin network free of cellular material) contrasting with the alive orange tissue (16 m depth, Chalkidiki peninsula, North Greece; photo by V. Gerovasileiou). (f) Almost completely dead individual of Spongia officinalis showing brown spongin network free of cellular material (10 m depth in Portofino, Liguria, Italy; photo by F. Betti). (g) Partially dead rhodophyte Mesophyllum sp. (10 m depth, Kalymnos Island, South Greece; photo by M. Sini). (h) Arbacia lixula exhibiting mortality with the loss of spines (3 m depth in Cap de Creus, NE Spain; photo by J. Garrabou). (i) Percentage (%) of MME records per phylum across the entire data set. (j) Number of taxa affected per phylum over depth. (k) Percentage (%) of MME records per phylum over depth. (l) Number of affected taxa per phylum and habitat category (barplot), and percentage (%) of affected habitat category across the overall mortality data set (pie chart).

MME severity

Four severity classes were considered based on the values of the percentage of affected colonies/individuals reported in each survey: No mortality (<10%), Low (≥10, <30%), Moderate (≥30, <60%), and Severe (≥60%) MME. The same categories were applied for species for which we monitored decline of cover (e.g., mussels).

Species and habitats affected

MMEs was assessed by field surveys mostly on species already affected in previous climate‐driven mortality events in the Mediterranean Sea (Cerrano et al., 2000; Garrabou et al., 2009, 2019). These species are mainly easily identifiable macrobenthic species, and often abundant at the surveyed locations. The majority of surveys identified benthic taxa (n = 50) to the species (90%) and genus levels (6%) while 4% were identified at level of broad taxonomic group (e.g., calcareous algae; Table S1).

Each MME record was assigned to three broad habitat categories, namely hard substrates, mobile substrates, and seagrasses, which represent the three principal aspects of the coastal Mediterranean seascape. Hard substrates refer to all types of natural or artificial rigid surfaces, mainly including rocks of varying physiognomic characteristics in terms of geological origin, inclination, and biotic features. Based on their overall structure and degree of exposure to natural or human disturbances, hard substrates were further subdivided into three distinct habitat subcategories: (a) the upper sublittoral zone, which is overall more exposed to light, thermal changes, wave action, and anthropogenic disturbances (for simplicity reasons in this subcategory we included the medio‐littoral zone since there were very few observations), (b) the lower sublittoral zone, which is generally characterized by more stable physical conditions, and within the Mediterranean basin it is mainly associated with the pre‐coralligenous and coralligenous assemblages (Ballesteros, 2006), and (c) marine caves, including both semi‐ and entirely submerged sea caves (Gerovasileiou & Bianchi, 2021). Mobile substrates include types ranging from mud, well sorted fine sand and gravel, to rhodolith and maërl beds, while the habitat category of seagrasses relates to observations of the endemic species P. oceanica.

Spatial and temporal mortality impact patterns

The surveys were mainly conducted in summer and fall seasons when the effects of climate‐driven mortality events usually occur (Garrabou et al., 2009). The surveys were conducted in 142 Monitoring areas (see Section 2.2.1 for the definition) encompassing a total of 352 surveyed locations (minimum of one and a maximum of 15 survey locations per Monitoring area and year), resulting in a total of 985 survey records. Observations covered all the eight Mediterranean ecoregions considered in this study (see above).

Depth related patterns of MMEs

The depth distribution of mortality occurrences was analyzed based on the lower and upper depths of mortality reported in each single record and assuming continuity of the phenomenon in this depth interval. The vertical profile was divided into levels of 5 m each, from the surface to the maximum recorded depth, and for each depth level the number of occurring records belonging to the different species and different severity classes (low, moderate, and severe) was counted. For instance, a record reporting Paramuricea clavata damaged from 0 to 16 m, was included in the counting in each of the following depth levels: 0–5, 5–10, 10–15, and 15–20 m. Records were analyzed year by year and for the entire study period. Species data were aggregated by phyla. In the plots, data were reported in terms of percentage of the total mortality records in the considered period. The main information gathered for each MME record is provided in Table S2.

Relationships between MHWs and MMEs

To explore to what extent the onset of the MMEs observed during the 2015–2019 were related with heat exposure (using both SST and in situ temperature data), we focused on the Northwestern Mediterranean since this ecoregion has most of the mortality records as well as the most comprehensive monitoring of both the biology and temperature series (65% mortality records and 100% of the monitoring areas with at least 10 years of in situ temperature time series) allowing for comprehensive coverage across years and depths.

We conducted a two‐level analysis. At the first level, we explored if the temporal patterns in surface heat exposure observed in the basin generally matched the patterns of ecological impacts observed during the 2015–2019 period. For this analysis, we first obtained the annual regional average of MHWs days calculated over the JJASON period with SST data per year corresponding to pixels (or the closest ones) in which mortality surveys were conducted (see above); and then jointly plotted the resulting MHW patterns with the annual mortality incidence at the ecoregion level by calculating the proportion of surveyed records displaying mass mortality impacts (including all different categories of impact Low, Moderate, or Severe).

For the second level, and to further explore the potential cause–effect relationship between heat exposure conditions and mortality incidence, we fitted two generalized linear models (GLMs) with binomial error structure and with either SST or in situ temperature data as the independent variables. In the first model, we used per pixel MHW days obtained from SST as the explanatory variable and the mass mortality incidence obtained at the Monitoring areas level as the dependent variable. In this case, MHWs days over the JJASON period were averaged across all pixels covering each Monitoring area (usually 1 to 10 pixels, but more than 20 for Cap de Creus and El Montgrí). For this analysis, only Monitoring areas with at least 2 years of monitoring over the 5‐year study period were selected, and within these areas, only years with at least three or more mortality records of recurrently monitored populations were considered to quantify the mortality incidence. These criteria allowed us to perform the analysis using the most consistent data in which mortality patterns could be unequivocally dated within the study period based on strong local ecological knowledge from the different research teams. This is critical to ensure an accurate assessment of the relationship between ecological data and the corresponding MHW patterns. Moreover, because the independent variable was surface temperature, we only included the shallowest mass mortality records in the analysis when multiple records were reported for the same species at the same location but at different depths. Finally, to avoid under estimation of MMEs, records from species that were not affected at any time during the study period across the region (e.g., Paracentrotus lividus) were not included in the analysis. Finally, the GLM was weighted by the sampling size in each year and monitored area to account for potential sampling biases.

A second GLM with in situ (across depth) data was conducted on monitored areas with both long‐term, in situ temperature monitoring (at least 10 years of in situ temperature data) and robust ecological data (at least three survey records for a given depth and year). Considering these criteria, only three monitored areas from the Catalan Sea could be included in the analyses. Within these areas, the number of in situ MHW days per year at each depth was considered as the explanatory variable and the corresponding mass mortality incidence at different depths was set as dependent variable. As for the first GLM, weights were also used to account for potential sampling biases across the three monitored areas, years and depths. The data sets used for the analysis performed are available at Garrabou et al. (2022).

RESULTS AND DISCUSSION

2015–2019: The warmest years on record in the Mediterranean Sea

Every location in the Mediterranean Sea has experienced considerable increases in annual mean temperatures in recent decades (Figure 1a). Over the entire Mediterranean basin, the 2015–2019 period has been the warmest since the beginning of satellite records in 1982 (Figure 1c), with average temperature now 1.2 ± 0.23°C (mean ± SD of absolute annual anomalies) higher compared with the 1982–1986 period. In fact, the average Mediterranean SST for three of the 5 years in the study period (2015, 2016, and 2018) were higher than in 2003, when the most damaging MHW on record occurred in the Mediterranean (Garrabou et al., 2009). However, different warming rates were observed across the Mediterranean ecoregions. In the Aegean and Levantine Sea ecoregions the average temperatures over the study period (2015 to 2019) were 1.4 ± 0.25°C and 1.3 ± 0.23°C warmer than the first 5 years of the satellite record (1982 to 1986; Figure 1a). Although the Mediterranean Sea is not warming as rapidly everywhere, the lowest average increase in temperature between the first and last 5‐year periods (Alboran Sea at 1.0 ± 0.32°C) is still greater than the global average. The average warming rate in the Mediterranean Sea from 1982 to 2019 was 0.38°C per decade, more than three times higher than the global average of 0.11°C per decade (IPCC, 2019).

Considering that SST in the Mediterranean Sea is increasing so rapidly, it is not surprising that there has been a dramatic rise in the occurrence and severity of MHWs. Nearly the whole Mediterranean basin experienced much more intense events during the 2015–2019 period than during the first 5 years of the satellite record (1982 to 1986; Figure 1b). During the first 5 years considered (1982–1986), there were very few MHWs with a category greater than I “Moderate,” (10% for II “Strong,” 0.05% for III “Severe,” and 0% for IV “Extreme”), whereas during 2015–2019, 99.99% of the surface of the basin was affected by at least category II “Strong” events, 55% by category III “Severe” events, and while they remained rare at only 3%, IV “Extreme” MHWs occurred on an annual basis (Figure 1d; MHW Category I “Moderate” not shown). For the full satellite record of 1982 to 2019, there have been only 11 years when extreme MHWs were detected, with five of these being the years of our study period (2015–2019). As with the differing rates of warming, Mediterranean ecoregions have also experienced different increases to the categories of MHWs detected therein. The Aegean Sea was one of the most impacted regions by MHWs, but it was also the region that experienced the most intense MHWs during the first 5‐year period (1982 to 1986), making the overall upward trend in MHWs here less steep than other ecoregions. The region that has seen the most dramatic increase in the category of MHWs has been the Levantine Sea, with the Tunisian Plateau the lowest increase (Figure 1b).

The annual increase in strong MHWs at the surface is also well supported by our results (Figure 1d). Less than 25% of the surface of the basin was affected by strong MHWs for most of the first half of the satellite record, with 2003 being a notable exception. Since 2012, the Mediterranean Sea has been experiencing strong MHWs over more than 75% of its surface, with 3 of the 5 years in the study period (2015, 2016, and 2018) having more than 90% coverage by strong or greater events. The spatial coverage of MHWs in the Mediterranean Sea since roughly 2010 has been higher than the global average and is increasing more rapidly (Figure 1d; Blunden & Boyer, 2020).

Changes in temperature at depth

The analysis, from the surface (derived from SST) and from 5 to 40 m of depth, for seven Northwestern Mediterranean locations over the 2010–2019 period (Figure 1a) evidenced that the rapid warming was not only for the surface, but extends down to the intermediate (20–25 m) and deep (>40 m) coastal habitats studied (Figure 1e). Warming was found highest below the surface, at 10 m depth (0.092 ± 0.019°C per year). It decreased with depth to 0.056°C per year at 40 m depth, which still remained 47% stronger than the average 1982–2019 Mediterranean SST warming taken as reference (see previous section) in line with previous studies (Kersting et al., 2013; Salat et al., 2019).

The occurrence of subsurface MHWs during the 2015–2019 period was also quite clear (Figure 1f). For these 5 years, a total of 73 MHW days were obtained at the surface, decreasing along the depth gradient to 47 ± 25 days at 20 m and to 32 days at 40 m.

2015–2019 MMEs: Species, habitats, spatiotemporal scales, and depth gradients

From a total of 985 field surveys conducted over the 2015–2019 period, 567 (58%) returned evidence of MMEs. These MMEs encompassed 50 taxa from 8 phyla (Figure 2; Table S1; Video S1), impacted different types of habitats (from hard to mobile substrates and seagrasses; Figure 2l) and were widespread across depths (0 to 45 m), ecoregions and years (Figures 2j,k and 3). Never before has the Mediterranean basin been affected by MMEs with such a number and diversity of affected taxa at such large spatial scales over five consecutive years (Garrabou et al., 2019).

FIGURE 3

Spatio‐temporal patterns of mass mortality event (MME) records and their severity across depths over the 2015–2019 period (a) Spatio‐temporal trends of MME records. For those areas that have been monitored from 1 to 3 years, the size of the circles represents the number of monitoring years, and their transparency is related to a mortality ratio on those areas. The mortality ratio (from 0; minimum to 1; maximum) is calculated by dividing the number of years exhibiting mortality in one area by the total number of monitoring years in that area. In the case of a mortality ratio equal to zero, circles have been colored in blue. For an extended version of this panel showing detailed temporal trends in areas with 1–3 years of monitoring, see Figure S2. For areas with more than 3 years of monitoring, pie charts have been used to represent the temporal trends in each of the 4/5 years of the study. (b) Severity of MME records across the entire data set. (c) Severity of MME records across depth. (d–h) Temporal trends in the severity of MME records across depth. The map lines in panel (a) delineate study areas and do not necessarily depict accepted national boundaries.

The number of taxa and phyla exhibiting mass mortality per year during the 2015–2019 period was on average 23 taxa and 7 phyla (ranging from 21 to 28 taxa and 7 to 8 phyla, respectively depending on the year). These values are much higher than reports for most previous years from 1978 to 2014 (Garrabou et al., 2019). In fact, only in 1999 (with 33 taxa from 6 phyla affected) and 2003 (with 26 taxa from 4 phyla affected) were similar or higher values exhibited (Cerrano et al., 2000; Garrabou et al., 2009; Perez et al., 2000). However, not all taxa and phyla were equally affected during the 2015–2019 MMEs (Figure 2i–k). Cnidarians were much more affected than any other group, with 14 impacted taxa and a number of observations accounting for more than 54% of the overall MME records. After Cnidaria, Bryozoa (8 taxa and 10.6% of MME records) was the second most affected taxonomic group, followed by Rhodophyta (3 taxa and 8.3% of MME records). These observed differences across taxonomic groups are in line with previous observations from local/regional MMEs (Garrabou et al., 2019), and are likely the consequence of the existence of inter and intra‐specific differences in tolerance to heat stress (Gómez‐Gras et al., 2019; Ledoux et al., 2015; Pagès‐Escolà et al., 2018; Savva et al., 2018). Interestingly, mortality was far from being restricted to shallow waters. Indeed, although the number of affected taxa gradually decreased with depth (from shallow waters down to >25 m depth with 30 to 8 taxa, respectively, Figure 2j,k), about a third of the overall MME records were reported from depths greater than 15–20 m, with 15–25 m being the most affected depth range. The depth range affected in the 2015–2019 MMEs was therefore similar to previously reported events from the Mediterranean (Garrabou et al., 2019). It is difficult to determine if this depth pattern of mortality is unique to the Mediterranean Sea considering the substantial gaps in knowledge in other temperate and tropical ocean regions, where data are mainly restricted to intertidal and shallow waters (e.g., Wernberg et al., 2016; Hughes et al., 2017, but see Frade et al., 2018).

As well as the wide range of species, phyla and depths affected, all the main Mediterranean coastal habitats experienced MMEs (Figure 2l). The majority of MMEs were associated with hard substrates, with 44 affected taxa and 88% of MME records, while seagrasses (6 taxa and 10% of MME records) and mobile substrates (6 taxa and 2% of MME records) were also severely impacted. Notably, the majority of MMEs associated with hard substrates occurred in shallow zones (0–25 m depth; 69%) and involved “coralligenous” animal‐dominated assemblages (49%), one of the high‐diversity communities in the Mediterranean (Ballesteros, 2006). A similar range of impacted habitats was reported during the 1999 and 2003 MMEs in the Mediterranean Sea (Garrabou et al., 2019), whereas reports of climate‐related MMEs from other ocean regions (temperate and tropical) usually involve just one or two broad types of habitat (e.g., coral reefs, kelp forests and/or seagrasses; Smith et al., 2021).

Overall the observed patterns point toward a future of winners and losers in the Mediterranean Sea in the face of climate change, which could likely lead to local, regional, or even pan‐Mediterranean ecological extinctions of species, widespread structural and compositional changes of ecological communities, and potential subsequent changes in ecosystem functioning, particularly where lost species are functionally unique (Bellwood et al., 2004; Bianchi et al., 2014; Harvey et al., 2022; Loya et al., 2001; Moullec et al., 2019). This could be the case for instance of the Mediterranean gorgonians (e.g., P. clavata, Corallium rubrum, Eunicella cavolini, and Eunicella singularis) which are among the most affected organisms in our study, comprising about 30% of MME records. These gorgonians, which dominate many diverse and abundant rocky communities (i.e., coralligenous assemblages), are considered to be functionally unique in the Mediterranean Sea as they provide the highly structurally complex, 3D habitats that are needed for many other associated species to thrive (Gómez‐Gras et al., 2021; Ponti et al., 2014, 2018; Verdura et al., 2019). This result suggests that Mediterranean gorgonians will likely suffer further declines in the upper range of their bathymetric distribution as climate change unfolds, which could have significant consequences for the functioning of Mediterranean benthic ecosystems (Gómez‐Gras et al., 2021) and subsequently, on the provision of associated services to human societies (Smith et al., 2021).

Our results also revealed several clear spatiotemporal patterns, indicating the pervasiveness of MMEs. In particular, from a total of 142 monitored areas covering all Mediterranean ecoregions and located tens to thousands of kilometers apart, the vast majority (78%) exhibited mass mortality of some taxa in at least one of the surveyed years (Figure 3a; Table S3). The severity of MMEs showed a geographic gradient, being highest in the eastern ecoregions (Levantine, Aegean, and Tunisian plateau) with more than 90% of mortality records corresponding to severe impact, and lowest in the Alboran Sea, where more than 50% of records correspond to low impact. Meanwhile, the records from the western and central ecoregions showed similar percentages (i.e., 33%) across the three severity classes (Figure 3a).

Regarding the temporal patterns, no single year out of the five studied had less than four ecoregions affected, with the year 2019 having the most widespread impact across seven ecoregions (Figure S2). At the ecoregion level, those located in the western part of the Mediterranean displayed impacts during all 5 years, followed by the Alboran and Levantine seas, which suffered MMEs during 4 years. The remaining ecoregions were impacted 2–3 years by MMEs (Figure 3a; Figure S2). Focusing at Monitoring areas level, from the 24 monitoring areas that were surveyed at least 3 years along the 5‐year study period, 79% suffered recurrent yearly impacts, whereas 42% of them experienced at least 4 years of cumulative MME records (Figure 3a).

Overall MME severity classes had a similar percentage of mortality records during the first 3 years. In 2018, the severity was lower, since almost 50% of MMEs were associated with low impacts. Inversely, 2019 had the highest impact, with 50% of records showing severe impact (Figure 3b). Patterns of severity along depth, in general, followed the mentioned yearly patterns, but with an observed decrease in severity with depth. Severe impacts were mainly observed in shallow and intermediate waters, while at deeper waters (>25 m), damage was mainly low or moderate (Figure 3d–h). However, the percentage of records defined as severe increased at all depth ranges in 2019. This vertical distribution of MME records may be mediated not only by the heat exposure (or driving factors) but also by the bathymetric range distribution and susceptibility of the affected species in the different geographic areas. Despite the observed variability in time and depth, it is important to underline that severe MMEs were recorded every year at all investigated depths; moreover a general trend of worsening of severities occurred at all depths throughout the study period (Figure 3d–h).

Because of the collaborative approach adopted, the sampling effort was not homogenous across monitoring areas, ecoregions and through the study period (Figure 3a; Table S3). For instance, Northwestern Mediterranean has 36‐times more monitoring areas than the Ionian ecoregion and within the Northwestern Mediterranean, only 42% of monitoring areas were sampled 3 years or more. However, considering the basin‐wide ocean warming trends (Figure 1), we suggest that increased sampling resolution across the different dimensions (spatial, temporal and taxonomic) would have only resulted in similar or even higher MME occurrences. That is, it is likely that many MMEs remain undetected, particularly in ecoregions and areas with lower research effort.

Overall, the observed recurrent sequence of five consecutive years with large‐scale, MMEs affecting multiple species is unprecedented in the Mediterranean basin and globally (Garrabou et al., 2019; Smith et al., 2021). During the last two decades, apart from the widespread MMEs reported in this study, the Mediterranean suffered only two major MMEs (1999 and 2003). Likewise, although recurrent mass bleaching and mass mortalities of corals and kelp forests have been observed along thousands of kilometers in the Atlantic, Pacific and Indian oceans (Hughes et al., 2017; Muñiz‐Castillo et al., 2019; Wernberg et al., 2016), only coral reefs have suffered a similar succession of recurrent large‐scale coral bleaching events during the 2014–2020, but not in five consecutive years period (Cheung et al., 2021; Hughes et al., 2017, 2021). Our results therefore highlight a potential new feature of the Mediterranean Sea as a climate change hotspot: the onset of recurrent large‐scale MMEs.

In the Northwestern Mediterranean Sea between 40% and 75% of surveyed populations were affected yearly by mortality during the 2015–2019 period (Figure 4b). This extensive onset of mortality events was concomitant to the extreme heat exposure observed during that period (Figure 4b). Indeed, the mortality incidence (i.e., proportion of surveyed populations affected by mortality) tracks well the annual heat exposure over the study period with 2016 and 2018 showing the lowest (about 40% of affected populations and 50 MHW days) and the highest (about 75% affected populations and 80 MHW days) values of mortality incidence and MHW days, respectively (Figure 4a,b). This result is in line with previous studies showing that MHWs are associated with the onset of MMEs in the Mediterranean Sea (Cerrano et al., 2000; Garrabou et al., 2009; Perez et al., 2000) as well as in other ocean regions (e.g., Hughes et al., 2017; Smale & Wernberg, 2013; Wernberg et al., 2016). Beyond this correlative evidence, establishing the relationship between the biological responses of marine biodiversity and the different levels of heat exposure is challenging and has been rarely addressed (Cheung et al., 2021; Hughes et al., 2021). Indeed, the high variability of responses observed among different species and populations and across spatial and temporal scales, as well as the lack of empirical data sets on these extreme (rare) events has undermined our ability to understand this relationship.

FIGURE 4

Relationship between heat exposure (marine heatwave [MHW] days) and mortality incidence in the Northwestern Mediterranean ecoregion during 2015–2019. (a) Map of the Northwestern Mediterranean ecoregion showing the location of the monitored areas included in the analysis; black dots indicate all monitored areas in the Northwestern Mediterranean Sea (used in the regional analysis shown in panel b), yellow dots indicate only the monitoring areas considered for the analysis shown in panel c and red triangles indicate the areas with long‐term, in situ temperature monitoring used for the in‐depth analysis shown in panel d (see Section 2.3 for further details). (b) Bars and points show, respectively, the yearly mean number of MHW days and mortality incidence (proportion of records showing mortality) observed across the Northwestern Mediterranean basin. Panels c and d show respectively the relationship between heat exposure (yearly average of MHW days during the JJASON period) at the surface (sea surface temperature) or across depth (from 5 to 40 m) and the corresponding mortality incidence in the studied monitored areas, years, and/or depths. The lines show the predicted values of the generalized linear models and their confidence interval (95%). The size of the points is proportional to the sampling size (number of populations for all the species surveyed in each monitored area, year and/or depth; minimum number = 3, maximum number = 50). The map lines in panel (a) delineate study areas and do not necessarily depict accepted national boundaries.

Here, we demonstrate that the mortality incidence is significantly related to the MHW days reported both from the satellite sea surface (Figure 4c; p value <.001; Table S4) and in situ temperature data along a depth gradient (Figure 4d; p value <.05; Table S4). Nonetheless, one can observe a high variability in the level of mortality impact for similar heat exposure. Considering for instance the range from 25 to 50 MHW days per year based on SST, the mortality incidence varies in the Monitored areas between 20% and 100% (Figure 4c). Similar variable responses were observed when MHW days based on in situ temperature data were considered (Figure 4d). The main components explaining this variability could be linked to (i) properties of the MHWs such as timing, duration, maximum and cumulative intensity (Elzahaby et al., 2021; Hobday et al., 2016), (ii) taxonomy (differences among species and populations, including their specific thermal niche, thermotolerance and physiological status; e.g., Coma et al., 2009; Arizmendi‐Mejía et al., 2015; Crisci et al., 2017; Pagès‐Escolà et al., 2018; Gómez‐Gras et al., 2019) and (iii) the ecological memory on the recurrence of MHWs over the same geographic areas (Figure 3), which may alter the response of populations re‐exposed to MHWs (Hughes et al., 2021; Turner et al., 2020). In fact, disentangling the factors behind this relationship including also the complex eco‐evolutionary feedbacks between the changing thermal environment and marine biodiversity is urgently needed to anticipate the consequences for ecosystems and human well‐being (Fuller et al., 2020; Hughes et al., 2021; IPCC, 2019; Turner et al., 2020).

Our results also show that the MHW days metric based on in situ temperatures may be problematic when comparing against MMEs because 0 MHW days often correspond to high mortality incidence at different depths (Figure 4d). This mismatch however is almost certainly due to the underestimation of heat stress conditions when using MHW in situ data. Firstly, because Mediterranean coastal waters are characterized by dynamic stratification regimes with short term (daily) large temperature changes (e.g., 10°C, Figure S1), which may prevent thermal conditions to overcome the 90th percentile temperature threshold for five or more continuous days; the criteria used for estimating MHW days (Bensoussan et al., 2010, 2019; Hobday et al., 2016). Secondly, because only the most recent 10 years of in situ data were available for the creation of the 90th percentile temperature thresholds (~2010–2019), which, given the current warming trend, implies that the temperature thresholds for the detection of in situ MHWs are much higher than those used for the surface data based on a longer climatology period (e.g., 1982–2011; Schlegel et al., 2019). Overall, this means that many more in situ MHW days could have likely been detected if a longer historical record of in situ temperature data would have been available. Similarly, many days with extreme hot conditions may have been omitted because of not meeting the 5 consecutive days criteria (i.e., heat spikes sensu Hobday et al., 2016). Thus, while MHW days are generally considered representative for chronic heat stress exposure (e.g., Oliver et al., 2018; Smale et al., 2019), the estimation of extreme hot days could also be desirable when attempting to quantify acute heat stress, especially along depth. Indeed, both approaches are relevant from a geophysical perspective and document different facets of the ongoing environmental changes. Defining more integrative heat stress indicators accounting for both exposure duration and intensity would also likely increase our predictive capacity on ecological impacts (Cheung et al., 2021; Hughes et al., 2021). Finally, to enhance our predictive capacity in the Mediterranean coastal areas, we need to develop our knowledge of the physical forcing of MHWs at the surface as well as at depth (Elzahaby et al., 2021).

CONCLUSIONS

The increase of frequency, intensity, and spatial scales of MHWs is driving major ecological changes in the marine ecosystems worldwide (Smith et al., 2021). Our results clearly indicate that the Mediterranean Sea is experiencing an acceleration of climate change impacts but the implications for socioecological systems are poorly understood. In fact, our ability to predict how MHWs will affect marine species and habitats is undermined by key gaps of knowledge. For instance, improving our understanding of the interactions among multiple stressors, the role of ecological memory and the importance of species' functional traits is urgently needed to identify the vulnerability of biodiversity to MHWs (Harvey et al., 2022; Turner et al., 2020). To efficiently tackle any of these challenges, increasing the resolution of empirical observation capacities over relevant spatial, temporal and taxonomic scales is essential. As we showed here, supporting international, collaborative research schemes offers an effective approach to this goal. As marine ecosystems are entering uncharted territories, reinforcing coordination and cooperation at regional, national, and international levels will undoubtedly help the endeavor of developing more effective management decisions in the face of ongoing climate change. Finally, our findings along with the two consecutive back to back bleaching events also observed in the Great Barrier Reef (Cheung et al., 2021; Hughes et al., 2021) might be worrisome signals that large‐scale mass mortality outbreaks are no longer the exception, but the new normal.

FUNDING INFORMATION

This paper was supported by Euromarine. Joaquim Garrabou acknowledges the funding by the “Severo Ochoa Centre of Excellence” (CEX2019‐000928‐S), the MCIU/AEI/FEDER [HEATMED; RTI2018‐095346‐B‐I00], Interreg‐Med Programme MPA‐Engage (1MED15_3.2_M2_337), the European Union Horizon 2020 research and innovation programme (Futuremares SEP‐210597628). Nuria Teixido acknowledges the French National Research Agency (4Oceans‐MOPGA grant, ANR‐17‐MPGA‐0001) and internal funds from the Stazione Zoologica Anton Dohrn. Gil Rilov was supported by the Israel Ministry of Environmental Protection. Hocein Bazairi, Free Spinosa, and Vasilis Gerovasileiou acknowledge the funding by the MAVA Fondation (MedKeyHabitats I Project) and the European Commission (Ecap‐MED II Project; projects implemented by UNEP/MAP‐RAC/SPA). Alfonso Ramos was supported the CIESM “Tropical Signals,” Stelios Katsanevakis and Maria Sini were supported by the Project “Coastal Environment Observatory and Risk Management in Island Regions AEGIS+” (MIS 5047038), implemented within the Operational Programme “Competitiveness, Entrepreneurship and Innovation” (NSRF 2014–2020), co‐financed by the Hellenic Government (Ministry of Development and Investments) and the European Union (European Regional Development Fund, Cohesion Fund), Stelios Katsanevakis, Maria Sini and Vasilis Gerovasileiou acknowledge the support of the MARISCA Project, co‐funded by 85% by the EEA GRANTS, 2009–2014, and 15% by the Public Investments Programme (PIP) of the Hellenic Republic. Ivan Guala and Daniele Grech thanks thank the support of the project “Pinna nobilis—ricerca per la sopravvivenza: un'iniziativa di Citizen Science per tracciare la mortalità di massa di Pinna nobilis in Sardegna” project and the L/7 grant (CUP 87G17000070002) funded by the Regione Autonoma Sardegna. Jean‐Baptiste Ledoux was funded by an assistant researcher contract framework of the RD Unit—UID/Multi/04423/2019—Interdisciplinary Centre of Marine and Environmental Research—financed by the European Regional Development Fund (ERDF) through COMPETE2020—Operational Program for Competitiveness and Internationalization (POCI) and national funds through FCT/MCTES (PIDDAC). This research was supported by the Strategic Funding UIDB/04423/2020 and UIDP/04423/2020 through national funds provided by the FCT—Foundation for Science and Technology and European Regional Development Fund (ERDF), in the framework of the program PT2020. Bernat Hereu and Cristina Linares acknowledge the support of the long‐term monitoring programme of the catalan Natural Parks, funded by the Departament de Territori i Sostenibilitat of the Generalitat de Catalunya. Cristina Linares acknowledges the support of the ICREA Academia programme. David Díaz acknowledges the support the research grant CTM2016‐77027‐R of the Programa Estatal de Investigación, Desarrollo e Innovación Orientada a los Retos de la Sociedad and Program of Marines Strategies of Spain funded by MITERD. Jamila Ben Soussi was partially funded by the Fondation Albert 2 Monaco (MIMOSA Project) and the Tropical Signals Program of CIESM. Giovanni Chimienti was supported by the Italian Ministry of Education, University and Research (PON 2014–2020, AIM 1807508–1, Linea 1), by the Ente Parco Nazionale del Gargano (Research agreement with CoNISMa N. 21/2018), and by the National Geographic Society (Grant EC‐176R‐18). Nathaniel Bensoussan acknowledges financial support from the European Commission through the programme “Caroline Herschell” in the context of the action “Developing Downstream applications and services on BIO‐PHYsical characterization of the seascape for COASTal management” (BIOPHYCOAST). Monica Montefalcone and Annalissa Azzola collected some of their data on MMEs in the frame of the project “Mare Caldo” funded by Greenpeace Italy. Núria Marbà acknowledges financial support from the Spanish Ministries of Economy and Competitiveness (CTM2012‐32603, CGL2015‐71809‐P) and Ministerio de Ciencia, Innovación y Universidades (RTI2018‐095441‐B‐C21). Diego K. Kersting acknowledges support by the postdoctoral fellowship programme Beatriu de Pinós funded by the Secretary of Universities and Research (Government of Catalonia) and the Horizon 2020 programme of research and innovation of the European Union under the Marie Sklodowska‐Curie grant agreement No 801370.

Appendix S1.

Click here for additional data file.

Appendix S1.

Click here for additional data file.