Trait gradients inform predictions of seagrass meadows changes to future warming

Comparing populations across temperature gradients can inform how global warming will impact the structure and function of ecosystems. Shoot density, morphometry and productivity of the seagrass Posidonia oceanica to temperature variation was quantified at eight locations in Sardinia (western Mediterranean Sea) along a natural sea surface temperature (SST) gradient. The locations are spanned for a narrow range of latitude (1.5°), allowing the minimization of the effect of eventual photoperiod variability. Mean SST predicted P. oceanica meadow structure, with increased temperature correlated with higher shoot density, but lower leaf and rhizome width, and rhizome biomass. Chlorophyll a (Chl-a) strongly impacted seagrass traits independent of SST. Disentangling the effects of SST and Chl-a on seagrass meadow shoot density revealed that they work independently, but in the same direction with potential synergism. Space-for-time substitution predicts that global warming will trigger denser seagrass meadows with slender shoots, fewer leaves, and strongly impact seagrass ecosystem. Future investigations should evaluate if global warming will erode the ecosystem services provided by seagrass meadows.

Understanding future warming effects on foundation species, as marine macrophytes, is pivotal to predict their distribution and physical structure 40 , as temperature is thought the most important range limiting factor 41 . Seagrasses are valuable providers of coastal ecosystem services including, carbon sinks, nursery grounds, habitat, nutrient cycling, sediment stabilization, trophic transfer to adjacent habitats [42][43][44] and protection from erosion 45,46 . Posidonia oceanica (L.) Delile is a slow-growing seagrass, endemic to the Mediterranean, experiencing widespread decline due to multiple local anthropogenic stressors 47 . The abrupt decline experienced by P. oceanica from recent heatwaves 48 , however, has seriously questioned its persistence for the coming decades 40 . Due to its vulnerability in aquaria and slow growth, laboratory experiments have been limited and controversial. Nevertheless, plants from warm thermal environments were found to activate a suite of physiological 49 and molecular mechanisms [50][51][52] to tolerate simulated heatwave exposures, whereas phenological response to warming likely involves higher flowering 53 and denser meadows 54 . This is a space-for-time substitution, a method for studying slow ecological processes, where the relationships between ecological variables are studied at sites that are assumed to be at different stages of development 55 . This study is based on the assumption that plant functional traits vary along environmental gradients and potentially predict responses to environmental change. Thus, to examine the performance of P. oceanica to future temperature conditions, we measured shoot density, morphometry and productivity at eight locations in Sardinia (western Mediterranean Sea) along a natural gradient of water temperature. Despite similar latitude (minimum interference of photoperiod), the western locations are generally cooler than the eastern sites, with differences in SST comparable to climate change scenarios for the twenty-first century for the Mediterranean Sea (peaking at 2.6 °C in 2100 56,57 ), making this space-for-time substitution informative for projections of trait changes over the next decades. Chl-a, a proxy of light irradiance, was a further driver of seagrass structure. P. oceanica is currently in the EU Marine Strategy Framework Directive monitoring protocols 58 .

Results
Seagrass variability. Shoot density changed considerably between Sardinian coasts (Table 1 and Fig. 1) as well as leaf width which was larger on the west than on the east side (Table 1 and Fig. 2), although both variables were significant across locations and areas. All other morphometrical variables were significantly affected by location and area, except for necrotic leaf portion that was only dependent on the area (Table 1 and Fig. 2).
The reconstruction analysis showed that annual plant productivity changed between coasts only in terms of number of scales (remnant leaf sheats) and rhizome width, being lower on the east coast. Rhizome width and biomass were significantly dependent on the location, while all other variables were highly area dependent (Table 1 and Fig. 2).

Relationship between seagrass and environmental variables.
Multiple regressions retained only mean temperature in four models indicating that leaf width, number of scales, rhizome width and rhizome biomass were negatively related with mean temperature. Shoot density was related to mean temperature and Chl-a, as well as the number of scales and rhizome width (Table 2, Figs. 3 and 4). Specifically: I) increased shoot density was correlated with increased mean temperature, while an opposite trend was found for the leaf width, number of scales and rhizome width (Fig. 3) and II) reduced shoot density, number of scales and rhizome width were correlated with increased Chl-a (Fig. 4). The response variables where models retained Chl-a as the only explanatory variable, were the number of leaves and rhizome length, which increased and decreased, respectively, with increasing Chl-a (

Discussion
Posidonia oceanica morphometry and productivity were linked to the thermal environment. Increased temperature triggered higher shoot density, but lower leaf and rhizome width, fewer scales and lower rhizome biomass. Additionally, Chl-a was a temperature independent driver of the plant performance. Temperature strikingly affected shoot density, increasing gradually across the thermal gradient from 496.1 ± 21.6 to 829.9 ± 43.2 shoots/ m 2 (mean ± SE n = 12) at AHO and REI, respectively. Shoot density is the most common descriptor of P. oceanica meadows defining its conservation status (Marine Strategy Framework Directive) assuming that higher densities reflect lower human influence and better marine water conditions. However, the density classes distinguished by previous authors (reviewed by 59 ), ignore natural environmental variation. Since our data were collected unaffected from local anthropogenic disturbances, our results highlight that thermal environment is critical factor in determining plant shoot density, providing evidence of the need to refer the seagrass density classes to the mean temperature environment.
Our results revealed a strong spatial association between plant traits and temperature across a gradient suggest that future warming is predicted to produce denser P. oceanica meadows. This finding is corroborated by long-term correlative data revealing that shoot density is a plant trait that varies with thermal environment 54 , providing evidence that the plant would rearrange (increasing the number of modules) the meadows structure with warming (Fig. 5). The fact that Chl-a is inversely related to the meadow density will sharpen this pattern, as this influence is disentangled from temperature effects and because both drivers work in the same direction, enhancing shoot density and potentially producing synergistic effects. In fact, numerical models of future Chl-a due to anthropogenic climate change, generally suggest a decrease in globally integrated primary productivity driven by a reduction in supply of macronutrients [60][61][62][63] . Nevertheless, predicting meadow structure based on the relationship between spatial pattern of plant traits and the environment assumes that the seagrass traits could change proportionally to climate change 54 , although the species may respond to finer-scale changes in environmental variables that cannot be predicted using averages 64,65 .
Regarding mechanisms regulating the Chl-a-shoot density interaction, our data support the hypothesis that different light conditions due to the phytoplankton density (not nutrient availability) are involved, although manipulative experiments are needed. In fact, evidence of reduction of P. oceanica shoot density with depth are commonly gained [66][67][68] , supporting the hypothesis that light extinction is pivotal 69 . observed that seagrasses growing in low light reduce shoot density and above-ground biomass as an acclimation response to reduce selfshading within the canopy.
Shoot density changes induced by the climate change, however, will involve other phenological traits, such as leaf width and number of leaves. Their dependence on shoot density has been interpreted as the result of selforganization to shading [70][71][72][73][74] . Reducing the size of ramets to attenuate intraspecific competition is a common pattern in clonal plants 73,74 . Productivity of P. oceanica was not directly dependent on shoot density, but it seems that it will be contrastingly affected by the temperature and Chl-a, so that predicting the number of scales and rhizome width in coming decades is not obvious and likely dependent on the strength of their associations. Therefore, the prediction about the productivity that can be made on the trait gradients (trait variation along environmental gradient) regards the decrease in rhizome biomass and length affecting the plant robustness through decades.
Future changes in temperature and Chl-a, may drive P. oceanica morphometry and productivity patterns that will affect the ecosystem services that seagrass meadows currently provide. Quantification of seagrass services, however, have never been provided on a structure-specific basis 75,76 and we believe this might become a relevant issue. Indeed, in a future warmer Mediterranean Sea, where summer mean SST increase will likely peak 2.9 °C and 2.7 °C for the end of the century on the east and west Sardinia coasts, respectively 56 , P. oceanica leaf canopy, structured by higher shoot density with bundles of a lower number of leaves smaller in width, can create a different habitat and associated community. Similarly, whether the reduction in rhizome width and biomass has consequences on both the vulnerability of plants to storms and Carbon storage remains unanswered. www.nature.com/scientificreports/ This study shed light on how seagrass systems could respond to climate change, independently of the effects of extreme events (such as heat waves), as the latter undoubtedly affect deleteriously the seagrass structure with die-offs 47,77,78 . Nevertheless, the extent the phenotypic gradients of the seagrass systems depend on acclimation versus adaptation processes should be measured. However, the analysis of processes involved in phenotypic plasticity and the possibility that such plastic responses might be adaptive is complex for both the long-life cycles and slow growth of most of the seagrasses that impede manipulative experiments and trans-generation assessments 79 . Further space-for-time substitutions to predict functional traits changes due to global warming in seagrasses are necessary. Future trait gradients analysis should consider wider thermal range to sharpen our prediction and establish how closely the highest mean temperature used in the model stands are to the tolerance limit of the seagrass.   Fig. 6) where differences in water conditions are evident. The western coastline receives Atlantic waters directly through the Western Mid-Mediterranean Current and is also influenced by coastal upwellings 80 . In contrast, the eastern coast is affected by the warm Algerian Current 81 . Seagrass meadows unaffected from local anthropogenic disturbances (e.g. harbour, fish farming, and urbanisation) were sampled in eight different locations (Fig. 6), with a hierarchical design: for both coasts of Sardinia, four locations were selected (Alghero = AHO, Bosa = BOS, Penisola del Sinis = SIN, and Gonnesa = GON for the west and Capo Comino = COM, Cala Gonone = CGO, Arbatax = ARB, and Costa Rei = REI for the east) from 40°34' to 39°15'N. At each location, three areas 100 s of m apart were randomly selected and sampled at a depth of 10 m.

Environmental data. For each location the SST for the years 2010-2019 were obtained by the Group for
High Resolution Sea Surface Temperature (GHRSST) daily, 1 km resolution SST (G1SST) dataset produced by JPL NASA (https:// coast watch. pfeg. noaa. gov/ erddap/ gridd ap/ jplMU RSST41. html) as a proxy of 10 m subtidal temperature 82 . Moreover, 1 Day Composite, 4 km resolution Chlorophyll-a data from NASA's Aqua Spacecraft (https:// coast watch. pfeg. noaa. gov/ erddap/ gridd ap/ erdMH 1chla 1day. html) were extracted for the same years. For the warm season 1st May-31st October (the period of the largest differences between the two coasts), daily SST and Chl-a data were averaged through years (Fig. 7) and the mean, maximum and variance for both variables were calculated (Table 4).

Seagrass data collection.
From 20th June to the 10th July 2020 the density of Posidonia oceanica shoots was estimated using 40 × 40 cm quadrats haphazardly placed within meadows (n = 4) and 20 orthotropic shoots were collected at each area. A total of 480 shoots were collected, transported to the laboratory and stored frozen. Sampling was non-lethal and followed the guidelines approved by the Marine Strategy Framework Directive (EC 2008) for the monitoring program. P. oceanica shoots were deposited as voucher specimens at the University of Sassari Herbarium (SS, collection 2000/, ID number: SS#14159-SS#14166).
The leaf length, leaf width, number of leaves and necrotic leaf portion were measured following 83 to estimate P. oceanica shoot morphometry. Furthermore, the age reconstruction technique based on the cyclic annual variation of the sheath thickness 84 was used to estimate shoot productivity through years: therefore, the number Data analysis. For each P. oceanica variable (shoot density, leaf length, leaf width, number of leaves, necrotic leaf portion, number of scales, rhizome elongation, rhizome width and rhizome biomass) a three-way anova was run to test the effect of 'Coast' (C, west vs east), 'Location' (L, 4 levels) random nested in C, and ' Area' (3 levels) random nested in L. Cochran's test was used to test variance homogeneity. With the aim of finding a relationship between the P. oceanica and the explanatory variables (mean temperature, maximum temperature, temperature variance, mean Chl-a, maximum Chl-a and Chl-a variance, Table 4), we ran separate multiple linear regression models for each P. oceanica response variables. No linear regression was run on leaf length since it is largely affected by herbivore pressure, and it cannot be evaluated unless controlled experiments are performed 85 . Data exploration followed 86 : outliers were inspected with Cleveland dotplots (and removed in four cases) and normality with histograms and Q-Q plots. Rhizome biomass was   www.nature.com/scientificreports/ and variance inflation factors (VIFs) were calculated. Several significant correlations were found, particularly, mean temperature, maximum temperature and temperature variance were correlated to each other, as well as mean Chl-a, maximum Chl-a and Chl-a variance. Thus, only mean temperature and mean Chl-a (the variables with VIFs < 3) and their interaction were considered in the analyses, even though the results obtained for each of them can be extended to all the correlated descriptors. The explanatory variables used in the final model were chosen with a backward selection process 80 . Model validation was run calculating and plotting: (I) standardized residuals against fitted values to assess homogeneity; Figure 5. Summary of the results. Shoot density, morphometry and productivity was measured to examine the performance of the seagrass to the thermal gradient. Warming will trigger denser seagrass meadows with slender shoots (lower rhizome width and biomass) and fewer leaves (scales).  www.nature.com/scientificreports/ (II) histogram of the residuals to verify normality; (III) residuals against each explanatory variable that was used in the model; (IV) residuals against each explanatory variable not used in the model. At the end, the model was assessed for influential observations using the Cook distance function. Correlations between P. oceanica shoot density and all the other plant variables were explored at the scale of area to identify eventual plant traits that might derive from a compensatory performance of the plant to temperature and Chl-a. Thus, following the same methodological approach, another multiple linear regression was run to identify the relationship between shoot density and the other response variables. Since rhizome width was correlated to leaf width and rhizome elongation was correlated to rhizome biomass, the model was run using leaf width, number of leaves and scales and rhizome biomass as predictors. All the analyses were run in R Core Team 87 , using the package MASS 88