Trait Gradients Predict Global Warming Changes in Seagrass Meadows

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 quantied 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 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.


Introduction
Global warming is expected to have profound consequences on biodiversity and functioning of major systems on Earth [1,2]. The impact of temperature increase has been measured over the past two decades [3][4][5][6], but understanding how this physical forcing affects ecosystems is unclear, particularly in the sea [7][8][9]. This, however, is critical for predicting the consequences of global warming and identify mitigation and restoration actions.
Much of experimental temperate reef ecology is focused on elucidating how temperature increases will impact the physiology, tness and distribution of organisms. Two main approaches are being employed to examine warming effects: (I) experiments with arti cial heating such as mesocosms [10][11][12][13] and (II) monitoring the response of organisms to temporal or spatial variation in temperature, across years [14][15][16] or latitude [17][18][19][20]. Each of these approaches has advantages and drawbacks. Manipulative experiments may examine responses to temperature or patterns not yet under natural conditions, such as intense, long lasting heat waves [21][22][23]. Experiments are typically done at small spatiotemporal scales and often ignore covarying abiotic conditions including light availability [24], UV irradiation and acidi cation [25,26], or biotic effects such as predation [27,28]. Conversely, comparing populations across sites with varying temperatures, such as latitudinal gradients, can provide information about the role of warming on the structure and function of future ecosystems, but it is often di cult to disentangle temperature from other covarying effects, such as photoperiod, light quality and quantity [29]. Moreover, marine sea surface temperature (SST) is commonly linked to chlorophyll-a (Chl-a), with high-temperature locations having low-nutrient availability and Chl-a [30,31] and high light attenuation [32,33]. Problems between laboratory and eld results are not surprising, since temperature, nutrients and irradiance effects may be cumulative or antagonistic depending on the species and system. Therefore, uncertainties with warming effects on marine biota are also indirectly due to co-variation between SST and Chl-a. While there are latitudes where these patterns are predictable, regional anomalies are also found especially where upwelling occurs [34]. Nevertheless, SST increase does not necessarily imply decreasing Chl-a, suggesting that complex processes, such as advection, de ne sea water conditions [34]. Further variability comes from natural variation in species physiological, morphological and life-history attributes among populations, as there is evidence of adaptation to spatial temperature gradients in many organisms and at different scales [35][36][37][38][39]. Species phenotypic gradients presumably can re ect patterns of genetic differentiation and local adaptation, making additional data potentially necessary to estimate how much of observed phenotypic differences are due to plastic responses versus adaptive differentiation between populations.
Understanding future warming effects on foundation species is pivotal to predict their distribution and physical structure [40]. Seagrasses are valuable providers of coastal ecosystem services including, carbon sinks, nursery grounds, habitat, nutrient cycling, sediment stabilization, trophic transfer to adjacent habitats [41][42][43] and protection from erosion [44,45]. Posidonia oceanica (L.) Delile is a slow-growing seagrass, endemic to the Mediterranean, experiencing widespread decline due to multiple local anthropogenic stressors [46]. The abrupt decline experienced by P. oceanica from recent heatwaves [47], 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 [48] and molecular mechanisms [49][50][51] to tolerate simulated heatwave exposures, whereas phenological response to warming likely involves higher owering [52] and denser meadows [53].
We assumed that plant functional traits vary along environmental gradients and potentially predict responses to environmental change. 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. Western locations are generally cooler than eastern sites, despite similar latitude. This gradient allowed examination of the response of P. oceanica to a wide range of SST, with minimum interference of photoperiod. Chl-a, a proxy of nutrient availability and light irradiance, was a further driver of seagrass structure. P. oceanica is currently in the EU Marine Strategy Framework Directive monitoring protocols [54].

Study locations and design
This study was done on the western and eastern coasts of Sardinia (Italy, western Mediterranean Sea, Fig.1) where differences in water conditions are evident. The western coastline receives Atlantic waters directly through the Western Mid-Mediterranean Current and is also in uenced by coastal upwellings [74]. In contrast, the eastern coast is affected by the warm Algerian Current [75].
Seagrass meadows far from anthropogenic disturbances were sampled in eight different locations (Fig. 1), 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 100s of m apart were randomly selected and sampled at a depth of 10m.

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://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURSST41.html) as a proxy of 10 m subtidal temperature [76]. Moreover, 1 Day Composite, 4 km resolution Chlorophyll-a data from NASA's Aqua Spacecraft (https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMH1chla1day.html) were extracted for the same years. For the warm season 1 st May-31 st October (the period of the largest differences between the two coasts), daily SST and Chl-a data were averaged through years (Fig. 2) and the mean, maximum and variance for both variables were calculated (Table 1).

Seagrass data collection
From 20 th June to the 10 th 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.
The leaf length, leaf width, number of leaves and necrotic leaf portion were measured following [77] to estimate P. oceanica shoot morphometry. Furthermore, the age reconstruction technique based on the cyclic annual variation of the sheath thickness [78] was used to estimate shoot productivity through years: therefore, the number of leaves (by counting the scales), rhizome elongation, rhizome width and biomass per year were measured on each shoot (after drying rhizomes for 48h at 60°C).

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 nding 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 1), we ran separate multiple linear regression models for each 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 [79]. Data exploration followed [80]: outliers were inspected with Cleveland dotplots (and removed in four cases) and normality with histograms and Q-Q plots. Rhizome biomass was square root transformed. Collinearity between continuous explanatory variables was inspected with pairplots, and variance in ation factors (VIFs) were calculated. Several signi cant 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, 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 nal model were chosen with a backward selection process [80]. Model validation was run calculating and plotting: I) standardized residuals against tted values to assess homogeneity; 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 in uential 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 [81], using the package MASS [82].

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

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 3, Figs 5 and 6). Speci cally: 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. 5) and II) reduced shoot density, number of scales and rhizome width were correlated with increased Chl-a ( Fig. 6). 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 (Table 3, Figs 5 and 6).
Finally, the regression model indicated that shoot density was negatively related to the number of leaves and leaf width ( Table  4).

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 de ning its conservation status (Marine Strategy Framework Directive) assuming that higher densities re ect lower human in uence and better marine water conditions. However, the density classes distinguished by previous authors (reviewed by [55]), ignore natural environmental variation. Since our data were collected far from 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 revealing 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 nding is corroborated by long-term correlative data revealing that shoot density is a plant trait that varies with thermal environment [53], providing evidence that the plant would rearrange (increasing the number of modules) the meadows structure with warming (Fig. 7). The fact that Chl-a is inversely related to the meadow density will sharpen this pattern, as this in uence 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 [56][57][58][59]. Nevertheless, predicting meadow structure based on the trait-spatial gradient association assumes that the seagrass traits could change proportionally to climate change [53], although the species may respond to ner-scale changes in environmental variables that cannot be predicted using averages [60,61].
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 [62-64], supporting the hypothesis that light extinction is pivotal.
[65] observed that seagrasses growing in low light reduce shoot density and aboveground biomass as an acclimation response to reduce self-shading 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 self-organized to self-shading [66][67][68]. Reducing the size of ramets to attenuate intraspeci c competition is a common pattern in clonal plants [69,70].
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 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. Quanti cation of seagrass services, however, have never been provided on a structure-speci c basis [71,72] and we believe this might become a relevant issue. Indeed, in a future warmer Mediterranean Sea 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.
This study shed light on how seagrass systems could respond to climate change, although the extent the phenotypic gradients 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 [73]. 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.  Figure 1 Study locations and areas along the Sardinian coasts. Left-hand map shows locations on the west (in blue) and east (in red) coasts: AHO=Alghero, BOS=Bosa, SIN=Penisola del Sinis, GON=Gonnesa, COM=Capo Comino, CGO=Cala Gonone, ARB=Arbatax, REI=Costa Rei. Right-hand inset maps show location of each study area within each location. Map produced with QGIS 3.16 software (https://qgis.org/en/site/forusers/download.html). Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries.
This map has been provided by the authors.

Figure 2
Mean variation from 1st May to 31st October in a) temperature (°C), in blue the western locations and in red the eastern, and b) Chl-a (mg/m3), in green the western locations and in yellow the eastern). The y axis of the latter plot is log2 scale.  Posidonia oceanica. Mean (+SE) morphometry (left) and productivity (right) variables. Morphometry: # of leaves/shoot, leaf width (cm), length (cm), and necrotic leaf portion (%). Productivity: # of scales/shoot*yr, rhizome elongation (cm/yr), rhizome width (cm/yr) and rhizome biomass (g/yr) across locations, in blue the western and in red the eastern. For each location data of the three areas are shown (n=20).

Figure 6
Plots from the multiple regression model of Posidonia oceanica response variables vs. mean Chl-a (mg/m3).

Figure 7
Summary of the results.