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 influenced 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 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. 2) and the mean, maximum and variance for both variables were calculated (Table 1).
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.
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 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 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 pair-plots, 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, 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; 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 [81], using the package MASS [82].