Resilience of seagrass populations to thermal stress does not reflect regional differences in ocean climate.

The prevalence of local adaptation and phenotypic plasticity among populations is critical to accurately predicting when and where climate change impacts will occur. Currently, comparisons of thermal performance between populations are untested for most marine species or overlooked by models predicting the thermal sensitivity of species to extirpation. Here we compared the ecological response and recovery of seagrass populations (Posidonia oceanica) to thermal stress throughout a year-long translocation experiment across a 2800-km gradient in ocean climate. Transplants in central and warm-edge locations experienced temperatures > 29°C, representing thermal anomalies > 5°C above long-term maxima for cool-edge populations, 1.5°C for central and < 1°C for warm-edge populations. Cool-edge, central and warm-edge populations differed in thermal performance when grown under common conditions, but patterns contrasted with expectations based on thermal geography. Cool-edge populations did not differ from warm-edge populations under common conditions and performed significantly better than central populations in growth and survival. Our findings reveal that thermal performance does not necessarily reflect the thermal geography of a species. We demonstrate that warm-edge populations can be less sensitive to thermal stress than cooler, central populations suggesting that Mediterranean seagrasses have greater resilience to warming than current paradigms suggest.

Introduction

Ocean warming is having profound effects on benthic marine ecosystems across the globe (Hoegh‐Guldberg & Bruno, 2010). Foundation species of seaweeds, seagrasses and corals can be particularly susceptible to warming and have undergone extensive thermal stress (Hughes et al., 2018; Smale, 2020), mortality (Marba & Duarte, 2010) and range contraction (Wernberg et al., 2016) over recent decades. At the same time, temperature effects on populations across a species geographical range are not homogeneous and depend on factors such as a population's thermal range position, exposure to warming, a species evolutionary history and adaptive capacity (Angilletta, 2009; Somero, 2010; Donelson et al., 2019).

Accurate predictions of the thermal sensitivity of species is critical to anticipate and manage the impacts of climate change. Currently, broad‐scale predictions of thermal sensitivity typically rely on metrics of thermal performance such as a species‐realised thermal niche or extrapolation of local thermal performance measures to other locations, as proxies of how close an organism may be living to its thermal limits. Both these approaches typically overlook the prevalence and importance of local adaptation or phenotypic plasticity in thermal performance between populations, despite having a potentially profound influence on when and where climate change impacts might occur (Donelson et al., 2019). For example, if individuals from different populations display similar optimal and upper thermal limits, then cool‐edge populations will have lower sensitivity to warming than warm‐edge populations in which ambient temperatures are already close to upper thermal limits (Bennett et al., 2019). This is an inherent assumption of species distribution models that rely on the realised distribution of species to predict the likelihood of extirpation or range shifts under climate change (Araújo & Peterson, 2012). If, conversely, thermal performance differs between populations – through local adaptation or phenotypic plasticity – then the thermal safety margin between ambient conditions and an individual’s upper thermal limit may not reflect a species thermal geography or remain relatively constant between populations (Bennett et al., 2015). Determining whether thermal performance varies between populations and how this relates to local climatology or a species distribution is critically important for our capacity to anticipate and manage climate change impacts. For the vast majority of species, within‐species patterns of thermal sensitivity remain unknown, highlighting a fundamental knowledge gap in climate change ecology.

In conjunction with direct physiological impacts, ocean warming is leading to changes in species distributions resulting in new assemblages of species coming together and competing for space and resources within an ever‐changing physical environment (Lenoir et al., 2020). As species redistribute around the world, novel interactions are emerging that may either accelerate or buffer ecosystems from change under different environmental contexts (Gilman et al., 2010; Lurgi et al., 2012). Incorporating this complexity into models is challenging based on current approaches (Pagès et al., 2018). For example, thermal performance studies are routinely conducted in controlled laboratory environments that are invaluable for understanding complex physiological and evolutionary processes, but often miss key ecological processes and trade‐offs that may determine an organism’s success in the wild (Buñuel et al., 2020). While the ecological realism of experimental systems is improving (Ullah et al., 2018), calibrating laboratory‐based experiments with real world multispecies ecological interactions remains a challenge. At the other end of the spectrum, direct impacts of thermal stress on natural ecosystems from marine heatwaves and long‐term warming offer invaluable insights into the whole‐of‐ecosystem response to warming (Smale et al., 2019). However, such events are irregularly distributed and historically uncommon (Oliver et al., 2018). In addition, the ecological consequences of marine heatwaves can be catastrophic, undermining the health of the ecosystem for which the information is needed (Wernberg et al., 2016).

In this study we took an intermediate approach and conducted a 12‐month translocation experiment of the seagrass Posidonia oceanica across its geographical range. Translocation experiments have the benefit of directly comparing the performance of geographically distant populations under common, relatively natural ecological settings. We conducted our study across a steep 5°C gradient in average annual temperatures from the western to eastern Mediterranean Sea. The Mediterranean has experienced rapid warming between 0.25 and 0.65°C dec−1 over the past 3 decades, 2–3 times faster than the average global ocean (Marbà et al., 2015). Marine heatwaves have impacted central populations of P. oceanica in recent decades (Diaz‐Almela et al., 2009; Marba & Duarte, 2010) and previous studies have predicted that the species could face functional extinction by 2050 (Jordà et al., 2012; Chefaoui et al., 2018). Empirical observations of natural meadows have reported 28°C as the upper threshold for P. oceanica survival, a finding supported by controlled experiments of thermal performance (Savva et al., 2018). However, these findings have primarily come from the western Mediterranean, where summer temperatures rarely exceed 28°C. By comparison, P. oceanica in parts of the eastern Mediterranean, regularly experience summer temperatures above 28°C (Chefaoui et al., 2017), similar to heatwave conditions attributed to losses in the western basin. High thermal sensitivity in western ‘central’ populations, therefore, raises questions about the sensitivity of the eastern ‘warm‐edge’ populations to temperature. The aim of this experiment was to compare whether warm‐edge populations are indeed more sensitive to climate warming than central and cool‐edge populations due to their thermal geography or whether adaptation or acclimatisation to local conditions modifies patterns of thermal sensitivity across a species range.

Materials and Methods

Experiment locations and climate

Trans‐Mediterranean translocation of P. oceanica fragments took place between Catalunya (Spain), Mallorca (Spain) and Cyprus in July 2018 and were monitored until July 2019 (Fig. 1). Sea surface temperature (SST) data for each transplant site were based on daily SST maps with a spatial resolution of ¼°, obtained from the National Center for Environmental Information (NCEI, https://www.ncdc.noaa.gov/oisst; Reynolds et al., 2007). These maps have been generated through the optimal interpolation of Advanced Very High Resolution Radiometer (AVHRR) data for the period 1981–2019. Underwater temperature loggers (ONSET Hobo pro v.2 Data logger, Onset Computer Corp., Bourne, MA, USA) were deployed at the transplant sites in Catalunya, Mallorca and Cyprus and recorded temperatures hourly throughout the duration of the experiment (1 yr). To obtain an extended time series of temperature at transplant sites, a calibration procedure was performed comparing logger data with SST from the nearest point on SST maps. In particular, SST data were linearly fitted to logger data for the common period. Then, the calibration coefficients were applied to the whole SST time series to obtain corrected‐SST data and reconstruct daily habitat temperatures from 1981 to 2019. Local climate data were also compared with the global thermal distribution of P. oceanica to assess how representative of the thermal distribution of the species were experimental sites (Supporting Information Methods S1). Collectively, seawater temperatures from the three locations spanned the 1st to 99th percentile of temperatures observed across the global thermal distribution of P. oceanica. Catalan, called here ‘cool‐edge’ populations, comes from the north‐western end of Posidonia’s distribution and fell within the 1st percentile of both winter minimums and summer maximum temperatures. Cyprus, called here ‘warm‐edge’ populations, are situated in the eastern end of Posidonia’s distribution and fall within the 99th percentile of winter minimums and summer maximum temperatures. Mallorcan, called here ‘centre’ populations, sit geographically toward the cooler edge in the western basin, but climatically toward the warmer edge, in the 80–90th percentiles of temperatures experienced across Posidonia’s distribution (Fig. 1).

Fig. 1

(a) Map of source and recipient transplant sites across the Mediterranean Sea. Colours indicate the source and direction of transplants in Catalunya (blue), Mallorca (yellow and orange) and Cyprus (red). (b) Illustrated example of a transplanted Posidonia oceanica fragment, characterised by an apical shoot and multiple vertical shoots. (c) Daily temperature anomalies experienced by P. oceanica transplants, with respect to the maximum annual temperature recorded in the site of origin (averaged between 1981 and 2019). Dashed vertical line, highlights switch between positive (stress) and negative (recovery) thermal anomalies. Maximum annual temperatures recorded in Catalunya, 25°C; Mallorca, 29.1°C; and Cyprus, 29.3°C. T0–T4 indicate the sampling period. Daily temperatures experienced by transplants relative to their thermal regimes of origin in (d) Catalunya, (e) Mallorca and (f) Cyprus. Shaded areas represent daily temperatures above thermal regime of origin. Different shades of blue, yellow, orange and red within (d–f), illustrate temperature profiles in the different recipient sites. Grey line represents the long‐term daily temperature, averaged between 1981 and 2019. SST, sea surface temperature.

The focus of the experiment was to examine differences in thermal performance of P. oceanica populations in the context of climate warming and marine heatwaves. Transplantation therefore took place toward warmer climates and procedural controls were conducted within each source location, resulting in six source‐to‐recipient treatments (Fig. 1). Reciprocal transplants toward cooler climates were not performed, and therefore potential differences in performance in relation to cool limits are beyond the scope of this study. Initial collection of P. oceanica, handling and transplantation was carried out simultaneously by coordinated teams in July 2018 (Table S1). Each recipient location was subsequently resampled four times throughout the experiment, in August–September 2018 (T1), October 2018 (T2), April 2019 (T3) and May–June 2019 (T4; Table S1).

Between 60 and 100 fragments were collected for each of the six treatments (i.e. cool–cool, cool–centre, cool–warm, centre–centre, centre–warm and warm–warm). A fragment was defined as a section of P. oceanica containing one apical shoot connected with c. five vertical shoots and 10–15 cm of rhizome with intact roots. Collection occurred at two sites within each location, separated by c. 1 km. Within sites, collections were conducted between 4 and 5 m depth and were spaced across the meadow to minimise the dominance of a single clone and damage to the meadow. Upon collection, fragments were transported in shaded seawater for up to 1 h back to the nearest laboratory.

Handling methods

In the laboratory, fragments were placed into holding tanks with aerated seawater, at ambient temperature and a 14 h : 10 h, light : dark cycle. All shoots were clipped to 25 cm length (from meristem to the tip of the longest leaves), to standardise initial conditions and reduce biomass for transportation. For transport by plane or ferry between locations, fragments were packed in layers within cool boxes. Each layer was separated by frozen cool packs wrapped in wet tea towels (rinsed in sea water). All fragments spent 12 h inside a cool box irrespective of their recipient destination, including procedural controls (i.e. cool–cool, centre–centre and warm–warm) to simulate the transit times of the plants travelling furthest from their source location (Fig. 1a). On arrival at the destination, fragments were placed in holding tanks with aerated seawater at ambient temperature as described above in their recipient location for 48 h, before field transplantation.

Measurement methods

At 1 d before transplantation, fragments were tagged with a unique number and the rhizome of each fragment was cable tied to a U‐shaped peg. Morphological traits for each fragment were measured and included: (1) length of the longest apical leaf, width and number of leaves, (2) total number of bite marks on leaves of three vertical shoots per fragment, (3) number of vertical shoots, (4) leaf count of three vertical shoots per fragment, and (5) overall horizontal rhizome length. A subset (n = 10) of fragments per treatment was marked prior transplantation to measure shoot growth. To do this, all shoots within a single fragment were pierced using a hypodermic needle. Two holes were pierced side‐by‐side at the base of the leaf/top of the meristem.

Transplant methods

All transplant sites were located at a 4–5 m depth in areas of open dead matte surrounded by P. oceanica meadow. In Mallorca and Cyprus, fragments were distributed between two sites, separated by c. 1 km. In Catalunya, the lack of suitable dead matte habitat meant that all fragments were placed in one site. Pegs, adjoined to each fragment were pushed in by hand, ensuring no damage to the plant and leaving the rhizome and roots in direct contact with the substrate. Fragments were planted along parallel transects at 50‐cm intervals and with a 50‐cm gap between parallel transects (Fig. S1). Different treatments were mixed and deployed haphazardly along each transect.

Resampling methods and herbivory

On day 10 of the experiment, a severe herbivory event was recorded at both warm‐edge translocation sites. Scaled photographs of all fragments were taken at this time to record the effects of herbivory on transplants. At the end of each main sampling period (T0–T1, T1–T2 and T3–T4), all pierced fragments were collected and taken back to the laboratory to measure shoot growth. At T1, T2 and T3, additional sets of fragments (n = 10 per treatment) were marked using the piercing method to record growth in the subsequent time period. In addition, at T1 and T3, n = 20 shoots within the natural meadow at each site were marked to compare growth rates between the native meadow and transplants. Underwater shoot counts and a scaled photograph were taken to record fragment survivorship, shoot mortality, bite marks and shoot length among all remaining fragments within each site and sampling time.

In the laboratory, morphological measurements (See Measurement Methods) were repeated on the collected fragments and growth of transplant and natural meadow shoots was measured. Growth (shoot elongation, cm d−1) of the marked shoots was obtained by measuring the length from the base of meristem to marked holes of each leaf (new growth) of the shoot and dividing the leaf elongation per shoot by the marking period (in days). For each shoot, total leaf length (cm per shoot) and the number of new leaves was recorded. The rate of new leaf production (new leaves per shoot d−1) was estimated by dividing the number of new leaves produced per shoot and the marking period. New growth was dried at 60°C for 48 h to determine carbon and nitrogen content of the leaves, and carbon to nitrogen (C : N) ratios. Carbon and nitrogen concentrations in the new growth leaf tissue were measured at the beginning of the experiment and at each subsequent time point for each treatment. Nutrient analyses were conducted at the Unidade de Técnicas Instrumentais de Análise (University of Coruña, Spain) with an elemental analyser FlashEA112 (ThermoFinnigan, San Jose, CA, USA).

Underwater photographs of shoots were analysed using ImageJ software (https://imagej.net). Maximum leaf lengths on each shoot in warm‐edge transplant sites (cool–warm, centre–warm and warm–warm) were recorded for the initial (day 10) herbivore impact, T1, T2 and T3 and related to transplant nutrient concentrations. Herbivore impact was estimated as the proportional change in length of the longest leaf relative to initial length at T0.

Thermal stress

Long‐term maximum temperatures were recorded as the average of annual maximum daily temperatures in each transplant site, averaged between years from 1981 to 2019. Maximum thermal anomalies were calculated as the difference between daily temperatures in a recipient site over the course of the experiment and the long‐term maximum temperature in the source site for each corresponding population. ‘Heat stress’ and ‘recovery’ growth periods of the experiment were defined as T0–T2 (July–October) and T2–T4 (November–June), respectively, corresponding to periods of positive and negative maximum thermal anomalies. Thermal anomalies experienced by the different transplant treatments were plotted using the geom_flame function in the heatwavesr package (Schlegel & Smit, 2018) in R (v.3.6.1, 2019).

Statistical analyses

Differences in fragment survivorship between source–recipient treatments were analysed using a chi‐squared (χ) test at the end of the heat stress period (T2) and the end of the recovery period (T4). Differences in (1) Shoot mortality per fragment, (2) herbivory impact, (3) C : N of leaves, (4) growth, (5) new leaf production and (6) total leaf length were each compared between treatments (n = 6) within discrete time periods using analysis of variance (ANOVA). Normality was tested using the Shapiro–Wilks test and Q–Q plots, and homogeneity of variance was tested using Levene’s test and residuals plot. No transformations were applied. Tukey’s honest significant difference (HSD) analysis was used to examine differences between groups for significant main effects. To assess the relationship between thermal performance, source location and (1) maximum absolute realised temperature and (2) maximum thermal stress anomalies, linear mixed effects models were applied using the lme4 package in R (Bates et al., 2015). Source location and temperature metric were treated as fixed effects and ‘recipient location’ and ‘time point’ were treated as random effects, to account for differences in performance across locations and seasonality. Separate models were conducted for each temperature metric (i.e. maximum temperature (Tmax) or heat stress anomaly) and performance metric (i.e. growth, new leaf production and total leaf length) combination. We started with a full model based on a priori hypotheses on the inclusion of terms and higher order interactions. We identified the minimum adequate model by a stepwise removal of nonsignificant terms using likelihood ratio tests of the model with the effect in question against the model without the effect in question. Visual inspection of residual plots did not reveal any obvious deviations from homoscedasticity or normality and therefore no transformations were conducted. P‐values were obtained by likelihood ratio tests of the full model with the effect in question against the model without the effect in question. All analyses were conducted in R (v.4.0.3, 10 October 2020).

Results

Realised thermal anomalies

The western Mediterranean Sea experienced marine heatwave conditions during the 2018 summer, coinciding with the translocation experiment. Maximum daily cool‐edge temperatures were 24.9°C, 1.3°C above long‐term summer maximum (averaged 1981–2019). Maximum daily temperatures in central sites were 29.1°C, also 1.3°C above long‐term summer maximum and the second highest temperature on record. Maximum daily temperatures at P. oceanica’s warm‐edge in the eastern Mediterranean were more typical at 29.3°C and 0.18°C above long‐term maximum.

As a result of the marine heatwave in the western Mediterranean, thermal anomalies experienced by cool‐edge and central transplants were similar in central and warm‐edge locations. The highest anomalies were experienced by cool–centre and cool–warm fragments, which experienced maximum anomalies of 5.4°C and 5.7°C above long‐term maximum temperatures in the cool edge, respectively. Positive thermal stress anomalies lasted for 130 d for cool–warm fragments and 102 d for cool–centre fragments over the course of the experiment (Fig. 2; Table S1). Centre–warm and centre–centre fragments experienced maximum anomalies of 1.53°C and 1.3°C, respectively, for a duration of 37 d and 33 d above long‐term summer maxima in P. oceanica’s range centre. Warm–warm transplants experienced equally high temperatures to cool–warm and centre–warm transplants, but the smallest relative anomaly (0.18°C for 2 d) above long‐term summer maxima. The coolest daily temperatures recorded throughout the course of the experiment ranged from 10.5°C in the cool edge, 13.9°C in the centre and 16.2°C in the warm‐edge.

Fig. 2

Patterns of shoot length and nutritional quality of seagrasses over time. (a) Shoot length patterns in August 2018 highlight the effects of a selective feeding event in Cyprus where herbivores selectively overgrazed cool–warm transplants (marked by an asterisk). (b) Carbon : nitrogen ratios of Posidonia oceanica leaves in the different treatments throughout the course of the experiment. C : N ratios in cool–cool plants (blue dashed line) were only measured in July 2018 and April 2019. (c) Differences in herbivory on cool–warm transplants compared with the adjacent centre–warm and warm–warm transplants.

Selective herbivory in transplants with higher nutritional quality

At 10 d after initial deployment, transplants were surveyed revealing selective overgrazing on cool–warm P. oceanica (Fig. 2). Shoot lengths of cool–warm transplants were reduced by 74.5 ± 4% across both warm‐edge transplant sites, from 22.3 ± 3 cm to 5.7 ± 9 cm. By contrast, centre–warm and warm–warm plants that were interspersed with cool–warm plants, lost significantly less biomass, from 25 ± 3% and 27.9 ± 5% of initial length, respectively (ANOVA, F 2,194 = 37.64, P < 0.001; Fig. 2a). Selective grazing patterns were consistent between two warm‐edge sites separated by 1.6 km. Selectivity patterns in grazing rates reflected differences in nutritional quality between seagrass populations and regional differences in herbivore assemblages. The C : N ratio of cool‐edge P. oceanica leaves was 18.1 ± 1.4 (mean ± 1SE), significantly lower (i.e. higher nutritional quality) at the start of the experiment than C : N ratios for central 34.7 ± 1.4 or warm‐edge 52.3 ± 1.8 populations (ANOVA, F 2,8 = 97.76, P < 0.001; Fig. 2b). The morphology of the bite scars observed on the consumed seagrass in Cyprus suggested that herbivorous fish (Siganus rivulatus, S. luridus, Sparisoma cretense) and potentially sea turtles (Chelonia mydas) were responsible for the intense grazing rates, although relative contributions of each species could not be quantified (Fig. S2).

By week 6, significant differences in shoot length were still observed between treatments in Cyprus (ANOVA, F 2,245 = 87.69; P < 0.001), but differences between treatments were reduced. Similarly, the C : N ratio in leaf tissue of the transplants started to equilibrate with local conditions and by October the C : N ratio was not significantly different between treatments within a recipient location (Fig. 2b). Transplants in cool‐edge and central locations displayed no evidence of overgrazing by herbivores in September and were significantly longer than all warm‐edge treatments (ANOVA, F 5,311 = 139.1, P < 0.001) (Fig. 2a). By the end of the heat stress period in October, however, all six treatments had reduced shoot lengths, with only subtle differences between them (Fig. 2a). Surviving leaves were, however, intact with no evidence of overgrazing by herbivores, suggesting that seasonal shedding, not herbivory was responsible for losses at cool‐edge and centre sites (Fig. 2a).

Survivorship

Despite experiencing extreme thermal anomalies and grazing rates at the start of the experiment, cool–warm transplants displayed high survivorship (number of surviving rhizome fragments) with no significant differences observed between treatments at the end of the thermal stress period (χ2 = 10.65, df = 5, P = 0.06; Fig. 3). Survivorship then declined in all treatments over winter months (between October and April), potentially due to dislodgement from storms. Notably, however, both centre–centre and centre–warm treatments displayed significantly lower survivorship than all cool‐edge and warm‐edge populations by the end of the experiment (χ2 = 58.07, df = 5, P < 0.001; Fig. 3a).

Fig. 3

Patterns of survivorship (a) and shoot survival (b) of Posidonia oceanica transplants over the course of the experiment. Transplant fragment survival is represented as the proportion of transplants remaining at the end of each growth period, relative to the expected number of transplants (i.e. initial conditions minus the fragments removed for growth measurements). Shoot survivorship is the mean change in the number of shoots per transplant (±1SE), relative to initial conditions. Blue circles and whiskers represent the fragments from Catalunya (cool edge), orange circles and whiskers represent fragments from Mallorca (central) and red represents fragments from Cyprus (warm edge).

The number of living shoots per fragment of rhizome was relatively stable over time, with slight declines over the period of thermal stress (Fig. 3b). Warm–warm and cool–centre treatments displayed significantly lower rates of shoot mortality than central and other cool‐edge treatments at the end of the thermal stress period (ANOVA, F 5,388 = 8.118, P < 0.001). By the end of the recovery period, changes in shoot numbers on the surviving fragments had stabilised and cool–centre fragments had recovered to display no net change from the beginning of the experiment.

Growth patterns

The number of leaves per shoot remained relatively stable through time and retained an average of 4.8 ± 0.13 leaves in cool–cool treatments, 4.3 ± 0.75 in centre–centre treatments and 4.3 ± 0.16 leaves per shoot in warm–warm treatments by the end of the experiment. Centre–warm treatments were the exception and declined from 4.0 ± 0.29 leaves per shoot at the beginning, to 2.5 ± 0.79 leaves per shoot by the end of the experiment. Underlying this pattern of stability, however, were marked differences in leaf production and turnover between treatments. New leaf production displayed a positive relationship with maximum temperature and significant differences between source populations, characterised by a higher rate of new leaf production in treatments sourced from the cool‐edge population, than the centre or warm‐edge populations (Table S2; Fig. 4a). In relation to heat stress anomalies, new leaf production displayed a strong positive relationship and no difference between source populations (Table S2; Fig. 4b), suggesting local thermal adaptation or acclimatisation within source populations. P. oceanica growth rates (shoot elongation) and total leaf length per shoot displayed nonsignificant relationships with maximum temperatures and heat stress anomalies (P > 0.05; Fig. 4d–f,h). More probably, source populations were the best predictor of growth and total leaf length over the course of the experiment (Table S2). At the end of the thermal stress period, growth rates were significantly different among transplant treatments (ANOVA, F 5,112 = 3.8, P = 0.003; Fig. 4f). Cool–cool and warm–warm treatments displayed the highest growth rates despite experiencing up to 6°C difference in temperature, and centre–centre and centre–warm treatments displayed significantly lower growth rates (Fig. 4f). By the end of the recovery period, growth rates had increased in all treatments except centre–warm shoots that remained significantly lower than cool‐edge (Tukey’s HSD, P = 0.005) and warm‐edge controls (Tukey’s HSD test, P = 0.018).

Fig. 4

Thermal performance of Posidonia oceanica transplants across the Mediterranean Sea. The y‐axes indicate performance metrics of new leaf production (a–c) and shoot elongation (d–f). The x‐axes indicate temperature metrics of absolute maximum temperature (a, d), maximum thermal stress anomaly (b, e) and a comparison between heat stress period (July–October 2018, grey) and recovery periods (November–June 2019, white) (c, f). Boxes and whiskers indicate mean ±  1SE rates per treatment. The shape of the box illustrates the time period in which measurements were recorded. The colour of the box represents the different transplant treatments and procedural controls as per previous figures.

Growth rates in natural meadows displayed similar patterns to transplants, albeit with higher rates over the same heat stress and recovery periods (Fig. S3). During the thermal stress period, growth rates were significantly higher in cool‐edge meadows than warm‐edge meadows and both cool‐edge and warm‐edge meadows displayed significantly higher growth rates than central meadows (ANOVA, F 2,66 = 17.12, P < 0.001). New leaf production was highest in central meadows and significantly lower in cool‐edge and warm‐edge meadows during the thermal stress period (ANOVA, F 2,66 = 13.56, P < 0.001).

Discussion

Improving predictions for how species and communities will be affected by climate change and extreme events is an urgent challenge. Here we examined how populations of the functionally important seagrass P. oceanica could be affected by heat stress, by conducting a translocation experiment across its range. Cool‐edge, central and warm‐edge populations displayed clear differences in thermal performance when grown under common conditions, but patterns contrasted with a priori expectations based on their thermal geography. Cool‐edge populations performed equally well to warm‐edge populations under common conditions and significantly better than central populations in terms of growth rates and survival. The strong performance of cool‐edge populations was remarkable given they experienced thermal anomalies of > 5°C above long‐term summer maxima, compared with 1.5°C for central populations and < 1°C for warm‐edge populations. Moreover, cool‐edge transplants suffered and subsequently recovered from severe overgrazing at the warmest site that removed > 75% of aboveground biomass, highlighting marked resilience of cool‐edge populations to combined effects of thermal stress and herbivory, two pervasive impacts of climate change in marine ecosystems.

Effects of thermal stress on seagrass populations

Patterns of thermal performance in P. oceanica varied among populations and deviated from expectations with respect to range position, highlighting the challenges of extrapolating climate change vulnerability predictions from local observations. Among cool‐edge, central and warm‐edge populations, only measures of new leaf production supported models of local adaptation/acclimatisation and no metrics consistently supported expectations of thermal‐niche conservatism between populations. The positive relationship between new leaf production and thermal anomalies indicates higher rates of leaf turnover (i.e. shorter lifespan) under warmer conditions, potentially reflecting a coping strategy for thermal stress within fragments. Cool‐edge and warm‐edge populations did display similar rates of growth and survival under warm‐edge (29°C) conditions, suggestive of niche conservatism, however central populations exhibited comparatively poor performance, in contrast with expectations. Moreover, patterns of poor thermal performance by central populations at the warm edge were consistent with results among centre–centre procedural controls, which experienced similar absolute temperatures to warm‐edge treatments. Growth patterns from the natural meadows over the same period further supported these findings, showing lower growth rates and higher new leaf production in centre populations relative to cool and warm‐edge meadows following the thermal stress period, reinforcing the idea that experimental findings were not an experimental artefact.

Unlike centre populations, the high resilience of cool‐edge P. oceanica to thermal stress was consistent with documented responses of congeneric Posidonia species to a natural heatwave of similar magnitude. In 2011, Western Australia experienced the most intense marine heatwave on record (Hobday et al., 2016) and Posidonia australis displayed relatively minor impacts despite experiencing thermal anomalies in excess of 5°C and absolute temperatures reaching 30°C (Strydom et al., 2020). In contrast with P. australis, other seagrasses (i.e. Amphibolis antarctica) experienced severe losses (Fraser et al., 2014) and kelp forests (Ecklonia radiata) suffered a 100‐km range contraction (Wernberg et al., 2016) in response to the 2011 heatwave. While the response by P. oceanica to large anomalies in the current study was relative to cool‐edge (cf. warm‐edge in P. australis) conditions, the similarity between species to survive temperatures up to 30°C highlights a marked latent capacity to withstand thermal stress in these edge populations. Latent tolerance to extreme events potentially reflects the warm evolutionary origins of the genus and the high resilience that Posidonia has displayed to climatic upheaval over the past 60 million years (Myr) (Por, 2009; Bianchi et al., 2012).

The idea of P. oceanica as resilient to heat stress may appear at odds with the prevailing literature, which has demonstrated (Diaz‐Almela et al., 2009; Marba & Duarte, 2010) and projected (Jordà et al., 2012; Chefaoui et al., 2018) high sensitivity of P. oceanica to warming. However, the results of the current study are consistent with these previous findings and help to place them in a broader geographical context. Central populations of P. oceanica were severely affected by thermal anomalies up to 1.5°C above long‐term summer maxima in the current study. Anomalies of this magnitude are slightly cooler than those experienced during the 2003 Mediterranean heatwave, which caused a steep increase in shoot mortality in central populations (Marba & Duarte, 2010). Our findings therefore support previous evidence that P. oceanica in central populations is living close to its thermal limits, but suggest that populations elsewhere may indeed have greater tolerance of high temperatures. Cool‐edge P. oceanica, for example, had relatively strong survival and recovery from conditions > 5°C above the maximum summer temperatures that it typically experiences in its natural environment, suggesting a comparatively low sensitivity of cool‐edge populations to extirpation in coming decades.

While field‐based evidence of performance could not separate cool‐edge and warm‐edge populations, the maximum realised temperature in our experiment was 29.3°C. Therefore, we cannot rule out the possibility that differences in performance between populations occur at higher temperatures. Indeed, laboratory‐based experiments have suggested that warm‐edge populations can survive temperature up to 36°C, whereas cool‐edge populations and central populations of P. oceanica display upper thermal limits of c. 30°C (Bennett et al., In review). Such experimental evidence would suggest that cool‐edge populations were close to their upper threshold during the current experiment, whereas warm‐edge populations retained a larger thermal buffer. Previous laboratory experiments, however, did not document the subsequent recovery of P. oceanica following thermal stress. While cool–warm and cool–centre treatments clearly showed signs of thermal stress in the current study, our results suggested that they have greater capacity to endure and recover from heat stress than central populations.

The reason why central populations of P. oceanica are more severely affected by high temperatures than cool and warm‐edge populations is unclear. The genetic structure of P. oceanica is characterised by strong separation between the eastern and western basins, consistent with a vicariance event during the last glacial maxima (Arnaud‐Haond et al., 2007). Within the western basin, genetic connectivity is relatively low for both the cool‐edge and central experimental collection locations, but both these locations are directly connected to the same central node within the regional metapopulation (Rozenfeld et al., 2008). Central and cool‐edge populations, therefore, share more genetic and climatological similarities over geological timescales than cool‐edge and warm‐edge populations (Chefaoui et al., 2017), whereas central and warm‐edge populations share more contemporary similarities in terms of selection pressure on upper thermal limits. Despite this, central populations exhibited the greatest sensitivity to temperature, across the distribution of P. oceanica forcing a re‐evaluation of the vulnerability and management of this species in response to climate change.

Selectivity and recovery of seagrass from herbivory

In addition to direct impacts of temperature on P. oceanica, high rates of grazing by herbivores on cool–warm transplants highlight the importance of indirect impacts of warming on ecosystem function. The selective patterns of grazing at warm‐edge sites were consistent with the contrasting patterns of nutritional quality of transplants from different locations. These results are consistent with previous studies that have demonstrated higher rates of herbivory on nutrient‐enriched plants (Prado & Heck, 2011; Campbell et al., 2018). Interestingly, increased grazing on nutrient‐rich seagrass was only observed in warm‐edge, not central sites, where low nutrient meadows were also present. This pattern was probably due to differences in herbivore assemblages between locations. Tropical range extending fishes such as the rabbitfish S. rivulatus and S. luridus are the dominant herbivores in Cyprus (J. Santana‐Garcon et al., pers. comm.) but are currently absent from central and cool‐edge locations in the western Mediterranean (Daniel et al., 2009). Feeding scars on cool–warm transplants suggested that these warm‐affiliated fishes and sea turtles were responsible for the observed overgrazing in Cyprus. While overgrazing of transplants by herbivores was not observed in centre and cool‐edge locations, herbivory pressure by cool‐affiliated species in the western Mediterranean is already relatively high (Pagès et al., 2012). Range extension by tropical herbivores into the western Mediterranean could, therefore, exacerbate the impacts of herbivory on nutrient‐rich seagrasses in the region and needs to be monitored carefully.

Conclusions

The combination of warming temperatures, extreme events and species redistributions has the potential to fundamentally change marine ecosystems in the Anthropocene. Identifying general principles, developing models and testing predictions about how these impacts will unfold across species distributions is critical to anticipate and adaptively manage ecosystems over the coming decades. Here we demonstrate that cool‐edge, central and warm‐edge populations of seagrass differ in their response to temperature, however these responses are not necessarily linear or conform to a priori expectations. Moreover, the highly selective effects of warm‐affiliated herbivores on cool‐edge populations, highlight how species redistributions and species interactions can have unexpected and consequential effects on ecological performance under climate change. Our findings make a strong case for ongoing empirical comparisons of thermal performance for functionally important species and highlight the complexity of inferring thermal sensitivity from a single location. In doing so these findings provide hope for greater resilience to warming in Mediterranean seagrasses than previously recognised.

Author contributions

SB and NM conceived the study. SB, TA, DK, CA, JB, XB, PK, GR, JS‐G, IS, AV and NM conducted field work. LC and SB analysed the herbivory data. GJ and SB analysed the temperature data. SB analysed the transplant data and wrote the manuscript with contributions from all authors.

Fig. S1 Example site map of transplant experimental setup.

Fig. S2 Different types of bite scars observed on seagrass leaved following herbivory event in Cyprus.

Fig. S3 Growth rates, new leaf production and total leaf length of Posidonia oceanica shoots growing in adjacent natural meadows at the time of the transplant study.

Methods S1 Thermal distribution measurements.

Table S1 Summary of recipient locations and sampling times when measurements on Posidonia oceanica transplants were made.

Table S2 Results of linear mixed effects model for new leaf production, shoot growth and total leaf length.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.

Click here for additional data file.