Thermal Performance Reects Evolutionary Legacies of Seaweeds and Seagrasses Across a Regional Climate Gradient

Comparative patterns in thermal performance between populations have fundamental implications for a species thermal sensitivity to warming and extreme events. Despite this, within-species variation in thermal performance is seldom measured. Here we compare within-species variation in thermal performance across the Mediterranean Sea, with between-species variation within communities, for two species of seagrass (Posidonia oceanica and Cymodocea nodosa) and two species of seaweed (Padina pavonica and Cystoseira compressa). Experimental populations from four locations representing approximately 75% of each species global distribution and a 6ºC gradient in summer temperatures were exposed to 10 temperature treatments (15ºC to 36ºC), reecting median, maximum and future temperatures. Thermal performance displayed the greatest variability between species, with optimal temperatures differing by over 10ºC within the same location. Within-species differences in thermal performance were also important for P. oceanica which displayed large thermal safety margins within cool and warm-edge populations and small safety margins within central populations. Overall, experimental upper thermal limits reected genus-level realised thermal limits, more than realised species-limits or maximum local temperatures. Our ndings suggest patterns of thermal performance in Mediterranean seagrasses and seaweeds retain deep ‘pre-Mediterranean’ evolutionary legacies, resulting in unexpected patterns of vulnerability to warming within benthic marine communities.


Introduction
Ocean warming is having profound effects on benthic marine ecosystems across the globe 1 . Foundation species such as seaweeds, seagrasses and corals can be particularly susceptible to warm temperatures and have undergone extensive thermal stress 2 , mortality 3 , and range contraction 4 over the past decade alone. While a lot of these impacts have been concentrated in warm edge populations 4 , cooler, central populations can also be highly sensitive to extirpation 3,5 , highlighting that not all species or locations are equally susceptible to warming. The difference in temperature between an organisms upper thermal limits and its upper environmental temperatures (herein, thermal sensitivity), can vary considerably both within a species and between species due to differences in thermal physiology and through geographic differences in thermal regimes 6 . This variation presents a taxonomically and geographically diverse range of responses to warming and raises fundamental questions about the ecological and evolutionary drivers of thermal performance patterns.
From a physiological perspective, expectations for how thermal sensitivity differ among populations across a species geographical range, depends on the variation in thermal performance between individuals and populations. At one extreme, if individuals display similar thermal limits irrespective of range position, then thermal safety margins will decrease toward the warm edge of a species range, creating a gradient of increasing thermal sensitivity from cool to warm edge populations 6 . This is a common assumption used in many species distribution models, for example, that rely on the realised distribution of species to predict the likelihood of extirpation or range shifts under climate change 7, but see 8,9 . If, on the other hand, thermal limits differ between populations -through local adaptation or acclimatisation to local thermal regimes -then thermal sensitivity can remain constant across a species range 5 . This, is an inherent assumption used to estimate the thermal sensitivity of reef building corals, for which the magnitude and the duration of local thermal anomalies are strong predictors of thermal bleaching 10 . In reality, most species may not have perfectly conserved, or perfectly adapted thermal niches, but lie somewhere within this spectrum. For the vast majority of species, within-species patterns of thermal sensitivity remain unknown, forcing climate change predictions to rely on measurements of thermal performance from a single location 11 or the realised thermal distribution of a species 12 .
Understanding the variability in thermal performance between populations and verifying the extent to which thermal performance re ects local environmental conditions, realised species distributions, or other factors, is crucial in order to accurately predict climate change impacts and adaptively manage ecosystems in the Anthropocene.
The Mediterranean Sea presents a fascinating region to study within-species and between-species patterns in thermal performance and sensitivity. Seaweeds and seagrasses are the dominant habitat forming species throughout the region, where they harbour high biodiversity and provide essential ecosystem services 13 . Many of the species that characterise rocky reefs (e.g. Padina pavonica) and seagrass meadows (e.g. Posidonia oceanica) in the Mediterranean Sea today, have deep evolutionary histories in the region dating back over 65 million years 14,15 . What's more, these species have survived a remarkably dynamic climatic history, spanning from tropical origins in the Tethys Sea during the Cretaceous, drying out of the basin during the Messinian salinity crisis 6Mya; 16 , to repeated glaciations during the Pleistocene (2.6-0.012 Mya). These upheavals have shaped the regions unique biodiversity to include species with tropical and boreal phylogeographic origins 13,17 . The Mediterranean Sea settled to its current climatic conditions approximately 2000 years ago 13 , characterised by a steep 6ºC climatic gradient from west to east. Despite their evolutionary resilience, seaweeds and seagrasses are highly susceptible to contemporary warming in the Mediterranean. Marine heatwaves, for example, have caused severe impacts to central populations of P. oceanica, leading to predictions that the species will face functional extinction by 2100 18,19 . The high sensitivity in central populations raises questions about the sensitivity of populations at the warm edge of the species distribution in the eastern Mediterranean.
Moreover, the survival of Mediterranean species throughout large climatic uctuations over evolutionary time scales, juxtaposed with the high sensitivity that has been observed to contemporary climate warming, raises fundamental questions about the in uence of evolutionary versus contemporary climatic conditions in shaping the thermal performance of species.
In this study we compare the thermal performance of four pan-Mediterranean species of seagrasses and seaweed, from four populations of each species spanning eight degrees in latitude and 30 degrees in longitude across the Mediterranean Sea. We then examine how thermal performance relates to local environmental temperatures at each location and species-and genus-level realised thermal distributions -as indicators of their evolved thermal a nity. Speci cally, we ask: 1) How does the thermal performance within a species differ across its distribution and 6ºC gradient in ocean temperatures? 2) Do species with similar realised thermal distributions display similar thermal optima and upper thermal limits? 3) How does the variability in thermal performance between different species from the same location compare to the variability between geographically separate populations of the same species?
We discuss these patterns in the context of the evolutionary origins of these species and our capacity to predict their contemporary vulnerability to climate change.
Thermal performance: growth Optimal and critical upper temperature limits of Mediterranean macrophytes varied both between populations and species (Fig. 2). Overall, 11/15 experimental populations displayed signi cant normal growth distributions in response to temperature. Of these, 10/11 signi cant distributions were recorded in experiments from Catalonia, Mallorca and Crete, whereas experimental populations from Cyprus displayed variable growth responses to temperature (Fig. 2, Table S5). For P. oceanica, optimal growth temperatures were recorded at 19 ± 1 ºC for Catalonia and 20 ± 0.7 ºC for Mallorca and Crete, while upper limits were recorded at 30 ºC for Catalonia, Mallorca and Crete. Growth patterns in Cyprus samples displayed a non-signi cant relationship with temperature when all data was included. By removing 28ºC and 32ºC outliers, the growth relationship with temperature improved (p = 0.055) and was characterised by a 27 ± 13ºC optimal temperature (Fig. 2, Table S5). P. oceanica from Cyprus also recorded the highest growth rates of any location during the experiments, averaging 0.2-0.4 cm day − 1 including in the warmest treatments 30-36ºC.
C. nodosa recorded higher optimal and lethal temperature limits than P. oceanica in all locations, except Cyprus. Optimal growth temperatures were 29 ± 0.5ºC, 28 ± 0.4ºC and 34 ± 5ºC in Catalonia, Mallorca and Crete, respectively ( Fig. 2E-G). Critical upper limits were estimated at 41ºC, 37.5 ºC and 48ºC for Catalonia, Mallorca and Crete, respectively, based on best-t growth models. C. nodosa in Cyprus displayed a non-signi cant growth relationship with temperature, despite recording relatively high growth rates across all temperature treatments. C. compressa and P. pavonica displayed wider growth curves than P. oceanica and C. nodosa. Optimal growth temperatures were 21 ± 1ºC and 19 ± 1.5ºC for C. compressa and 15.5 ± 0.5ºC and 18 ± 0.4ºC for P. pavonica, in Catalonia and Mallorca, respectively. Upper thermal limits for C. compressa and P. pavonica were recorded at 32ºC in Catalonia, and 30ºC for C. compressa and 32ºC for P. pavonica in Mallorca. Growth curves for C. compressa and P. pavonica in Crete and Cyprus were non-signi cant (Table S5).

Metabolic rates
At the population level, patterns in gross primary productivity (GPP) and respiration (R) were relatively at or variable for 8/15 and 11/15 populations, respectively, meaning that optimal or critical limits could not be de ned. Among the populations where optimal GPP could be de ned, optimal temperatures ranged between 28 ± 4ºC for Mallorca and Cyprus populations of C. nodosa and C. compressa and 32 ± 8ºC for Crete and Cyprus populations of P. pavonica (Fig. S3, Table S6). Metabolic rates were primarily heterotrophic (i.e. GPP < R) across temperatures in the upper portion of each populations thermal range.
Patterns in thermal safety margins across species range.
Thermal safety margins generally declined toward the warm range edge of species distributions, re ecting conserved thermal limits between populations (Fig. 3). P. oceanica in Cyprus was a notable exception to this pattern, whereby thermal safety margins were higher in Cyprus as a result of higher thermal tolerance at the warm range edge. Thermal safety margins of C. nodosa declined from cooler to warmer populations, albeit with larger thermal safety margins than expected, based on the realised thermal distribution of the species (Fig. 3). C. nodosa in Crete also displayed higher thermal safety margins relative to range position, than the other populations. P. pavonica displayed contrasting patterns in thermal safety margins between growth metrics. Net growth in 2/4 populations where upper thermal limits could be identi ed, displayed declining thermal safety margins toward warmer locations, consistent with expectations of the species realised distribution. In contrast, thermal safety margins based on GPP increased from cooler to warmer populations.
Similarity between thermal performance and local, specieslevel and genus-level realised temperatures.
Experimentally determined optimal temperatures of seagrass and seaweed populations most closely resembled median local habitat temperatures of the experimental population and were typically higher than median temperatures observed across a species' or genus' realised distribution (Fig. 4). Median local environmental temperatures were − 2.1 ± 1.9 ºC cooler than the thermal optima for P. oceanica, when averaged among locations, whereas thermal midpoints of species and genus-level thermal distributions were, on average, -4.6 ± 2.3ºC and − 6.2 ± 2.4 ºC cooler than thermal optimum of P. oceanica, respectively (Fig. 4A). In contrast, upper thermal limits resembled genus-level maximum realised temperatures more than maximum temperatures in the local habitat (Fig. 4B). Upper thermal limits for P. oceanica for example, were, on average, 4.6 ± 1.14ºC above local thermal maxima, compared to only 0.98 ± 0.16ºC above maximum realised temperatures for the species and genus. All other species displayed similar patterns, albeit with larger differences between fundamental and realised limits (Fig. 4B).

Discussion
Quantifying spatial and taxonomic variability in thermal sensitivity is critically important in order to identify vulnerable areas and triage management and conservation efforts in response to climate change.
Here we quanti ed the thermal performance of four dominant habitat forming species of seaweed and seagrass, from four locations across the Mediterranean Sea. We found major differences in thermal performance between species as well as subtle, albeit important differences between populations within species. Interestingly, optimal temperatures most closely resembled local habitat temperatures, whereas upper thermal limits re ected genus-level realised thermal limits. This suggests that thermal performance may have evolved to be optimal under typical local conditions, while upper thermal limits may have been conserved to tolerate rare extreme events. An important exception to this general pattern, was the high thermal tolerance of seagrasses in Cyprus. Cyprus represents the warm distribution edge of many species in the Mediterranean Sea including P. oceanica and yet, displayed the lowest thermal sensitivity due to higher thermal tolerance limits than observed in other locations across the Mediterranean. Observed differences in thermal performance between species and within species provide important context about the potential vulnerability of benthic marine ecosystems to warming and contributes to a general understanding about the underlying drivers of thermal sensitivity.
Major differences in thermal performance between species was the standout nding from this study. This nding is signi cant, because of the close similarities these species share in their realised thermal distributions. Realised thermal limits are a commonly used metric of thermal sensitivity 7  The differences in thermal tolerance between species is also remarkable given the shared selection pressures that have acted on these species. Following the nal closure of the Suez Isthmus, the Mediterranean Sea underwent repeated climatic disturbances that fundamentally affected biodiversity in the basin 17,24 . At the beginning of the Pleistocene, approximately 2.6 Mya, the Mediterranean Sea underwent rapid cooling, leading to mass extinction of much of its tropical-a liated biota 13 . Glacial and interglacial periods then resulted in alternating waves of cool-a liated and tropical species colonizing the Mediterranean from the Atlantic, only to become extinct or have their distribution restricted during the subsequent cool or warm cycle. Evidence of this can be seen in the contemporary distribution of many species (e.g. Sparisoma cretense) with disjunct Atlantic-Levantine distributions due range contractions from the western Mediterranean Basin during glaciations. Despite these profound climatic disturbances, P. oceanica, C. nodosa, C. compressa and P. pavonica remain dominant throughout the Mediterranean, suggesting a high degree of plasticity in all four species. Moreover, differences in thermal tolerance between species persist in spite of strong thermal selection pressures over several million years.
While species-level differences in thermal performance were most conspicuous, more subtle differences in thermal performance were also observable between some populations. Most notably was the thermal performance of P. oceanica in Cyprus, which did not display a decline in growth at higher temperatures as observed in other populations. This nding contrasts with the prevailing paradigm that thermal sensitivity will be highest in warm edge populations. Warm edge populations, by de nition face the highest habitat temperatures of any population across a species range, which if thermal limits are conserved between populations, mean they also have the highest sensitivity to warming. In the case of P. oceanica we observed a u-shaped pattern in thermal safety margins, with central populations living closest to their upper thermal limit and cool and warm range edge displaying the highest thermal safety margins.
This nding provides context for the paradigm that P. oceanica is highly sensitive to warming 3,18 . This idea emerged primarily through studies conducted in the western Mediterranean (e.g. Mallorca) where temperatures above 28ºC have driven severe declines in shoot density of P. oceanica 3 . Our ndings support this previous work insofar as thermal limits of P. oceanica were estimated to be 30ºC, very close to the maximum temperatures that previously led to population declines. However, our ndings also help to explain how P. oceanica can thrive at its warm distribution limit in Cyprus where average daily summer temperatures are already 29ºC and peak daily summer temperatures frequently exceed 30ºC.
With respect to metabolic rates, warming resulted in a reduction of net production in all species except for C. nodosa. A decrease of NP with warming indicates that climate change could result in a change to heterotrophic metabolism, becoming a CO 2 source and an oxygen sink, further aggravating global warming. A net reduction of NP owing to a steeper increase in respiration rates than in gross primary production has been predicted by the Metabolic Theory of Ecology 25 owing to the fact that activation energies for autotrophs are half that of heterotrophs 26,27 . Previous studies have shown a higher increase of respiration than primary production for planktonic communities in response to warming [26][27][28][29][30][31] , whereas for benthic-dominated communities this prediction has been questioned 32 . Our results demonstrate that macrophyte species can also increase their respiration rates faster than primary production with warming, resulting in heterotrophic metabolic rates at high temperatures.
Finally, experimental artefacts need to be considered in regard to their potential contribution to the variable thermal performance results. Thermal performance under experimental conditions is contingent on the experimental design, acclimation rates and seasonality, among other factors 33 . As such, thermal performance measured under experimental conditions needs to be considered as a relative, not absolute estimate of the ecological reality that might occur in a natural ecosystem. Nevertheless, by using standardised methodologies of collection, transportation, acclimation and experimentation in the same experimental system, with identical conditions between experiments, methodological variation was minimised between experiments.
Seasonal timing was one potential source of variation between experiments. While all experiments took place around summer, timing of experiments ranged between June -September, due to logistical constraints. For seaweeds, poor condition of C. compressa from Crete at the time of collection (July) was likely due to seasonal senescence and resulted in it being removed from the analysis. For P. oceanica, it is unlikely that seasonality was responsible for the poor growth relationship to temperature in Cyprus. Growth rates of P. oceanica from Cyprus were the highest recorded in any experiments, suggesting that the poor growth relationship was not an artefact of poor performance. Performance of seagrasses may, however, have been in uenced by experimental conditions as a result of their clonal growth physiology.
Seagrasses are clonal plants connected by rhizomes and are dependent on the connectivity between ramets for nutrient translocation 34 . The consistently negative NP observed across temperatures for seagrasses may be indicative that these species are less suited to growing in fragments under experimental conditions than seaweeds which source their nutrients from the water column.
Our ndings reveal the thermal performance of four dominant habitat forming species across their geographical range in the Mediterranean Sea. Between species differences in optimal and upper thermal limits were pronounced despite all species sharing similar geographical distributions. Within species thermal performance was relatively conserved between populations across a 6ºC climate gradient in the Mediterranean, with the clear exception of P. oceanica at its warm range edge. Our ndings highlight nuanced inter and intra-speci c patterns in thermal performance that would be overlooked through a reliance on realised distributions or measurements from a single population. Contemporary patterns in thermal performance for Mediterranean macrophytes potentially re ect deep evolutionary legacies of species and local differences in the historical and contemporary marine climates. Population-speci c patterns in thermal performance have important implications for the distribution and abundance and conservation of Mediterranean benthic habitats in response to climate change.

Field collections
Thermal tolerance experiments were conducted on two seagrass species (P. oceanica and Cymodocea nodosa) and two seaweed species (Cystoseira compressa and P. pavonica) from four locations spanning 8 degrees in latitude and 30 degrees in longitude across the Mediterranean (Fig. 1

Species description and distribution
The species used in this study are all common species throughout the Mediterranean Sea, although differ in their biological traits, evolutionary histories and thermo-geographic a nities (Fig. S1). Posidonia oceanica is endemic to the Mediterranean Sea with the all other congeneric species found in temperate Australia 14 . The distribution of P. oceanica is restricted to the Mediterranean, spanning from Gibraltar in the west to Cyprus in the east and north into the Aegean and Adriatic seas 35 (Fig. S1A). Genetic diversity is highest around the centre of its distribution, owing to a secondary contact zone around the strait of Sicily, following vicariance during the last glacial maxima 36 . P. oceanica is a slow growing clonal species with peak growth in spring (April-June) and irregular owering events following peak summer temperatures (Aug-October) 37 .
Cymodocea nodosa distribution extends across the Mediterranean Sea and eastern Atlantic Ocean, where it is found from south west Portugal, down the African coast to Mauritania and west to Macaronesia 38 (Fig. S1B). Congeneric species of C. nodosa are found in tropical waters of the Red Sea and Indo-Paci c, suggesting origins in the region at least prior to the closure of the Suez Isthmus, approximately 10Mya. Genetic diversity of C. nodosa is highest in the warm margins of its distribution and lowest in cooler regions 38 . C. nodosa is a fast-growing species with below-ground seed release. Seeds have negative buoyancy and poor dispersal 39 .
Like C. nodosa, Cystoseira compressa has a distribution that extends across the Mediterranean and into the eastern Atlantic, where it is found west to Macaronesia and south to northwest Africa (Fig. S1C). The genus Cystoseira has recently been reclassi ed to include just four species; C. compressa, C. humilis, C.
foeniculacea and C. aurantia 40 . Congeneric Cystoseira spp. all have warm-temperate distributions from the Mediterranean to the eastern Atlantic. C. compressa inhabits shallow rocky shores with peak growth and reproductive phenology in Spring (April-July) 41 . It has relatively large zygotes with distribution range less than one meter. However, the C. compressa thallus has aerocysts, enabling buoyancy and raft formation, facilitating the distribution of fertile branches over several kilometres 41 .
The distribution of Padina pavonica is conservatively considered to resemble C. nodosa and C. compressa, spanning throughout the Mediterranean and into the eastern Atlantic. P. pavonica was previously thought to have a global distribution, but molecular analysis of the genus has found no evidence to support this 15 . Instead it has been suggested that P. pavonica was potentially misclassi ed outside of the Mediterranean, due to morphological similarity with congeneric species 15 . Padina is a monophyletic genus with a worldwide distribution from tropical to cold temperate waters 15 . Most species have a regional distribution, with few con rmed examples of species spanning beyond a single marine realm sensu 42 . Recent molecular and morphological assessment of Padina in the Mediterranean Sea, identi ed two new cryptic species (P. ditristromatica and P. pavonicoides), formerly considered to be P. pavonica 43 . These species can only be distinguished morphologically by the number of cell layers making eld identi cation impossible. Given the close relatedness, shared distribution and di culty to distinguish between these species, we herein refer to them collectively as Padina pavonica, but acknowledge the possibility that the three species were used unwittingly in experiments. We considered the poleward distribution limit of P. pavonica to be the British Isles 50ºN 44 , although this identi cation has not been con rmed by molecular analysis.

Sample collection and acclimation and aquarium setup
Sample collections were conducted at two sites, separated by approximately 1 km, within each location. Collections were conducted between 1-3 m depth and were spaced across the reef or meadow to try and minimise relatedness between shoots or fragments. Upon collection, fragments were placed into a mesh bag and transported back to holding tanks in cool, damp, dark conditions following 5 . Fragments were kept in aerated holding tanks in the collection sites at ambient seawater temperature and maintained under a 14:10 light-dark cycle until transport back to Mallorca. Prior to transport, P. oceanica shoots were clipped to 25 cm length (from meristem to tip), to standardise initial conditions and reduce biomass for transport. For transport back to Mallorca, fragments were packed in layers within cool-boxes. Frozen coolpacks were wrapped in damp tea towels (rinsed in seawater) and placed between layers of samples. On arrival at the destination, samples were returned to holding tanks with aerated seawater and a 14:10 lightdark cycle.
Experimental design: thermal performance experiments Following 48 hrs under ambient (collection site) conditions, samples were transferred to individual experimental aquaria, which consisted of a double layered transparent plastic bag lled with 2 L of ltered seawater (60 µm) following 45 . 16 experimental bags were suspended within 80L temperaturecontrolled baths. In total, ten baths were used, one for each experimental temperature treatment. Bath temperatures were initially set to the acclimatization temperature (i.e. in situ temperatures) and were subsequently increased or decreased by 1°C every 24 hours until the desired experimental temperature was achieved. Experimental temperatures were: 15,18,21,24,26,28,30,32,34 and 36°C (Table S2). For each species, four replicate aquarium bags were used for each temperature treatment with three individually marked seagrass shoots or three algal fragments placed into each bag. For P. oceanica, each marked plant was a single shoot including leaves, vertical rhizome and roots. For C. nodosa, each marked individual consisted of a 10 cm fragment of horizontal rhizome containing three vertical shoots. Individually marked seaweeds contained the holdfast, and 4-5 fronds of P. pavonica (0.98 ± 0.06 g FW; mean ± SE) or a 5-8 cm fragment for C. compressa (3.67 ± 0.1 g FW; mean ± SE). Once the targeted temperatures were reached in all of the baths, experiments ran for 14 days for the algal species and 21 days for seagrasses to allow for measurable growth in all species at the end of the experiment. Experiments were conducted inside a temperature-controlled chamber at constant humidity and air temperature (15°C). Bags were arranged in a 4x4 grid within each bath, enabling four species/population treatments to be run simultaneously. Bags were mixed within each bath so that one replicate bag was in each row and column of the grid, to minimise any potential within bath effects of bag position. Replicate bags were suspended with their surface kept open to allow gas exchange and were illuminated with a 14h light:10h dark photoperiod through uorescent aquarium growth lamps. The water within the bags were mixed with aquaria pumps. The light intensity throughout each bag was measured via a photometric bulb sensor (LI-COR) and ranged between 180-258 µmol m − 2 s − 1 . The temperature in the baths was controlled and recorded with an IKS-AQUASTAR system, which was connected to heaters and thermometers. The seawater within the bags was renewed every 72 hrs and salinity was monitored daily with an YSI multi-parameter meter. Distilled water was added when necessary to ensure salinity levels remained within the range of 36-39 PSU, typical of the study region. Carbon and Nitrogen concentrations in the leaf tissue were measured at the end of the experiment for triplicates of the 24ºC treatment for each species and location (Fig. S2) at Unidade de Técnicas Instrumentais de Análise (University of Coruña, Spain) with an elemental analyser FlashEA112 (ThermoFinnigan).

Growth measurements and statistical analyses
Net growth rate of seagrass shoots was measured using leaf punching-technique 46 . At the beginning of the experiment seagrass shoots were punched just below the ligule with a syringe needle and shoot growth rate was estimated as the elongation of leaf tissue in between the ligule and the mark position of all leaves in a shoot at the end of the experiment, divided by the experimental duration. Net growth rate of macroalgae individuals was measured as the difference in wet weight at the end of the experiment from the beginning of the experiment divided by the duration of the experiment. Moisture on macroalgae specimens was carefully removed before weighing them. Patterns of growth in response to temperature were examined for each experimental population using a gaussian function: , where k = amplitude, µ = mean and σ = standard deviation of the curve. Best t values for each parameter were determined using a non-linear least squares regression using the 'nlstools' package 47 in R 48 . 95% CI for each of the parameters were calculated using non-parametric bootstrapping of the mean centred residuals. The relationship between growth metrics and the best-t model was determined by comparing the sum of squared deviations (SS) of the observed data from the model, to the SS of 10 4 randomly resampled datasets. Growth metrics were considered to display a signi cant relationship to the best-t model if the observed SS was smaller than the 5th percentile of randomised SS. Upper thermal limits were de ned as the optimal temperature + 2 standard deviations (95th percentile of curve) or where net growth = 0. Samples that had lost all pigment or structural integrity by the end of the experiment were considered dead and any positive growth was treated as zero.

Metabolic rates
Net production (NP), gross primary production (GPP) and respiration (R) were measured for all species from the four sites for ve different experimental temperatures containing the in-situ temperature during sampling up to a 6ºC warming (see SM Table S3 for details). Individuals of the different species were moved to methacrylate cylinders containing seawater treated with UV radiation to remove bacteria and phytoplankton, in incubation tanks at the 5 selected temperatures. Cylinders were closed using gas-tight lids that prevent gas exchange with the atmosphere, containing an optical dissolved oxygen sensor (ODOS® IKS), with a measuring range from 0-200 % saturation and accuracy at 25ºC of 1% saturation, and magnetic stirrers inserted to ensure mixing along the height of the core. Triplicates were measured for each species and location, along with controls consisting in cylinders lled with the UV-treated seawater, in order to account for any residual production or respiration derived from microorganisms (changes in oxygen in controls was subtracted from treatments). Oxygen was measured continuously and recorded every 15 minutes for 24 hours.
Changes in the dissolved oxygen (DO) were assumed to result from the biological metabolic processes and represent NP. During the night, changes in DO are assumed to be driven by R, as in the absence of light, no photosynthetic production can occur. R was calculated from the rate of change in oxygen at night, from half an hour after lights went off to half an hour before light went on (NP in darkness equalled R). NP was calculated from the rate of change in DO, at 15 min intervals, accumulated over each 24 h period. Assuming that daytime R equals that during the night, GPP was estimated as the sum of NP and R. To derive daily metabolic rates, we accumulated individual estimates of GPP, NP, and R resolved at 15 min intervals over each 24 h period during experiments and reported them in mmol O 2 m − 3 day − 1 . A detailed description of calculation of metabolic rates can be found at Vaquer-Sunyer et al. 49 .

Sea temperature measurements and reconstruction
Underwater temperature loggers (ONSET Hobo pro v2 Data logger) were deployed at collection sites in Catalonia, Mallorca, Crete and Cyprus and left to record hourly temperatures for 12-18 months between Page 13/17 (SST) have been used. In particular, we use daily SST maps with a spatial resolution of 1/4°, obtained from the National Center for Environmental Information (NCEI, https://www.ncdc.noaa.gov/oisst 50 ). These maps have been generated through the optimal interpolation of Advanced Very High Resolution Radiometer (AVHRR) data for the period 1981-2019. The closest points to the collection sites have been extracted from the satellite images and the resulting time series calibrated against the in-situ data.
Corrected-SST data was used to reconstruct daily habitat temperatures from 1981-2019.

Thermal distribution measurements
We estimated the realised thermal distribution for the four experimental species and the genus level distribution for each species by downloading occurrence records from the Global Biodiversity Information Facility (GBIF.org (11/03/2020) GBIF Occurrence Download). Occurrence records from GBIF were screened for outliers and distributions were veri ed from the primary literature 15,35,38,40,43,51 and Enrique Ballesteros (pers. comms) (Fig. S1). Mean, 1st and 99th percentiles of daily SST's were downloaded for each occurrence site for the period between 1981-2019 using the SST products described above (Table   S4). Thermal range position of species at each experimental site were standardised by their global distribution using a Range Index RI; 52 . Median SST at the experimental collection sites were standardized relative to the thermal range observed across a species realized distribution, using the equation Declarations manuscript with contributions from all authors.

Data availability statement
All data for this study have been deposited in Dryad (https://doi.org/10.5061/dryad.d2547d81r) and will be made publicly available upon publication.

Ethical approval
Collections of plant material for this project were conducted in accordance with institutional, national, and international guidelines and legislation. Permissions for collections of plant material were obtained from local governing bodies.