Recurrent extreme climatic events are driving gorgonian populations to local extinction: low adaptive potential to marine heatwaves

Extreme climatic events (ECEs), such as marine heatwaves (MHWs), are a major threat to biodiversity. Understanding the variability in ecological responses to recurrent ECEs within species and underlying drivers arise as a key issue owing to their implications for conservation and restoration. Yet, our knowledge on such ecological responses is limited since it has been mostly gathered following “single-event approaches” focused on one particular event. These approaches provide snapshots of ecological responses but fall short of capturing heterogeneity patterns that may occur among recurrent ECEs, questioning current predictions regarding biodiversity trends. Here, we adopt a “multi-event” perspective to characterize the effects of recurrent ECEs and the ecological responses in Paramuricea clavata, a Mediterranean temperate coral threatened by MHWs. Through a common-garden experiment repeated three consecutive years with the same individuals from three populations, we assessed the respective roles of environmental (year effect), genetic (population effect) and phenotypic (population-by-environment interactions effect) components in the ecological response to recurrent heat stress. The environmental component (year) was the main driver underlying the responses of P. clavata colonies across experiments. To build on this result, we showed that: i) the ecological responses were not related to population (genetic isolation) and individual (multilocus heterozygosity) genetic make-up, ii) while all the individuals were characterized by a high environmental sensitivity (genotype-by-environment interactions) likely driven by in-situ summer thermal regime. We confront our experimental results to in situ monitoring of the same individuals conducted in 2022 following two MHWs (2018, 2022). This confirms that the targeted populations harbor limited adaptive and plastic capacities to on-going recurrent ECEs and that P. clavata might face unavoidable population collapses in shallow Mediterranean waters. Overall, we suggest that biodiversity forecasts based on “single event” experiments may be overly optimistic and underscore the need to consider the recurrence of ECEs in assessing threats to biodiversity.


Introduction
Extreme climate events (ECE) linked to climate change such as heatwaves pose significant challenges for biodiversity (Jacox et al., 2020;Maxwell et al., 2019;Parmesan et al., 2000;Pinsky et al., 2019;Ummenhofer & Meehl, 2017).ECEs have been associated with the increased frequency of mass mortality events (MMEs), accelerating species demographic decline (Leung et al., 2017;Smale et al., 2013) and questioning the future of biodiversity (Smale et al. 2019).Besides demographic decline, field surveys revealed heterogeneity in the patterns of ecological responses to ECEs across taxonomic (species, populations, and individuals), spatial (regions) and temporal (years) scales (van Bergen et al., 2020;Hughes et al., 2003;Pansch et al., 2018).Yet, to date, the potential effects of recurrent ECEs and the related temporal variability in ecological responses has been poorly considered (but see Ahrens et al., 2021;Brown et al., 2023;Brown & Barott, 2022;Husson et al., 2022;Logan et al., 2014;Regan & Sheldon, 2023).Indeed, most of our knowledge on ecological responses to ECEs has been gathered studying one particular event following a "single-event approach" (Altwegg et al., 2017).This approach provides a snapshot of the ecological responses but falls short of revealing the consequences of recurrent ECEs on biodiversity.Consequently, the need to develop a "multi-events" perspective has been recently acknowledged to improve predictions on species abilities to face ECEs (Bailey & Van De Pol, 2016).
Ecological responses to ECEs, considered here as phenotypes, are shaped by the interplay among 'genetic ', 'environmental' and 'plastic' components (Merilä & Hendry, 2014).This interplay has been investigated using common garden experiments (Malyshev et al., 2016) and long-term time series (Regan & Sheldon, 2023).The 'genetic' component relies on the standing genetic variation (Dixon et al., 2015), shaped by the balance between evolutionary forces (e.g.local adaptation, genetic drift; Bay & Palumbi, 2014).The 'environmental' component includes biotic (e.g.species interactions) and abiotic (e.g.temperature) factors that characterize a specific habitat and influence ecological responses (Scheiner, 1993).Finally, 'plastic' components result from the interaction of genetic and environmental components.
This "population-by-environment" and "genotype-by-environment" interactions underlie plasticity in phenotypic response at the population and individual levels (Matesanz et al., 2019).
Individual phenotypic plasticity, considered as a single genotype expressing different phenotypes function of the environment, can provide a short-term buffer allowing organisms to immediately face ECEs before genetic adaptation occurs (Chevin et al., 2010;Reusch, 2014).
Marine heatwaves (MHWs), known as discrete periods of anomalously warm water (Frölicher et al., 2018;Smith et al., 2023), are some of the most challenging ECEs for marine diversity, threatening ecosystem's structure and functioning (Smale et al., 2019) from kelp forests (Arafeh-Dalmau et al., 2020) to coral reefs (Genin et al., 2020;Gómez-Gras et al., 2021a;Hughes et al., 2021).In the last two decades, the Mediterranean Sea has been recurrently impacted by MHWs driving mass mortality events (MMEs) impacting multiple phyla of benthic macro-invertebrates (Cramer et al., 2018;Garrabou et al., 2022).In this particular case, as in other marine ecosystems, field surveys (Garrabou et al., 2001) and long-term monitoring studies (Gómez-Gras et al., 2021b;Montero-Serra et al., 2019) combined with "single event" experiments in controlled conditions (Crisci et al., 2017;Gomez-Gras et al. 2022), provided insights into the heterogeneity of ecological responses within species impacted by MHWs.Yet, whatever the system considered, the impact of genetic, environmental, and plastic components underlying the differential ecological responses remain poorly understood, particularly, in the context of recurrent MHWs (but see Hughes et al., 2021).
We aim to advance the characterization of the impacts of recurrent ECEs on withinspecies diversity and to infer the respective roles of the genetic, environmental, and plastic components in the ecological responses.We adopt a multi-events perspective combining experiments and in-situ survey focusing on the Mediterranean scene and on the red gorgonian Paramuricea clavata (Risso, 1826).This habitat-forming octocoral is a well-suitable candidate given the impacts of MHWs on shallow populations monitored for 20 years (Garrabou et al., 2021, see below).We repeated during three consecutive years (2015,2016,2017) a commongarden thermal stress experiment at fine spatial scale (populations separated by < 1 km), in which we controlled for different aspects of genetic (same genotypes tested) and environmental (same experimental set-up) components.Specifically, we; i) estimated the percentage of heterogeneity in the ecological responses, respectively explained by the taxonomic (individual and populations) and temporal (years) variabilities; ii) tested the significance of the genetic (population effect), environmental (year effect) and plastic (population-by-years effect) components on the ecological responses.We discussed the obtained results in the light of iii) the populations and individuals' genetic make-up (i.e.genetic drift and heterosis); iv) environmental sensitivity analyses looking at genotype-by-environment interactions and v) insitu summer thermal regimes.We vi) contrasted the ecological responses from the experiments to the ecological responses reported from a field survey of the same individuals conducted in 2022 after two MHWs.
Compared to previously published single event approaches, our results point toward a lack of adaptive potential, whether plastic or genetic, of the shallow populations of P. clavata to the recurrent MHWs.We suggest to cautiously consider predictions of biodiversity trends based on "single event" approaches.
"Single event" experiments using common garden set-ups identified the thermal risk zone of this species for temperatures over 25º C (Crisci et al 2017).These experiments have characterized ecological responses and underlying processes in populations from contrasted environments (depth) and separated by distances of tens to thousands of kms (Arizmendi-Mejía et al., 2015b;Crisci et al., 2017).Altogether, the ecological importance of P. clavata as a habitat-forming species and the extensive available ecological information make this model highly relevant to evaluate the roles of the genetic, environmental, and plastic components on the ecological responses to recurrent ECEs.

Identification of the colonies:
Ninety adult colonies (>50 cm) were randomly chosen around 15 m depth from three different sites separated by hundreds of meters at Medes Islands, Spain (42°02'60.00"N 3°12'60.00"E): Pota del Llop (N=30), La Vaca (N=30) and Tascons (N=30, Figure 1a).Each colony was marked during September of 2014 using plastic tags with a unique ID (Figure 1b).From every marked colony, two apical fragments of 10 cm were collected between the 21 st of September and 22 nd of October of 2015, 2016 and 2017.These fragments were retrieved in 2 L bags of water and immediately transported (maximum transportation time 2h) in coolers to the Aquarium Experimental Zone (ZAE) of the Institut de Ciències del Mar (ICM-CSIC, Barcelona, Spain).

Common garden set-up repeated in 2015, 2016 and 2017:
Upon arrival, colony fragments were mechanically fixed to experimental plates as described in Crisci et al. (2017) (Figure S1).All colonies were acclimated for one week in an open water system with 50 μm sand-filtered running seawater at 17-18°C.No mortality signs and tissue necrosis were detected during this period in any of the colonies and/or years.Moreover, all sampled colonies showed active polyps during feeding.After the acclimation period, we conducted the common garden experiments as described in Crisci et al. (2017).Specifically, each fragment was divided into two fragments, one for control (18°C) and one for the heat stress treatment (25°C).For the heat stress treatment, temperature was increased from 18 to 25 °C over a period of 3 days and maintained at 25°C for the next 28 days.Colonies were fed during the entire experimental course (see Supplementary information for further details about the set up).
Phenotypic response: survival analysis, individual fitness, and modeling the response to thermal stress Several descriptors were then used to statistically compare the ecological responses (differences in thermotolerance) among colonies and populations.The percentage of tissue necrosis (extent of injury) by colony was monitored visually every day.Ecological impacts on P. clavata were considered as "low" when 10-30% necrosis was observed, "moderate" when necrosis was >30-60% and "severe" for >60% necrosis, following the impact classification of the T-MEDNet mass mortality database (e.g.Garrabou et al., 2019).In addition, we estimated for each population and year: i) the daily extent of injury per colony (% of tissue necrosis), ii) the daily percentage of affected colonies (>10% of tissue necrosis), and iii) the daily percentage of dead colonies (100% tissue necrosis).
To model the response to thermal stress across years and populations, we computed a Principal Component Analysis (PCA) with individuals and years as cases and the percentage of necrosis per day as variables.A first exploratory PCA explained 50.29% of the variance and was strongly correlated with the percentage of necrosis from day 10 th to 28 th (Figure S2a).To improve the fitness of the data, days 1 to 9 were removed from the data set.Thus, the second PCA explained 75% of the variance for the remaining days (Figure S2b).All PCAs were performed with PCA() function from FactoMineR R package (R Core Team, 2022; Lê et al., 2008).Scores from the second PCA were considered as a proxy for individual fitness as they were correlated to the baseline percentages of necrosis.These scores were employed as response variable for two linear models: 1) including 'population' and 'year' as predictor variables (i.e.fixed factors); 2) using 'individual' as a random factor, added to the previous predictor variables, and affecting only the model intercepts.Available data were not sufficient to fit all the extra number of parameters (i.e.coefficients), thus the effect of the factor 'individual' was not tested in the slope of the model.We consider the 'individual' effect in the intercept as a different baseline of resistance to thermal stress.Linear models were fitted using lm() and lmer() functions from the stats and the lme4() R-packages (Fox & Weisberg, 2019).
To meet the assumptions of residual normality and homoscedasticity, we transformed the response variable with the Box-Cox transformation, implemented with the boxcox() function from MASS R package (Venables & Ripley, 2002).These assumptions were tested with the Shapiro-Wilks and Levene tests using Shapiro.test() and leveneTest() from car R-package (Fox & Weisberg, 2019).The best model was selected based on AIC (Akaike information criterion) using the function anova() (R Core Team, 2022) and a post-hoc Tukey test with the glht() function from multcomp R-package (Hothorn et al., 2008).We estimated the percentage of the contribution of each factor to the total variance computing the difference in log-likelihood between models with and without each factor (individual, population, and year) pondered by the degrees of freedom.The random intercepts obtained for every individual in the linear model 2 (hereafter "nec-int") were used as a proxy for individual fitness to test for heterosis (see below).

Genetic components: population structure and individual heterosis.
DNA extractions, genotyping protocols, microsatellite characteristics, quality check, and analyses of genetic diversity are described in Supplementary information.We characterized the genetic diversity and structure of the three sampling sites by genotyping 87 individuals collected in at least two of the three years with 14 microsatellites.
We conducted a Discriminant Analysis of Principal Components (DAPC; Jombart et al., 2010) from adegenet R-package (Jombart et al., 2010), by implementing the function find_clusters() and specifying a maximum number of clusters with max.n.clust = 3.A maximum number of 100 PCs was chosen and the lowest value of Bayesian Information Criterion (BIC) was applied to estimate the number of clusters.GENEPOP 4.1.4(Rousset, 2008) was used to compute the overall and pairwise FSTs estimators from Weir & Cockerham (1984).Genotypic differentiation between sites was tested using an exact test (Raymond & Rousset, 1995) with default parameters.In small and isolated populations, inbreeding depression can negatively impact individual fitness constraining the response to ECEs (Fitzpatrick & Reid, 2019).Accordingly, we estimated the genetic differentiation proper to each site by calculating the site-specific FST and 95% High Probability Density Intervals in GESTE (Foll & Gaggiotti, 2006).The sitespecific FST estimates the relative impact of genetic drift on the differentiation of the considered site relative to the remaining ones.
At the individual level, a positive correlation between heterozygosity and fitness-related traits, known as heterosis, has been reported in some species (David, 1997).We looked for heterosis in the response to thermal stress testing the correlation between the values of nec-int as proxy for individual fitness and the standardized individual multilocus heterozygosities (sMLH) computed using the R package InbreedR() (Stoffel et al., 2016).The slope of the linear model [lm(nec-int ~ sMLH)] was compared to its null distribution obtained with a Monte-Carlo permutation test with 10,000 permutations.

Environmental component: summer thermal environments
We analyzed and compared the in situ thermal regime at 15m depth experienced by the colonies before sampling.Temperature data for Medes islands was obtained from the T-MEDNet database, which follows a continuous temperature series since 2004(von Schuckmann et al., 2019)).We assessed the recent summer local thermal regime calculated over the 3-months period of June, July and August, prior to collection in September.We consider the timing and magnitude of the daily temperature cycles and the exposure to warm conditions; thus we considered the averages of maximum temperatures during the warmest periods of the year, number of days with high temperatures and the ecological threshold T23 (i.e.number of extreme heat days at ≥23ºC).In addition, we considered positive temperature anomalies as the number of marine heat spikes (MHS) above the Inter-annual percentile 90 th (iT90 threshold) lasting less than five days, while prolonged discrete periods of anomalously warm water surpassing the iT90 percentile for more than 5 days, were considered as MHWs.Presence of MHWs were detected with Python module marineHeatWaves (Hobday et al., 2016), while impact categories registered at Medes islands were assigned following the classification of the T-MEDNet database (Garrabou et al., 2019;Hobday et al. 2018).MHWs categories were set as significant temperature foldchanges from the Inter-annual mean temperature and the iT90 threshold.Thus, "Low" classifications correspond to the temperature ranges between the interannual mean temperatures and iT90, while "Moderate" are considered as a 1-2-fold, "Strong" as a 2-3-fold and "Severe" as a 3-4-fold, respectively.

Plastic component: 'genotype-by-environment' interactions and environmental sensitivity
We estimated the variability in the 'genotype-by-environment' interactions by characterizing the environmental sensitivity of each genotype following Falconer & Mackay (1996).We computed three environmental values corresponding to the yearly mean phenotypes (i.e.mean PCA scores over individuals for each year), considering the 66 genotypes present during the three years.We then plotted each individual phenotype (PCA score) against the environmental value for each year and we computed the regression slope, which is considered as an estimator of the environmental sensitivity of the genotype (Falconer & Mackay, 1996).We used this plot to classify the sensitivity of the genotypes in three categories adapted from Bonacolta et al. (2024).Resistant genotypes were expected to show low intercepts in the first year and approximated null slopes (i.e.low and constant level of necrosis in the 3 experiments).Hypersensitive genotypes were expected to show high intercepts in the first year and null slopes (i.e.high and constant level of necrosis in the 3 experiments).Finally, sensitive genotypes were expected to show low intercepts in the first year and positive slopes (i.e.increasing level of necrosis through time).

Field survey of necrosis rates following 2018-2022 MMEs:
Medes Islands were impacted by two MHWs in 2018 and 2022 with associated mass mortality events (Garrabou et al., 2022;Zentner et al., 2023).We surveyed by scuba diving the percentage of tissue necrosis in situ for the same individuals used in the experiments.This survey was done in October 2022.Differences in mean tissue necrosis of colonies between populations were tested using a parametric one-way ANOVA, followed by post hoc Tukey's HSD tests.

Phenotypic responses of P. clavata during the three common garden experiments
Signs of tissue necrosis were observed for all populations in the three years.In 2015 and 2016, the mean extent of injury was of moderate impact with values below the 60% at the end of the experiment (day 28 th ; mean ± SE): 32.14% ± 7.25 and 38.33% ± 6.67 for La Vaca; 27.03% ± 6.20 and 35.16% ± 7.09 for Pota del Llop and 37.89% ± 7.26 and 41.2% ± 6.91 for Tascons (Figure 2).On the contrary, in 2017, colonies showed severe impacts with >60% of average tissue necrosis earlier by day 14 th in all populations, and all colonies died by day 18 th (100% of tissue necrosis) with the exception of Pota del Llop which reached 100% of tissue necrosis by day 24 th (Figure 2).

Individual fitness: modelling the response to thermal stress
The linear model 2 including the random factor 'individual' was retained (lowest statistically significant AIC = 710.88;df = 7 [Chi-square = 15.5, p < 0.001]).The 'individual' random factor had a significant effect, suggesting that each individual has a different baseline of resistance to necrosis.Regarding the fixed factors of the model, the deviance test showed that only the 'year' factor was significant (Table S2a), while post-hoc Tuckey test showed that significant differences were due to the year 2017 (Table S2b).The factor 'year' was contributing to 95.01% of the variance of the data, followed by the random factor 'individual' with a contribution of around 11 times less (4.1%).Finally, 'population' and 'population-byyear interaction' were non-significant and showed the lowest contributions: 0.5% and 0.39%, respectively.

Environmental component: thermal environments
Recent thermal history patterns, considered as June, July and August, revealed similar mean ± SD temperatures for the three years: 21.3ºC ± 1.6 in 2015, 21.6ºC ± 1.0 for 2016 and 21.8ºC ± 1.3 for 2017 (Table S7).Concomitantly, extreme heat days (T23) were detected in all years during the summer season (Table S7).For 2015 and 2016, at 15m depth, maximum summer temperatures reached 24.8 and 23.5ºC, respectively, and a low number of total extreme heat days exposure to ≥23ºC was recorded (N=12 for 2015 and N=2 for 2016; Table S7, Figure S6ab).The years 2015 and 2016 revealed several periods of anomalous high temperatures during the summer season, but no MHWs (Figure 4, Table S8).Interestingly, the year 2017 reported the highest maximum temperatures of 24.9 ºC with a total of 19 days of exposure at extreme temperatures (≥23 ºC), surpassing the thermal limit of P. clavata (Table S7, Figure S6).In addition, unlike for 2015 and 2016 where no MHW occurred, two MHWs occurred in 2017 with strong, and severe classifications from June to July and several heat spikes (Figure 4).The mean maximum temperature for these MHW was 24.1ºC (Table S8).and 90 th percentile (iT90, red dotted line).Days below the red solid line were considered as "cool days", whereas days above were considered as "warm days".Following Hobday et al.

Plastic component: 'genotype-by-environment' interactions and environmental sensitivity
Regarding the sensitivity analyses, environmental values (yearly mean PCA scores) were -0.08 ± 1.46 for 2015, 0.34 ± 1.55 for 2016, and 3.56 ± 0.55 for 2017.All individuals but one show positive slopes ranging between 0.16 and 1.65.Environmental sensitivity varied by a factor of 10 among individuals (Figure 5).Following our classification, genotypes were shared among hyper-sensitive and sensitive categories with no resistant genotypes (Figure 5).

Discussion
We combined replicated experiments, population genetics, and an in-situ field survey to reveal the prominent influence of the environmental component (likely yearly variations in summer thermal regime) on the heterogeneity in ecological responses to thermal stress in P. clavata.
The low influence of genetic and plastic components combined to the high environmental sensitivity of the tested genotypes point toward a dramatically low adaptive potential to recurrent MHWs.This "multi-event" perspective strengthens the recent call to carefully consider predictions on biodiversity evolution based on single-event experiments.
The heterogeneity in the ecological response to thermal stress is mostly driven by the thermal regime during summer.
The three populations of Paramuricea clavata showed high levels of tissue necrosis (moderate to severe mortality) at the end of each experiment confirming the species sensitivity to thermal stress (Crisci et al., 2017;Gómez-Gras et al., 2022).The ecological responses were similar between populations in 2015 and 2016 (moderate mortality) compared to 2017, in which colonies died before the end of the experiment (severe mortality).More than 95% of the variance in the ecological response was explained by the factor year (environmental component), while factors individual, population (genetic component), and the population-byyear interaction (plastic component) explained only 4.1%, 0.5%, and 0.39%, respectively.
Factors individual and year (environmental component) were significant, which suggests different baselines of resistance to thermal stress among individuals and confirms the major environmental effect, mostly driven by the year 2017.This environmental effect was refined by the sensitivity and thermal regime analyses.First, the environmental value (mean phenotype) for 2017 was two to three times higher than the values for 2015 and 2016.Then, we reported positive regression slopes for almost all genotypes supporting an increased negative impact of environmental conditions (environmental sensitivity) from 2015 to 2017.Summer conditions during 2017 showed the largest number of MHS and MHWs compared to 2015 and 2016.Consequently, we posit that colonies of P. clavata were driven close to their physiological limits by the 2017 extreme summer conditions, which may have hampered any adjustment to thermal stress, whether genetic or plastic, during the experiment.
The relative impact of environmental, genetic, and plastic components in differential responses to ECEs has been screened in different species.For example, ubiquitous population-byenvironment interactions (plastic component) have been detected in 172 species of plants (Matesanz & Ramírez-Valiente, 2019), but are lacking in others (Shao et al., 2022).Single event experiments with tropical corals identified local adaptation (e.g.Thomas et al., 2022) and adaptive plasticity (e.g.Drury et al., 2022) as drivers of differential bleaching responses.Here, we found relatively similar and high levels of necrosis among populations with a prevailing impact of the environmental component.These findings contrast to our previous studies based on "single event" experiments at larger spatial scales where differential ecological responses (distinct necrosis levels) were observed among populations in some cases (Arizmendi-Mejía, et al., 2015a;Arizmendi-Mejía et al., 2015b;Bonacolta et al., 2024;Crisci et al., 2017;Gómez-Gras et al., 2022).Two non-exclusive hypotheses relying on the spatial and temporal features of the experiments might explain this apparent discrepancy.The experiments conducted to date in P. clavata have considered a wide range of geographic distances, from local to regional (Arizmendi-Mejía et al., 2015b;Crisci et al., 2017) and interregional (Bonacolta et al., 2024;Gómez-Gras et al., 2022) scales.These experiments demonstrated population heterogeneity in ecological responses to thermal stress triggered by different drivers (e.g.genetic isolation, microorganisms).Considering that the impact of the genetic component on ecological responses can vary over spatial distances (e.g.Galloway & Fenster, 2000;Joshi et al., 2001), we hypothesize that the lack of significant population effect observed here can be related to the fine spatial scale of the experiment, which flattened the differences between populations.Contrary to previous experiments, we targeted three populations from similar habitats at the same depth range in a close spatial proximity (hundreds of meters).The potential for contrasted genetic make-up at such fine spatial scale is low (but see Ledoux et al., 2015;Richardson et al., 2014) as supported by the comparable levels of genetic isolation of the three populations (overlapping population-specific FSTs).Re-analyzing the different experiments in P. Clavata accounting for ecological or spatial distance between populations should allow to go further in this hypothesis.
Considering a temporal perspective, the discrepancy among experiments in P. clavata can result from an intensification of extreme climatic events from summer 2009 (Crisci et al., 2017) to summer 2019 (Gómez-Gras et al., 2022) and summer 2022 (this study; Rovira et al. 2024), which could have driven the colonies closer to their physiological limits in the later and warmer years.That is, the summer thermal regime previous to the first experiment (summer 2009(summer , Crisci et al., 2017) )  MHWs.Third, one of the strongest mortality events ever observed was reported in 2022 (Estaque et al., 2023, Rovira et al. 2024) corroborating the rise in MHWs in the Mediterranean.
Strengthening this hypothesis, the intensification of disastrous ECEs in the last decades is not peculiar to the Mediterranean (see Stillman, 2019).The frequency of bleaching events in tropical corals increased worldwide since 1980 (Hughes et al., 2018), with detrimental cumulative effects of heatwaves in the last ten years (Hughes et al., 2021).

What's next for shallow populations of P. clavata?
Our results question the persistence of the studied shallow populations of P. clavata.Both the experiments and field surveys conducted here supported a limited potential adaptability whether based on genetic or plastic component.Yet, we revealed some variability in the baseline ecological responses at the individual levels (4.1%).While this individual variability may have been seen as a hopeful "raw material" for adaptation to MHWs, the sensitivity analyses and the field survey showed how it was almost totally squeezed in 2017 (experiment) and 2022 (in the field).Dedicated studies combining population genomics, environmental, and mortality data are needed to further look for genome-environment associations, potential outlier loci involved in the differential responses and to estimate of the genomic offset of P. clavata.However, the lack of clear adaptive potential revealed here combined to the species life history traits (e.g. generation time >12 years, Coma et al., 1995) are such in a contrast with the current MHWs temporal dynamics that any evolutionary response seems compromised.
The potential for adaptation is also questioned in tropical corals in which candidate genetic loci identified to date only show relatively elusive influence on heat stress tolerance (Fuller et al., 2020, but see Matz et al., 2020).In the same collapsing line, the absence of population-byenvironment interactions suggest limited potential for evolutionary changes in adaptive plasticity (Sirovy et al., 2021).Recent studies point towards ecological memory, an increase in stress tolerance following previous exposure, as a key mechanism for coral acclimation to MHWs (Hackerott et al., 2021;Hughes et al., 2018).Yet, results are contrasted among species with a decrease in bleaching sensitivity following repeated heat stress in some species but not in others (Brown et al., 2023;Hughes et al., 2021).Our study allows first insights into P. clavata environmental memory.First, colony fragments used in a particular year were submitted to summer thermal conditions of the previous years.Yet, the worst ecological responses to thermal stress were observed during the last experiment in 2017 with high environmental sensitivity (positive slopes and lack of resistant genotypes in the sensitivity analyses).Then, most marked colonies showed necrosis during the field surveys following 2018 and 2022 MHWs events.These results question any increase in thermotolerance as expected with the ecological memory hypotheses and strengthen the limited potential for adaptation in P. Clavata.

Conclusion
As temperature and frequency of ECEs continue to rise (Hughes et al., 2021;Garrabou et al., 2022), the evolution of biodiversity is more than ever a central concern for society.Adopting a "multi-event" perspective that combined replicated common-garden experiments in aquaria and mortality surveys in the field performed on the same colonies, our study points toward an inevitable collapse of many of the shallow populations of P. clavata.This collapse would emerge from a low to non-existent adaptive response, whether driven by genetic or plasticity, combined to a high environmental sensitivity and a potential intensification of MHWs.
Considering the small spatial scale of our study, extrapolation at larger scale should be made cautiously.Yet, population collapses of P. clavata linked to recurrent MHWs have been observed in other Mediterranean regions (Garrabou et al., 2021;Gómez-Gras et al., 2021b).
Moreover, field surveys following the 2022 MHW event in this and in other regions reported terrific mortality rates.Other populations from the same region showed a similar proportion of total affected colonies (70%) with almost 40% of the tissue necrosis (Rovira et al. 2024).
Hundreds of km apart, populations until 20m depth displayed on average >80% of affected colonies and an increase by 142% of the degree of impact following the 2022 MHW compared to the previous MME in 2003 (Estaque et al., 2023).Worrying, this trend in P. clavata could likely be transposable to many of the Mediterranean habitat-forming and sessile species impacted by MHWs (Garrabou et al., 2022;Gómez-Gras et al., 2021a).We predict a shift in these species' upper distribution limits, which will lead to a simplification of associated benthic communities hampering potentially related ecosystem functions and services (Gómez-Gras et al., 2021a).
This study echoes two recent calls regarding the future of marine diversity in the context of extreme climatic events.First, the impacts of ECEs on biodiversity should be studied from a temporal perspective, which accounts for their recurrence (Hughes et al., 2021).In this line, we suggest that predictions of biodiversity evolution based on "single event" approaches should be considered cautiously as they can be overly optimistic.Then, while conservation and restoration actions should be able to slow, and/or to some extent reverse, locally the collapsing trend followed by many marine habitat-forming species, immediate action on greenhouse gas emissions remains the only way to protect these species globally.

Figure 2 .
Figure 2. Average tissue necrosis (mean extent of injury ± SE) of P. clavata colonies during the 28 days of exposure in common garden experiments for La Vaca, Pota del Llop, and Tascons in 2015 (green), 2016 (blue) and 2017 (yellow).Mortality severity is highlighted in light yellow (low), orange (moderate), and red (severe).

Figure 3
Figure 3. a) Individual membership probabilities are represented by a vertical line, where the different color segments indicate the individual proportion to each cluster (K=3) as estimated by the discriminant analysis of principal components (DAPC): Pota del Llop (light-blue), Tascons (orange) and La Vaca (brown).Mean membership probabilities are given above each colored cluster b) Scatter plot of the DAPC: each dot corresponds to one individual (N=87) from each of the three populations.

Figure 4 .
Figure 4. Daily mean temperature values recorded at 15m from June to September (Summer season) at Medes Islands with respect of the inter-annual climatological mean (red solid line)

Figure 5 .
Figure 5. Sensitivity analyses and genotype-by-environment interactions of P. clavata during 2015, 2016 and 2017 experiments.Individual ecological response (PCA scores) was plotted against each environment value (mean PCA score by year).The slope of the regression for each individual is considered as an estimator of the environmental sensitivity of the individual.Populations are displayed in colors: Pota del Llop (light-blue), Tascons (brown), and La Vaca (purple).
Intraspecific differences in the response to thermal stress: does the geographic scale matter or are recent summer conditions overwhelming P. clavata physiological capacities?
was less stressful than the 2017 summer thermal regime (this study) allowing some colonies to face the 2009 experimental stress while colonies were totally swept by the 2017 experimental stress.Three main points support this hypothesis.First, the environmental sensitivity analysis considering 2015, 2016, and 2017 experiments clearly show a decrease in the variation of phenotypic responses among the individuals and an increase of the yearly environmental value between 2015/2016 and 2017.Individuals that showed relatively low necrosis in the first two years were as strongly impacted in 2017 (positive regression slopes) as individuals showing high necrosis during the first two experiments (regression slopes ~ 0).Second, necrosis was observed in a vast majority of the marked colonies in the three populations (>70%) during the field survey following the 2018 and 2022