Aurélien Boyé
email: aurelien.boye@ifremer.fr

Ronan Becheler

Olivier Gauthier

Vincent Le Garrec

Christian Hily

Marion Maguer

Jacques Grall

Boyé Aurélien

Gauthier Olivier

Becheler Ronan

Le Garrec

Hily Christian

Maguer Marion

Grall Jacques

Drivers and limits of phenotypic responses in vulnerable seagrass populations: Zostera marina in the intertidal

Keywords: aquatic plant ecology, marine phanerogams, stress-gradient hypothesis, cumulative stressors, resilience, phenotypic traits, clonality, intraspecific variability, eelgrass, foundation species

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Seagrasses are marine angiosperms able to form extensive meadows that are critically important for the biodiversity and functioning of marine coastal ecosystems (Duffy, 2006).

They act as foundation species and thereby define many aspects of ecosystems (Ellison, 2019). In particular, they fulfil key roles in carbon sequestration and nutrient cycling, in coastal protection and erosion control, and in the support of commercial fisheries (Barbier et al., 2011). However, the rapid and profound degradation of coastal environments over the past decades [START_REF] Lotze | Depletion, Degradation, and Recovery Potential of Estuaries and Coastal Seas[END_REF] has put seagrasses under high pressure and led to massive shrinkages of their distribution worldwide [START_REF] Orth | Global Crisis for Seagrass Ecosystems[END_REF][START_REF] Waycott | Accelerating loss of seagrasses across the globe threatens coastal ecosystems[END_REF].

Temperate North-East Atlantic meadows have suffered severe losses over the past decades (Airoldi & Beck, 2007). Recently, rates of seagrass losses seem to be decelerating throughout Europe and recovery of fast-growing species have been documented in several locations (de los Santos et al., 2019). Yet, challenges to the recovery of European seagrass meadows still remain. In this region, seagrass diversity is limited, especially in coastal and estuarine areas where only two seagrass species are found: Zostera marina, which is the most common and abundant, and Zostera noltii [START_REF] Hily | The seagrasses of Western Europe[END_REF]Short et al., 2007). Although both species may co-occur, most meadows are monospecific stands [START_REF] Moore | Zostera: Biology, Ecology and Management[END_REF]Short et al., 2007). They are therefore acutely vulnerable to environmental stressors and susceptible to pandemic diseases [START_REF] Waycott | Accelerating loss of seagrasses across the globe threatens coastal ecosystems[END_REF]. In Western Europe, for instance, "wasting disease" outbreaks have led to the die-back of almost 90% of Z. marina populations in the 1930's [START_REF] Godet | Before and after wasting disease in common eelgrass Zostera marina along the French Atlantic coasts: A general overview and first accurate mapping[END_REF]. This remain the major driver of seagrass loss throughout Europe and signs of Z. marina recovery remain limited (de los Santos et al., 2019), impeded by the emergence of new challenges, such as the impact of warming waters on shallow meadows [START_REF] Krause-Jensen | Century-long records reveal shifting challenges to seagrass recovery[END_REF].

In these conditions, resilience of North-East Atlantic seagrass meadows cannot lean on community-based processes and almost entirely relies on the adaptive capacity and response plasticity of the two species present, and in particular of the most common one: Z. marina (Unsworth et al., 2015). Intra-specific phenotypic [START_REF] Maxwell | Phenotypic plasticity promotes persistence following severe events: Physiological and morphological responses of seagrass to flooding[END_REF] and genotypic [START_REF] Hughes | Seagrass genotypic diversity increases disturbance response via complementarity and dominance[END_REF]Reusch et al., 2005)

variability of Zostera

Accepted Article

This article is protected by copyright. All rights reserved species was shown to confer some resilience to the seagrass meadows. However, the capacity and limits of these response mechanisms remain to be fully determined, especially in the face of multiple interacting stressors that are likely to affect seagrass species differently depending on their life history strategies [START_REF] Kilminster | Unravelling complexity in seagrass systems for management: Australia as a microcosm[END_REF][START_REF] O'brien | Seagrass ecosystem trajectory depends on the relative timescales of resistance, recovery and disturbance[END_REF].

A continuum of ecological strategies exists across seagrass taxa, which generally follow a functional trade-off between two extremes [START_REF] Kilminster | Unravelling complexity in seagrass systems for management: Australia as a microcosm[END_REF][START_REF] O'brien | Seagrass ecosystem trajectory depends on the relative timescales of resistance, recovery and disturbance[END_REF].

On the one hand, K-strategists (e.g. Posidonia oceanica) tend to be large and persistent species that mostly invest in clonal reproduction and on resistance traits, such as large energy storage. They often have long resistance and recovery times. On the other hand, r-strategists are smaller, colonizing species (e.g. Ruppia maritima) that tend to rely more on sexual reproduction and recovery traits, such as seed banks, and less on resistance traits. They therefore decline but also recover more quickly. Zostera marina occupies an intermediate position along this resistance-recovery continuum [START_REF] Kilminster | Unravelling complexity in seagrass systems for management: Australia as a microcosm[END_REF][START_REF] O'brien | Seagrass ecosystem trajectory depends on the relative timescales of resistance, recovery and disturbance[END_REF]. Its investment between sexual and clonal reproduction (i.e. rate of clonality) can be extremely variable among and even within meadows (Becheler et al., 2014;[START_REF] Jarvis | Characterization and ecological implication of eelgrass life history strategies near the species' southern limit in the western North Atlantic[END_REF]. As such, Z. marina seems able to tune its strategy across a wide range of the functional trade-off between r-and K-extrema and both annual (opportunistic-like) and perennial strategies have been reported depending on environmental conditions [START_REF] Kim | Annual and perennial life history strategies of Zostera marina populations under different light regimes[END_REF]Van Lent & Verschuure, 1994a).

Zostera marina phenotypic variability represents a challenge from a management perspective to interpret the temporal trajectories of the meadows, identify the drivers behind their changes and adequately evaluate their status and resilience [START_REF] Kilminster | Unravelling complexity in seagrass systems for management: Australia as a microcosm[END_REF]. Variation of Z. marina phenotypic traits, density and reproductive strategy have been documented in relation to a wide range of environmental variables including substrate type, nutrient availability, temperature, salinity, depth, light availability and current regime [START_REF] Jarvis | Characterization and ecological implication of eelgrass life history strategies near the species' southern limit in the western North Atlantic[END_REF][START_REF] Kim | Annual and perennial life history strategies of Zostera marina populations under different light regimes[END_REF][START_REF] Moore | Zostera: Biology, Ecology and Management[END_REF]Van Lent & Verschuure, 1994b). However, these different drivers may exert similar or opposite pressures on seagrass traits (e.g. de Los Santos et al., 2010) and current knowledge of how these multiple responses are integrated and expressed in natural conditions remains sparse (Stockbridge et al., 2020;[START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF]. In particular, the compound effects of

Accepted Article

This article is protected by copyright. All rights reserved climate and hydrodynamic conditions, two major factors through which climate change is expected to affect coastal ecosystems [START_REF] Harley | The impacts of climate change in coastal marine systems[END_REF], with local drivers, such as substrate type and tidal variation, remain to be thoroughly explored in natural intertidal populations.

There is now compelling evidence that the strong physiological constraints associated with intertidal environments may be exacerbated by interactive effects with other environmental factors and anthropogenic pressures [START_REF] Lefcheck | Multiple stressors threaten the imperiled coastal foundation species eelgrass (Zostera marina) in Chesapeake Bay, USA[END_REF]. These interactive effects makes shallow seagrass populations particularly vulnerable in the face of climate changes worldwide and may be a key limit to their recovery in the North-East Atlantic [START_REF] Krause-Jensen | Century-long records reveal shifting challenges to seagrass recovery[END_REF]Unsworth et al., 2015). Hence, important experimental efforts have been devoted to understand the responses of shallow and intertidal Z. marina populations to cumulative stressors (e.g. [START_REF] Moreno-Marín | Additive response to multiple environmental stressors in the seagrass Zostera marina L[END_REF]Zimmerman et al., 2017). Yet, long-term and large-scale monitoring efforts are needed to confirm, generalize and upscale available experimental knowledge in order to guide conservation efforts [START_REF] Lefcheck | Multiple stressors threaten the imperiled coastal foundation species eelgrass (Zostera marina) in Chesapeake Bay, USA[END_REF][START_REF] Potouroglou | The roles of flowering, overwinter survival and sea surface temperature in the long-term population dynamics of Zostera marina around the Isles of Scilly, UK[END_REF][START_REF] Witman | Towards an integration of scale and complexity in marine ecology[END_REF][START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF]. Disentangling the relative and cumulative effects of multiple concomitant factors remains challenging in observational studies (e.g. [START_REF] Kaldy | Factors controlling Zostera marina L. growth in the eastern and western Pacific Ocean: Comparisons between Korea and Oregon, USA[END_REF] and long-term monitoring across large environmental gradients is necessary for this purpose (Verheyen et al., 2017). Currently, most published long-term studies of shallow Z. marina populations phenotypic responses have focused on a limited geographic area or on a single driver [START_REF] Kaldy | Factors controlling Zostera marina L. growth in the eastern and western Pacific Ocean: Comparisons between Korea and Oregon, USA[END_REF][START_REF] Kim | Annual and perennial life history strategies of Zostera marina populations under different light regimes[END_REF][START_REF] Potouroglou | The roles of flowering, overwinter survival and sea surface temperature in the long-term population dynamics of Zostera marina around the Isles of Scilly, UK[END_REF] while those covering large environmental gradients have lacked temporal replication and focused on gradients of sediment and tidal conditions [START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF].

To fill this gap, we study the phenotypic responses of intertidal Z. marina across a region (1) that is a highly diverse environmental mosaic [START_REF] Gallon | Regional and latitudinal patterns of soft-bottom macrobenthic invertebrates along French coasts: Results from the RESOMAR database[END_REF], covering most of the intertidal habitats of Z. marina in the North-East Atlantic (Hily et al., 2003), and (2) was shown to support high genetic and phenotypic diversity for Z. marina (Becheler et al., 2010). We use a regional monitoring of eight intertidal Z. marina beds over five nonconsecutive years along the 500 km coasts of Brittany (France). Monitored variables covered the density, morphology, above-and belowground biomass as well as epiphytic load and were complemented by genotypic data. The aims of this study are to (1)

Accepted Article

This article is protected by copyright. All rights reserved disentangle and quantify the relative importance of a wide range of potential drivers (water temperature, salinity, exposure regime, tidal emersion, sediment and climatic conditions) on Z. marina phenotypic variability, (2) identify potential traits trade-offs constraining Z. marina morphological and structural responses to environmental changes in the intertidal and (3) explore the potential mechanisms underlying the morphological and reproductive strategies of Z. marina to address their implications for conservation purposes.

MATERIAL AND METHODS

Study area and sampling of Zostera marina

Eight intertidal monospecific Zostera marina beds were monitored along the coasts of Brittany (France; Figure 1) for five non-consecutive years, in 2007 and from 2009 to 2012, using data collected in the context of the RÉseau BENThique (REBENT) monitoring programme (http://www.rebent.org). These eight meadows were chosen to encompass the spectrum of environmental settings in which intertidal Z. marina meadows can be found in Brittany: from sheltered bays with fine sediments and turbid waters to exposed areas with coarse sediments and fully marine conditions, through semi-opened habitats [START_REF] Hily | The seagrasses of Western Europe[END_REF][START_REF] Moore | Zostera: Biology, Ecology and Management[END_REF]. Details on the environmental conditions encompassed in this study can be found in Table S1, Figure S7, as well as in Boyé et al.

(2017) for substrate composition.

Sampling was consistently performed for all beds around the spring equinox of each year, between the end of February and the end of April, to limit the influence of the seasonal variability of the descriptors measured [START_REF] Martínez-Crego | Selection of multiple seagrass indicators for environmental biomonitoring[END_REF]. In order to account for within-meadow variability when interpreting differences among meadows (Balestri et al., 2003), three fixed sampling points distributed 200 metres apart were visited within each seagrass bed as shown with the Glénan meadow (Figure 1). At each sampling point, all shoots in two 0.05 m 2 quadrats were collected to measure densities, above and belowground biomass, and describe each shoot's morphology with measures of sheath height, leaves length and width as well as the number of leaves per shoot.

Associated drifting algae biomass were also measured in each quadrat. Additionally, ten shoots were randomly collected at each sampling point to estimate epiphytic load. These

Accepted Article

This article is protected by copyright. All rights reserved shoots were sampled as described above and measured for all variables characterising their morphology and above-and belowground biomass.

Morphological and structural traits measurements

Sheath height was measured from the first node to the separation mark of the leaves.

The length of each leaf was measured from the first node to the apex. The number of broken leaves was counted and expressed as a percentage of the total number of leaves found in each quadrat. One leaf of median length was used to estimate the leaf width for each shoot. Leaves, roots/rhizomes, and brown, red and green algae biomass were estimated as dry weight after 24 hours desiccation at 60°C for each quadrat. Epiphytes were collected using razor blades and their biomass were measured with the same procedure. The ratio between epiphyte and leaf biomass was then computed and averaged across the 10 shoots randomly collected at each sampling point. Total Z. marina leaf and root/rhizome biomass, Z. marina densities and algae biomass were scaled up and expressed per square metre for the two quadrats. To assess the relative investment of Z. marina between its aboveground and belowground parts, a leaves-toroots/rhizomes biomass ratio was calculated within each quadrat based on the total leaf and root/rhizome biomass measured in each. Mean leaf biomass and root/rhizome biomass per shoot were estimated by dividing the total biomass by the shoot densities within each quadrat. For all other variables (sheath height, leaf length and width, number of leaves per shoot, proportion of broken leaves), mean values were calculated for each of the two quadrats separately. Broken leaves were kept in the calculation of mean leaf length to reflect both the physiological and mechanical impacts of the seagrass environment. The two per quadrat values of each sampling point were then averaged to estimate all the variables at the level of the sampling point. All Z. marina morphological and structural traits and their units are summarised in Table1.

Estimation of the rate of clonality

The rate of clonality c, corresponding to the relative frequency of the descendants resulting from clonal reproduction within a population [START_REF] Marshall | Maintenance of genetic variation in apomictic plant populations[END_REF]Stoeckel & Masson, 2014), was inferred in six of the eight meadows in this study (all but Glénan and Sept-Îles). These six meadows were part of genetic surveys in 2009 and 2012 (Becheler

Accepted Article

This article is protected by copyright. All rights reserved et al., 2014All rights reserved et al., , 2010) ) where two quadrats, positioned in the central part of the Z. marina stand and separated by about tens metres, were sampled. Rates of clonality were inferred from the transition of genotypic frequencies between the two sampling years using the Bayesian approach ClonEstiMate (Becheler et al., 2017). This rate varies between 0 (purely sexual population) and 1 (for strict clonality).

Environmental variables

A 25 mm width and 50 mm long sediment core was collected at each sampling point for grain size distribution and organic matter content assessment. Sediments were dried in an oven (24 h at 60 °C) and separated into 15 fractions (<63 μm, 63, 80, 100, 125, 160, 200, 315, 500, 800, 1250, 2000, 3150, 5000 and >10000 μm) whose masses were measured. Fractions were afterwards grouped into gravel (> 2 mm), sand (63 µm to 2 mm) and silt and clay (<63 µm; [START_REF] Fournier | Analyses granulométriques, principes et méthodes[END_REF]. Organic matter content was estimated by mass loss on ignition at 450 °C for 5 hours.

Water temperatures, salinities and currents were obtained from the publicly available database MARC (https://marc.ifremer.fr/en) using the MARS3D model (2.5 km grain, 40 depth levels), previously validated in the study area [START_REF] Lazure | Development of a hydrodynamic model of the Bay of Biscay. Validation of hydrology[END_REF][START_REF] Poppeschi | Unraveling Salinity Extreme Events in Coastal Environments: A Winter Focus on the Bay of Brest[END_REF]. All variables were extracted daily for the years under study at midday near the sediment surface. Linear interpolation on neighbouring grids of the beds coordinates was used to estimate the variables at the site level (no difference for the three points within each site). Meadows being intertidal, and thus close to the land, it was necessary to shift the extraction coordinates away from the coast to avoid model edge effect (Figure S1).

When needed, salinity and temperature data from the literature and environmental monitoring programme such as the Service d'Observation en Milieu Littoral (SOMLIT; Roscanvel meadow) were used to adjust the extraction coordinates. Temperature and salinity HOBO U24-002-C data loggers were also deployed for one month within three meadows to ensure the adequacy between the MARS3D model data and in situ finescale measures. Overall, this model was shown to adequately reproduce coastal hydrology in the study area [START_REF] Lazure | Development of a hydrodynamic model of the Bay of Biscay. Validation of hydrology[END_REF] but while broad-scale variation of hydrology across meadows and years may be well represented, we note that fine-scale temporal variation and extreme events remain underestimated with these data (Figure

Accepted Article

S2; [START_REF] Poppeschi | Unraveling Salinity Extreme Events in Coastal Environments: A Winter Focus on the Bay of Brest[END_REF]. Extracted water temperature, salinity and current time series were characterised by their distribution over the six months preceding the sampling dates of each meadows (in order to match with the response timescale of seagrass structural and demographic traits; [START_REF] O'brien | Seagrass ecosystem trajectory depends on the relative timescales of resistance, recovery and disturbance[END_REF]Roca et al., 2016) with the maximum, minimum and the first, second (median) and third quartile.

The same characterisation was applied to tidal coefficients retrieved from the SHOMAR software (SHOM, 2014) in order to describe the tidal cycle during the six-month period preceding the sampling. The frequencies of spring (tidal coefficient superior to 100) and neap (inferior to 40) tides were also recorded. Tidal variations are common to all sites and reflect solely temporal changes over the whole study area. In addition, estimations of tidal range, exposure time and meadows depth are provided in Table S1 as indication of the overall tidal environment of the meadows. They were not used for the modelling of Zostera marina traits in the redundancy analysis (see Table 1). Indeed, all these meadows are in macrotidal environments (tidal range >4m; [START_REF] Whitfield | Ecosystem and Biotic Classifications of Estuaries and Coasts[END_REF] and exposure time and depth is similar across the meadows. Moreover, these meadows being intertidal, depth is well above reported depth limits for Zostera marina (Duarte et al., 2007;[START_REF] Krause-Jensen | Sea bottom characteristics affect depth limits of eelgrass Zostera marina[END_REF] and was not considered limiting here.

Meteorological conditions were also retrieved from the nearest meteorological stations and integrated over the six months preceding the sampling dates. The meteorological variables comprised the minimum and maximum air temperature and the maximum wind velocity recorded over the six months, the total number of frost days and the cumulative rainfall over this period. Furthermore, the monthly mean of temperature, the monthly mean of the daily range of temperature, and the monthly mean of wind velocity were averaged over the six months. A North-Atlantic Oscillation (NAO) index based on the months preceding sampling (December, January and February) was retrieved from the National Center for Atmospheric Research Staff (http://climatedataguide.ucar.edu, PCbased index accessed in March 2015).

Light is generally not expected to be strongly limiting for intertidal meadows [START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF] and was therefore not surveyed during the monitoring. However, given the hypothesis of potential light limitation advanced by [START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF] after the

Accepted Article

monitoring, photosynthetically active radiation (PAR) data were retrieved a posteriori from EuSeaMap models [START_REF] Populus | EUSeaMap. A European broad-scale seabed habitat map[END_REF]. These PAR data are estimated from MERIS satellite data and available at a 1km resolution as the average PAR (mol photons m -2 day -1 ) per month over the 2005-2009 period, thus covering the beginning of this study period. The values extracted at the meadows coordinates (Figure S3) were all well above the minimum light requirement thresholds reported in the literature for Z. marina (Bertelli & Unsworth, 2018). Thus, they were unlikely to explain any potential light limitation of the meadows under study, and were not included in the redundancy analysis modelling the response of Zostera marina to environmental changes (see Table 1).

Statistical analyses

Balanced datasets are desirable for statistical tests such as multivariate analysis of variance (MANOVA) so that the year 2008, initially monitored in the REBENT programme, was discarded due to missing data. For the same reason, missing organic matter data for two points of Roscanvel in 2007 and one point of Glénan in 2010 as well as root biomass for the three points of Sainte-Marguerite in 2007 were estimated using k-Nearest neighbour imputation. This prevented the removal of entire sites or years while giving neutral weights to these observations. Epiphytes data were missing for the whole year 2011 and thus could not be estimated in this way. Therefore, they were not included as an explanatory variable for the modelling of Zostera marina attributes in the redundancy analysis (see Table 1).

Environmental explanatory variables were divided into five groups to distinguish between the effects of 1) associated flora, 2) substrate conditions, 3) hydrological and climatic conditions, 4) exposure regime and 5) tidal cycles (Table 1). Within each subset, collinear variables were removed using variance inflation factors (VIF) with a threshold of 5 [START_REF] Legendre | Numerical ecology, 3rd English edition[END_REF], leading to the removal of the gravel fraction of the sediment and of some of the variables characterising the distributions of temperature, salinity, currents and tidal coefficient values.

To test the null hypothesis of no variation of Z. marina morphological and structural traits through space and time, a multivariate analysis of variance (MANOVA) was computed by redundancy analysis (RDA; Rao 1964) on the standardized biometric variables and

Accepted Article

This article is protected by copyright. All rights reserved tested by permutations [START_REF] Legendre | Distance-Based Redundancy Analysis: Testing Multispecies Responses in Multifactorial Ecological Experiments[END_REF][START_REF] Legendre | Numerical ecology, 3rd English edition[END_REF]. Sites and years were coded by Helmert contrasts [START_REF] Legendre | Statistical methods for temporal and space-time analysis of community composition data[END_REF] and homogeneity of multivariate dispersions was tested at the α = 0.05 significance level prior to this analysis (Anderson, 2006). This method allows testing and quantifying spatial and temporal variations along with their interaction, which estimates if temporal variations are similar across all sites; or expressed differently, if the spatial patterns are constant through time.

Sampling units in which Z. marina shared similar morphological and structural traits were identified using K-means clustering [START_REF] Hartigan | Algorithm AS 136: A k-means clustering algorithm[END_REF][START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF].

Clustering was performed on the standardized Z. marina data and optimal number of groups between 2 and n/2 (n=120 samples) was identified using the Caliński and Harabasz criterion (1974). Variation of Z. marina's characteristics in space and time were illustrated through the mapping of the groups so defined in an interaction map as proposed by [START_REF] Legendre | Community surveys through space and time: testing the space--time interaction in the absence of replication[END_REF]. Patterns of trait variation were also visualised using Principal Component Analysis (PCA) on the standardized variables. A particular focus was given to the co-variation between Zostera marina densities and its above-and belowground biomass, as it can be indicative of the processes underlying seagrass demography (Cabaço et al., 2013(Cabaço et al., , 2007;;[START_REF] Vieira | Biomass-density data analysis: A comment on Cabaço et al[END_REF][START_REF] Vieira | The biomass-density relationship in seagrasses and its use as an ecological indicator[END_REF].

Redundancy analysis (RDA) was used on the standardized Z. marina traits to relate the variation of Z. marina morphological and structural features to the environmental variables. Stepwise selection based on adjusted coefficient of multiple determination (R 2 adj ) was applied to the explanatory data within each of the five sets of environmental variables previously described (Table 1; Blanchet et al., 2008). Variation partitioning was used to quantify the individual and shared effects of each set of environmental variables [START_REF] Peres-Neto | Variation partitioning of species data matrices: estimation and comparison of fractions[END_REF]. A MANOVA was then performed on the residuals of the RDA model using all sets of selected environmental variables together to evaluate if all the spatial and temporal structures of the Z. marina attributes were explained by our environmental dataset or if we missed important variables constraining Zostera marina development.

Accepted Article

This article is protected by copyright. All rights reserved MARS3D hydrological data were extracting using Python. All other statistical analyses were conducted using R (R Core Team, 2018) and relied on the G2Sd [START_REF] Fournier | G2Sd: a new R package for the statistical analysis of unconsolidated sediments[END_REF], VIM (Templ et al., 2015) and vegan [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2.3-4[END_REF] packages. All the data used in the study, in addition to R scripts to reproduce the analyses and the figures, can be found online at http://doi.org/10.5281/zenodo.5018108.

RESULTS

Spatial and temporal variations of Zostera marina traits

Statistical evidence for a strong spatial and temporal structuring of Z. marina traits was found among the eight meadows and over the five years of the study (Table 2). Sites, years, and their interaction explained 72% of the total variation of Z. marina traits.

Temporal variations, although significant, were of lesser extent (5.9%) compared to the strong signal exclusively attributable to the sites, which explained 43.2% of the total variation of Z. marina traits. The significant interaction between space and time explained 22.9% and indicated that differences among meadows were not constant through the years or, in other words, that temporal changes of the meadows did not concern all sites equally. The site-specific dynamics that underlie this interaction term are illustrated by the spatio-temporal map of the two K-means clusters (Figure 2A).

Separating the sampling units into two groups maximized the Caliński-Harabasz criterion compared to more groups. Saint-Malo, Sept-Iles and Glénan meadows were almost constantly in a first group (Type 1) throughout the study whereas Sainte-Marguerite, Molène and Roscanvel were generally classified in a second group (Type 2). Arcouest and Callot meadows alternated between the two states and were the main contributors to the previously mentioned space-time interaction. Type 1 (Figure 2A; light grey square) was characterised by high shoot densities, low investment in aboveground parts compared to belowground system (both in terms of biomass and length/width of the shoots and leaves) and greater root systems at the square metre scale on average (Figure 2B and S4). Type 2 (Figure 2A; black diamond) comprised sampling units with, on average, lower Z. marina densities and a plant investment more orientated to aboveground parts with higher above to belowground biomass ratios, higher aboveground biomass per shoot, higher leaf and sheath lengths as well as higher leaf

Accepted Article

This article is protected by copyright. All rights reserved width (Figure 2B andS4). The number of leaves and the belowground biomass per shoot showed less variation between the two groups (Figure S4).

Overall, two main configurations emerged from K-means clustering: cases where meadows were constituted by developed shoots in low densities, and cases where meadows displayed higher densities but of smaller shoots. This conspicuous compromise between aboveground biomass per shoot and shoot density was further supported by the absence of any meadow combining high densities with high leaf biomass per shoot (Figure 3). As a result, a negative relationship between the two features of the plant arose at the regional scale. Such relationship was also suggested by the opposite direction theses variables displayed in the PCA (Figure S5). Either one or the other of these variables, but never both, seemed to be favoured according to the sites, which constitutes the main factor discriminating the meadows under study.

Rates of clonality inferred for 11 quadrats distributed across 6 meadows (Table S2) were highly variable, ranging from 0.1 to 0.8. The most clonal meadows were Sainte-Marguerite and Molène (assigned to Type 2) and Arcouest (that switched between the two types over the years) with assessed rates reaching up to 0.7 and 0.8. The lowest rates of clonality (0.1 < c < 0.3) were found in Roscanvel (also assigned to the Type 2)

and Callot (assigned either to Type 1 or 2 depending on the year). Intermediate values were found in Saint-Malo, the only representative of the Type 1 for which we have genotypic data (c values of 0.5 and 0.6).

Effects of environmental constraints on Zostera marina traits

Among the five sets of explanatory variables, all but the one describing the associated flora (F = 1.47, p = 0.14 using 9999 permutations) explained a significant amount of Z. marina trait variation (R 2 adj between 10.4 and 27.8%; p < 0.05). Stepwise selection was then used within the four significant sets to assess the most important forcing variables in explaining the patterns of variation of Z. marina morphological and structural traits.

Variables selected concerned sediment characteristics, seawater salinity and rainfall, atmospheric and seawater temperature, average wind and current conditions, and influence of tidal emersions (Table 1). Overall, the 13 selected variables significantly explained around 40% (R 2 adj ) of the total variation of Z. marina attributes (F = 7.1, p <

Accepted Article

This article is protected by copyright. All rights reserved 0.001 using 9999 permutations). Hydrological and climatic conditions (see Table 1) emerged as the prevailing factors governing Z. marina traits in variation partitioning, explaining 6.2% alone and 24.2% when considering all shared portions of variance (Figure 4). Exposure regime explained 15.1% of Z. marina variations, of which only 0.6% was shared with substrate characteristics. By contrast, exposure regime shared 8.8% of explained variance with hydrological and climatic variables. Substrate characteristics and tidal cycles (the latter only explaining pure temporal variations) equally explained 13.4% of Z. marina variations.

The variation explained by our set of environmental variables (Figure 5) matched to a large extent with the main patterns emphasized by the unconstrained ordination of Z. marina traits (Figure S5). In particular, sites were well discriminated by the RDA model, with the main predicted gradient corresponding to the two groups delineated by K-means clustering i.e. Type 2 (such as Roscanvel and Molène) versus Type 1 (such as Sept-Iles, Saint-Malo and Glénan; Figure 5A &B). Meadows belonging to Type 2 were associated with higher mud content, mean wind velocities, rainfall and median seawater temperatures (Figure 5C). Among Type 1 meadows, the RDA identified different environmental drivers depending on the meadows, suggesting that the same Type 1 trait combination occurs in different environments. The RDA also highlighted subtle trait differences within these Type 1 meadows according to the environment. On the one hand, higher root/rhizome biomass per square metre were found in the Sept-Iles meadow, this Type 1 being best explained by the higher average currents velocities this site is exposed to. Higher percentage of broken leaves, as well as higher epiphytic load (Figure S6) were found in Saint-Malo and Glénan on the other hand. The Type 1 trait combinations found in these two sites were best explained by the high temperature variability (number of days with frost and average daily range of temperature; Figure S7) and by the (mobile) sandy conditions experienced by the plant in these sites.

Characteristics displayed by Z. marina in these two meadows were also correlated with more spring tides (tidal coefficient superior to 100) and bigger tides (higher maximum tidal coefficient) over the six preceding months.

All spatial and temporal signals in the data were not extracted with the set of selected explanatory variables (Table 2, Residuals). Within the 59.8% of variance left unexplained,

Accepted Article

This article is protected by copyright. All rights reserved 46% (significant at α = 0.05 level) remained structured in space and time, representing 27.5% of the total variance (46% of 59.8%). Most of the spatial and temporal structure left unexplained concerned the interaction term, which remained significant within the residuals and represented 18.5% of total variance, while it represented 22.9% in the raw data. Year-to-year variations on the other hand were well accounted for by our statistical model as no pure temporal signal was found significant in the residuals. Likewise, pure spatial signal was relatively well explained by the selected variables. Indeed, residual sites variation only represented 7.7% of the original total variance, which, albeit significant, is 6 times lower that the initially detected spatial signal. Overall, if the RDA failed to accurately model the site-specific variations at the origin of the significant spacetime interaction (e.g. one-off variation occurring in a specific site; see Figure 2A), it explained most of the persistent spatial differences as well as the year-to-year variations of Z. marina phenotypic traits over the region.

DISCUSSION

Extent and drivers of Zostera marina trait variability

This study provides a comprehensive view of the extent and drivers of Zostera marina trait variation across Brittany. This region's 500 km coast is a highly diverse environmental mosaic (Boyé et al., 2017;[START_REF] Gallon | Regional and latitudinal patterns of soft-bottom macrobenthic invertebrates along French coasts: Results from the RESOMAR database[END_REF]. Although its macrotidal shores [START_REF] Whitfield | Ecosystem and Biotic Classifications of Estuaries and Coasts[END_REF] span a limited thermal gradient, generally lying within a suitable range for Z. marina growth (Table S1 & Figure S7; [START_REF] Lee | Effects of irradiance, temperature, and nutrients on growth dynamics of seagrasses: A review[END_REF][START_REF] Nejrup | Effects of salinity and water temperature on the ecological performance of Zostera marina[END_REF][START_REF] Wilson | Climate change projections reveal range shifts of eelgrass Zostera marina in the Northwest Atlantic[END_REF], this region is representative of the diversity of Z. marina intertidal habitats in the North-East Atlantic in terms of substrate, salinity and hydrodynamic conditions [START_REF] Hily | The seagrasses of Western Europe[END_REF][START_REF] Moore | Zostera: Biology, Ecology and Management[END_REF]. Brittany has been described as a hot-spot of genetic and phenotypic diversity for Z. marina (Becheler et al., 2010) and this is in agreement with the extensive variability in morphology, meadow structure and rates of clonality reported here at the regional and within-bed scales.

Overall, most of the persistent spatial differences and year-to-year variations of Z. marina traits were explained by our set of abiotic variables but the low resolution of the available environmental data likely explains the limited ability of our model to capture site-specific one-off variation (Fig S2; [START_REF] Gohin | Satellite and in situ monitoring of chl-a, turbidity, and total suspended matter in coastal waters: Experience of the year 2017 along the french coasts[END_REF][START_REF] Poppeschi | Unraveling Salinity Extreme Events in Coastal Environments: A Winter Focus on the Bay of Brest[END_REF]. In agreement with

Accepted Article

This article is protected by copyright. All rights reserved other North-East Pacific intertidal populations, substrate type and tidal variations were important predictors of Z. marina traits [START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF]. However, trait variation was primarily explained by temperature-related variables. Despite experimental evidence [START_REF] Gustafsson | Algal mats reduce eelgrass (Zostera marina L.) growth in mixed and monospecific meadows[END_REF], algae did not appear to have a significant effect at such broad scales. This hierarchy of drivers therefore supports the tight link existing between broad-scale climate variability and intertidal seagrass populations (Rasheed & Unsworth, 2011;Thom et al., 2014). Nonetheless, it also highlights a non-negligeable role of localscale factors such as local tidal variation that may reinforce climatic constraints [START_REF] Helmuth | Climate change and latitudinal patterns of intertidal thermal stress[END_REF]Thom et al., 2014). These are often neglected in broad-scale seagrass distribution models that mostly focus on temperature and light availability [START_REF] Jayathilake | A modelled global distribution of the seagrass Accepted Article This article is protected by copyright. All rights reserved biome[END_REF][START_REF] Wilson | Climate change projections reveal range shifts of eelgrass Zostera marina in the Northwest Atlantic[END_REF]. As for other foundation species (Firth et al., 2021), better accounting for local drivers may lead to more accurate projections of seagrass distribution changes. In the future, combining finer-scale in situ environmental data with higher-resolution models will also be needed to improve our understanding and predictions of seagrass changes in the face of more frequent and intense extreme events [START_REF] Oliver | Longer and more frequent marine heatwaves over the past century[END_REF][START_REF] Poppeschi | Unraveling Salinity Extreme Events in Coastal Environments: A Winter Focus on the Bay of Brest[END_REF].

Phenotypic strategies of Zostera marina in the intertidal: facilitation-maximizing versus competition-minimizing

Two main trait configurations ("Types") were delineated among the studied meadows.

While most meadows were consistently in the same configuration, some (e.g. Callot and Arcouest) changed over time. The two types are opposite ends of a continuum defined by a trade-off between density and size of shoots (Figure 3). Type 1, found for instance in Glénan, Sept-Iles and Saint-Malo, matches with the typical description of intertidal or shallow Z. marina meadows, namely, a combination of high densities, small shoots and leaves, and a biomass investment orientated towards below-rather than aboveground parts [START_REF] Krause-Jensen | Eelgrass, powerful predictive tool[END_REF][START_REF] Park | Photoacclimatory Responses of Zostera marina in the Intertidal and Subtidal Zones[END_REF]. This trait configuration appears particularly suited to intertidal life, as it is known to enhance moisture retention, thereby reducing thermal stress [START_REF] Maxwell | The fundamental role of ecological feedback mechanisms for the adaptive management of seagrass ecosystems -a review[END_REF], provide shade and photoprotection against high-light stress [START_REF] Park | Photoacclimatory Responses of Zostera marina in the Intertidal and Subtidal Zones[END_REF]Schubert et al., 2015), as well as greater resistance to the erosive forces and sediment instability typical of intertidal areas (Peralta

Accepted Article

This article is protected by copyright. All rights reserved et al., 2008). Type 1 is therefore generally considered as a facilitation-maximizing trait configuration [START_REF] Fischman | Optimizing coastal restoration with the stress gradient hypothesis[END_REF].

At the other end of the continuum, the Type 2 meadows of Sainte-Marguerite, Molène and Roscanvel have low densities of large shoots (in terms of length and aboveground biomass) and a biomass production mainly devoted to photosynthetic parts. Type 2 is more typical of deeper meadows [START_REF] Krause-Jensen | Eelgrass, powerful predictive tool[END_REF] as it maximizes photosynthetic capacity and light harvesting (Ferguson et al., 2016;[START_REF] Ralph | Impact of light limitation on seagrasses[END_REF]. Indeed, larger shoots and density together increase self-shading. Hence, as light becomes more limited with depth, seagrass density tends to decrease whereas shoot length and aboveground biomass per individual exhibit the opposite pattern (Beca-Carretero et al., 2019;Enríquez et al., 2019;[START_REF] Ralph | Impact of light limitation on seagrasses[END_REF]. This process is traditionally referred to as a "self-thinning" response from the plant (Yoda et al., 1963). It results from high intra-specific competition increasing ramet spacing through densitydependent mortality, density-dependent recruitment, or through density-dependent modular growth in clonal plant such as seagrass (Brun et al., 2006;de Kroon, 1993;[START_REF] Johnson | The role of sexual reproduction in the maintenance of established Zostera marina meadows[END_REF][START_REF] Marbà | Scaling of ramet size and spacing in seagrasses: Implications for stand development[END_REF][START_REF] Olesen | Biomass-density patterns in the temperate seagrass Zostera marina[END_REF][START_REF] Yang | Reversal of intraspecific interactions by an ecosystem engineer leads to variable seedling success along a stress gradient[END_REF]. Self-thinning alleviates competition among remaining shoots. Type 2 therefore corresponds to a competition-minimizing trait configuration [START_REF] Fischman | Optimizing coastal restoration with the stress gradient hypothesis[END_REF].

6.3.

Potential mechanisms underlying the two phenotypic strategies: an explanation through the 'stress-gradient hypothesis'

Here, we report a negative relationship between shoot size and density across meadows of similar depth as well as the unexpected presence in the intertidal of a Type 2 (competition-minimizing) strategy, where a Type 1 (facilitation-maximizing) is generally favoured [START_REF] Krause-Jensen | Eelgrass, powerful predictive tool[END_REF][START_REF] Park | Photoacclimatory Responses of Zostera marina in the Intertidal and Subtidal Zones[END_REF]. Based on similar observations in intertidal Zostera marina populations of the Salish sea, [START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF] hypothesized the role of unexpected light limitations in intertidal environments. Nonetheless, such negative relationships have been reported for seagrass species along various environmental gradients, such as physical disturbance, nutrient load, temperature and salinity (Barry et al., 2017;Ferguson et al., 2016;[START_REF] Mcdonald | Seagrass growth, reproductive, and morphological plasticity across environmental gradients over a large spatial scale[END_REF][START_REF] Peralta | Morphometric Accepted Article This article is protected by copyright. All rights reserved variations as acclimation mechanisms in Zostera noltii beds[END_REF]. More generally, the 'stress-gradient hypothesis' (SGH) predicts that facilitative

Accepted Article

This article is protected by copyright. All rights reserved effects prevail over competitive ones under harsh conditions, while competitive effects become more important under milder environmental conditions (Callaway & Walker, 1997). Therefore, Type 1 (facilitation-maximizing) is expected under any "harsh environment" while Type 2 (competition-minimizing) would be expected in any "mild environment", where growth is not constrained by the general environment but by the single most limiting resource [START_REF] Liebig | Die organische Chemie in ihrer Anwendung auf Agricultur und Physiologie[END_REF]. This is often light or nutrients for seagrass [START_REF] Lee | Effects of irradiance, temperature, and nutrients on growth dynamics of seagrasses: A review[END_REF][START_REF] Ralph | Impact of light limitation on seagrasses[END_REF].

In agreement with the SGH, Type 1 (facilitation-maximizing) was associated in our study to two sources of stress: (i) temperature/desiccation (Saint-Malo and Glénan) and (ii) hydrodynamics (Sept-Iles). Those are two of the major factors limiting the distribution of seagrass populations in shallow environments [START_REF] Koch | Beyond light: Physical, geological, and geochemical parameters as possible submersed aquatic vegetation habitat requirements[END_REF]. Type 1 was associated to maximum and minimum temperatures falling outside Z. marina optimal range (Figure S7; [START_REF] Lee | Effects of irradiance, temperature, and nutrients on growth dynamics of seagrasses: A review[END_REF][START_REF] Nejrup | Effects of salinity and water temperature on the ecological performance of Zostera marina[END_REF][START_REF] Wilson | Climate change projections reveal range shifts of eelgrass Zostera marina in the Northwest Atlantic[END_REF]. It was also linked to higher occurrences and bigger spring tides that may locally reinforce climatic constraints [START_REF] Helmuth | Climate change and latitudinal patterns of intertidal thermal stress[END_REF]Thom et al., 2014). These results are therefore in line with previous studies that described the self-facilitative feedbacks promoted by a Type 1 configuration in the face of thermal stress and strong erosive forces [START_REF] Maxwell | The fundamental role of ecological feedback mechanisms for the adaptive management of seagrass ecosystems -a review[END_REF][START_REF] Peralta | Consequences of shoot density and stiffness for ecosystem engineering by benthic macrophytes in flow dominated areas: a hydrodynamic flume study[END_REF].

In contrast, Type 2 (competition-minimizing) was found in meadows exposed to lower hydrodynamic (measured by median current velocities) and lower temperature variability.

Type 2 meadows were exposed to temperatures falling within the optimal range for Z. marina [START_REF] Lee | Effects of irradiance, temperature, and nutrients on growth dynamics of seagrasses: A review[END_REF][START_REF] Nejrup | Effects of salinity and water temperature on the ecological performance of Zostera marina[END_REF][START_REF] Wilson | Climate change projections reveal range shifts of eelgrass Zostera marina in the Northwest Atlantic[END_REF]. Hence, these beds did not seem to be constrained by any other environmental variables than those with a strong link with light availability. While in situ light availability was quantified in this study from broad-scale satellite data, these 4-year average PAR data available at 1km resolution were all above the minimum light requirement thresholds reported in the literature for Z. marina (Fig. S3; Bertelli & Unsworth, 2018). However, in support of [START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF] hypothesis, Type 2 was found in more turbid waters (see Figure 5): it is associated to muddier sediment, higher average wind velocities potentially promoting re-suspension, and higher precipitations, which are strongly linked with riverine inputs in the area (Tréguer et al. 2014). In combination, these environmental factors may

Accepted Article

contribute to reduce light availability, while increasing seagrass light requirements due to muddier conditions (Zabarte-Maeztu et al., 2021). However, such fine-scale variability in in situ water clarity is notoriously hard to retrieve from currently available satellite data [START_REF] Gohin | Satellite and in situ monitoring of chl-a, turbidity, and total suspended matter in coastal waters: Experience of the year 2017 along the french coasts[END_REF]. This interpretation of light as the most limiting factor of these Type 2 meadows is also congruent with their high above-to-belowground biomass ratio (de Los Santos et al., 2010;Ferguson et al., 2016). Indeed, seagrass investment in above-and belowground structures is proportional to resource supply. Belowground-parts are favoured for nutrient acquisition and storage whereas aboveground-parts help maximize photosynthetic capacity under nutrient-sufficient conditions but limiting light-levels [START_REF] Herbert | Phosphorus availability and salinity control productivity and demography of the seagrass Thalassia testudinum in Florida Bay[END_REF].

Overall, the distribution of Type 1 and Type 2 meadows in this study match with the predictions of the SGH. This is in agreement with previous studies supporting (1) the occurrence of self-thinning processes in Z. marina and other seagrass species (Brun et al., 2006;Enríquez et al., 2019;[START_REF] Marbà | Scaling of ramet size and spacing in seagrasses: Implications for stand development[END_REF][START_REF] Olesen | Biomass-density patterns in the temperate seagrass Zostera marina[END_REF][START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF], and ( 2) the congruence of Z. marina dynamics with the SGH [START_REF] Fischman | Optimizing coastal restoration with the stress gradient hypothesis[END_REF][START_REF] Yang | Reversal of intraspecific interactions by an ecosystem engineer leads to variable seedling success along a stress gradient[END_REF]. Predictions of the SGH allow to unify previous observations of seagrass responses across different stress gradients (e.g. Barry et al., 2017;Ferguson et al., 2016;[START_REF] Mcdonald | Seagrass growth, reproductive, and morphological plasticity across environmental gradients over a large spatial scale[END_REF][START_REF] Peralta | Morphometric Accepted Article This article is protected by copyright. All rights reserved variations as acclimation mechanisms in Zostera noltii beds[END_REF]. The SGH can therefore provide a useful simplifying scheme to predict seagrass complex contextdependent feedbacks [START_REF] Maxwell | The fundamental role of ecological feedback mechanisms for the adaptive management of seagrass ecosystems -a review[END_REF][START_REF] Yang | Reversal of intraspecific interactions by an ecosystem engineer leads to variable seedling success along a stress gradient[END_REF] and can help generalize seagrass responses to environmental changes (Figure 6), with several implications in terms of conservation and restoration strategies.

Implications for monitoring, conservation and restoration

As the proposed phenotypic typology is directly related to the SGH, it can be used to adapt restoration efforts to environmental conditions (see [START_REF] Fischman | Optimizing coastal restoration with the stress gradient hypothesis[END_REF]. In terms of monitoring, the convergence toward the Type 1 (facilitation-maximizing) in markedly different environments confirms that seagrass structural traits are sensitive to changes but are not stressor-specific indicators (Roca et al., 2016). Therefore, they might be used to give an integrative portrait of seagrass status but need to be coupled to other indicators, such as remote sensing data (cover and fragmentation), physiological,

Accepted Article

This article is protected by copyright. All rights reserved biochemical and genetic indicators, in order to identify stressor-specific responses across different spatial and temporal scales (Roca et al., 2016). Moreover, recent experimental work suggests that a salinity-stress gradient can result in a Stress-Induced Morphometric Response (SIMR) in several seagrass species with shoot density increasing with stress intensity up to a mortality-threshold (Collier et al., 2014). Here, the concordance of our results with the predictions from the SGH suggests that this SIMR can be generalized to a broad range of other stressors (Figure 6). As a consequence, shoot-proliferation should not be considered, alone, as a consistent indicator of good ecological status for seagrass species.

The design and efficiency of management actions are contingent on the reproductive strategy of the species and the transitory or persistent nature of the meadows [START_REF] Kilminster | Unravelling complexity in seagrass systems for management: Australia as a microcosm[END_REF]. Indeed, rates of clonality (r) deeply influence the resilience and evolution of seagrass [START_REF] O'brien | Seagrass ecosystem trajectory depends on the relative timescales of resistance, recovery and disturbance[END_REF]. The r ranges from a theoretical minimum of 0 for purely sexual population to a maximum of 1 for strictly clonal populations. Here, r ranged from 0.1 to 0.8, confirming the high reproductive variability of Z. marina, even across a 500km coastline (Becheler et al., 2014;[START_REF] Jarvis | Characterization and ecological implication of eelgrass life history strategies near the species' southern limit in the western North Atlantic[END_REF]. Together with the significant spacetime interaction highlighting substantial site-specific phenotypic variation, this high variability of r suggests that conservation efforts should focus on site-specific actions, accounting for the local specificities of the meadows, rather than broad scale measures.

Rates of clonality (r) are influenced by both investment in sexual reproduction and seedling success. Theoretically, high investment in sexual reproduction is expected for Type 1 meadows in stressful environments (Cabaço & Santos, 2010;[START_REF] Kim | Annual and perennial life history strategies of Zostera marina populations under different light regimes[END_REF].

Conversely, Type 2 meadows in stable conditions are expected to follow an Initial Seedling Recruitment (ISR) that should favour high clonality (Eriksson 1993), with low seedling success due to competition [START_REF] Johnson | The role of sexual reproduction in the maintenance of established Zostera marina meadows[END_REF][START_REF] Yang | Reversal of intraspecific interactions by an ecosystem engineer leads to variable seedling success along a stress gradient[END_REF].

Accordingly, the prevalence of persistent clones was highest in some Type 2 meadows (e.g. Molène) but unexpectedly high contribution of sexual reproduction was found in other Type 2 meadows (e.g. Roscanvel). Combination of muddy sediments and calm hydrodynamic conditions (Figure 5) that favour seed germination [START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF] may contribute to the low r of these meadows, but the drivers of these unexpected values warrant further research. The Type 1 Saint-Malo population showed intermediate r and

Accepted Article

This article is protected by copyright. All rights reserved was not the meadow with the highest contribution of sexual reproduction. This intermediate r could partly reflect the low survival of seedlings in this harsh intertidal environment (Valdemarsen et al., 2010;[START_REF] Yang | Relative Impacts of Natural Stressors on Life History Traits Underlying Resilience of Intertidal Eelgrass (Zostera marina L.)[END_REF], despite expectations from the SGH that facilitative effects could enhance seedling success in stressful conditions [START_REF] Yang | Reversal of intraspecific interactions by an ecosystem engineer leads to variable seedling success along a stress gradient[END_REF]. It matches with the mixed-annual reproductive strategy observed in some thermally stressed environment [START_REF] Jarvis | Characterization and ecological implication of eelgrass life history strategies near the species' southern limit in the western North Atlantic[END_REF], where only a fraction of the individuals are renewed through seed recruitment and coexists with persistent clones (Becheler et al., 2014). Overall, while variation of both Z. marina recruitment success [START_REF] Yang | Reversal of intraspecific interactions by an ecosystem engineer leads to variable seedling success along a stress gradient[END_REF] and phenotypic traits (this study) were shown to match predictions from the SGH, observed variation of clonality rates did not match expectations (Figure 6).

If observed r did not appear randomly distributed (extremes in Type 2 meadows, and intermediate in the Type 1 Saint-Malo), more studies are needed to clarify how reproductive modes and phenotypic types are linked across stress-gradients.

Overall, we highlight two distinct phenotypic types with a strong theoretical support from the 'stress-gradient hypothesis' (Callaway & Walker, 1997). While their link with reproductive strategies needs to be elucidated, these two types already provide a useful typology to guide conservation and restoration efforts [START_REF] Fischman | Optimizing coastal restoration with the stress gradient hypothesis[END_REF]. In particular, the trade-off between shoot size and density underlying these two meadow types seems to be one of the key constraints governing Z. marina phenological plasticity in intertidal environments. Indeed, the trait configuration allowing to deal with the strong hydrodynamics and climatic constraints associated with intertidal life (Type 1, facilitationmaximizing) seem incompatible with the responses usually involved in dealing with light limitation (Type 2, competition-minimizing) due to water quality degradation. This can explain the synergistic effects often reported on seagrass when increased temperature co-occur with reduced light availability and increased competition (Stockbridge et al., 2020). Here, we suggest, through predictions of the SGH, that this apparent incompatibility can be generalized to multiple stress gradients. However, how this incompatibility is resolved for different stressors combinations warrants further experimental work. For example, de los Santos et al. (2010) showed that light effects prevailed over hydrodynamic ones in intertidal Z. noltii. More experimental work is therefore needed to understand how and when the SGH can help explain and predict 1 for full description of each set). Only the variables selected by the stepwise procedure were used for this analysis. The redundancy analysis corresponding to this model is presented in Fig. 5. Variation is expressed as a percentage (%) of the total variation of Z. marina traits based on adjusted R 2 . Empty fractions are negative values and may be considered as zeros. Residuals are unexplained variations.

Accepted Article

Figure 5. Redundancy analysis of the standardized Zostera marina traits against the selected environmental variables (same model than in Fig. 4). The response variables and explanatory variables selected by the stepwise procedure are summarised in Table 1 along with their abbreviations and units. Scaling type 2. The first two axes represent 27.85% and 8.11% of the total variation respectively. A. Linear constraints scores for all samples (one point of one site at a given year) with the 95% confidence dispersion ellipses of each site. Within-site dispersions represent the temporal variability during the 5 years of the study, and trait variation among the three points sampled per year, that is explained by the model. B. Response variables i.e. Zostera marina traits. C. Explanatory variables i.e. the environmental variables selected.

Table 1. List of all the variables, with their abbreviations and units, used in this study. Only explanatory variables kept after the removal of collinear variables, identified using Variance Inflation Factor, are presented here. Stepwise selection procedure was performed within each five sets of explanatory variables (associated flora, substrate, hydrological and climatic conditions, exposure regime and tidal cycles) and selected variables are shown in the 5th column.

‡ Epiphytes data were not included as explanatory variables in the modelling of Zostera marina traits due to missing data for an entire year but they were kept as they served interpretation.

Figure 1 .

 1 Figure 1. Localities of the eight Zostera marina meadows (black stars) monitored by the RÉseau BENThique (REBENT) in Brittany in 2007 and from 2009 to 2012. Each site was sampled for Zostera marina phenotypic traits as well as for drifting algae biomass (assessed with two 0.05 m 2 quadrats). The biomass of associated epiphytes was quantified from 10 additional shoots haphazardly collected within the meadow. A core was also sampled for sediment characterisation. This survey was performed within each site at three points located 200 m apart, as shown for the Glénan meadow.

Figure 2 .

 2 Figure 2. A. Spatio-temporal map of the membership of each sampling unit to the groups defined by K-mean clustering performed on the standardized Zostera marina traits. The two groups are represented with different point shapes and intensities of grey. B. Estimated probability density function of key Z. marina traits differentiating these two groups with the vertical dashed-lines representing the mean of the variables for each group. Tick marks across the bottom of each plot indicate observed data points. See Table1for abbreviations and units and FigureS4for the distribution of all the traits measured.

Figure 3 .

 3 Figure 3. Relationship between the mean leaf biomass per shoot and the shoot density of Zostera marina in the 8 sites monitored during the 5 years of the study.

Figure 4 .

 4 Figure 4. Venn diagram illustrating the results of variation partitioning of the standardized traits of the Zostera marina meadows with respect to four sets of explanatory variables: substrate conditions, exposure regime, hydrological and climatic conditions, and tidal cycles (see Table1for full description of each set). Only the variables selected by the stepwise procedure were used for this analysis. The redundancy analysis corresponding to this model is presented in Fig.5. Variation is expressed as a percentage (%) of the total variation of Z. marina traits based on adjusted R 2 . Empty fractions are negative values and may be considered as zeros. Residuals are unexplained variations.

Figure 6 .

 6 Figure 6. Proposed generalisation of the phenotypic variation of Zostera marina across any environmental gradient based on the 'stress-gradient hypothesis' (SGH). The links between phenotypic types and reproductive strategies (bottom panel) warrants further research. Indeed, observed rates of clonality did not match expectations, but sample size was low (indicated between parenthesis), especially for Type 1 meadows (only one meadow with genotypic data).

Accepted Article

This article is protected by copyright. All rights reserved data from the REBENT monitoring programme (http://www.rebent.org) are available in the Quadrige database (http://envlit.ifremer.fr/resultats/base_de_donnees_quadrige) and in the database of the marine observatory of the IUEM (available upon request: https://www-iuem.univ-brest.fr/observatoire). The complete time series are available for each meadow at:  Saint-Malo: https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003702  Arcouest: https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003701  Sept-Iles: https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003700  Callot: https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003703  Sainte-Marguerite: https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003704  Molène: https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003705  Roscanvel : https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003706  Glénan : https://wwz.ifremer.fr/surval/Donnees/Graphes-30140#/sensor/60003707 11. REFERENCES Airoldi, L., & Beck, M. W. (2007). Loss, status and trends for coastal marine habitats of Europe. Biology: An Annual Review, 45, 357-405. doi: 10.1201/9781420050943.ch7 Anderson, M. J. (2006). Distance-based tests for homogeneity of multivariate dispersions.

Oceanography and Marine

Biometrics, 62 (1), 245-253. doi: 10.1111/j.1541-0420.2005.00440.x Balestri, E., Cinelli, F., & Lardicci, C. (2003). Spatial variation in Posidonia oceanica structural, morphological and dynamic features in a northwestern Mediterranean coastal area: A multi-

Accepted Article

This article is protected by copyright. All rights reserved scale analysis. Marine Ecology Progress Series, 250, 51-60. doi: 10.3354/meps250051 Barbier, E. B., Hacker, S. D., Kennedy, C., Koch, E. W., Stier, A. C., & Silliman, B. R. (2011). The value of estuarine and coastal ecosystem services. Ecological Monographs, 81 (2), 169-193. doi: 10.1890/10-1510.1 Barry, S. C., Jacoby, C. A., & Frazer, T. K. (2017). Environmental influences on growth and morphology of Thalassia testudinum. Marine Ecology Progress Series, 570, 57-70.

Beca-Carretero, P., Stanschewski, C. S., Julia-Miralles, M., Sanchez-Gallego, A., & Stengel, D. B. (2019). Temporal and depth-associated changes in the structure, morphometry and production of near-pristine Zostera marina meadows in western Ireland. Aquatic Botany, 155, 5-17. doi: 10.1016/j.aquabot.2019.02.003 Becheler, R., Benkara, E., Moalic, Y., Hily, C., & Arnaud-Haond, S. (2014). Scaling of processes shaping the clonal dynamics and genetic mosaic of seagrasses through temporal genetic monitoring. Heredity, 112 (2), 114-121. doi: 10.1038/hdy.2013.82 Becheler, R., Diekmann, O., Hily, C., Moalic, Y., & Arnaud-Haond, S. (2010). The concept of population in clonal organisms: mosaics of temporally colonized patches are forming highly diverse meadows of Zostera marina in Brittany. Molecular Ecology, 19, 2394-2407. doi: 10.1111/j.1365-294X.2010.04649.x Becheler, R., Masson, J. P., Arnaud-Haond, S., Halkett, F., Mariette, S., Guillemin, M. L., Valero, M., Destombe, C., & Stoeckel, S. (2017). ClonEstiMate, a Bayesian method for quantifying rates of clonality of populations genotyped at two-time steps. Molecular Ecology Resources, 17(6), e251-e267. doi: 10.1111Resources, 17(6), e251-e267. doi: 10. /1755-0998.12698 -0998.12698 Bertelli, C. M., & Unsworth, R. K. F. ( 2018). Light stress responses by the eelgrass, Zostera marina (L). Frontiers in Environmental Science, 6, 1-13. doi: 10.3389/fenvs.2018.00039 Blanchet, F. G., Legendre, P., & Borcard, D. (2008). Forward selection of explanatory variables.

Ecology, 89(9), 2623-2632. doi: 10.1890/07-0986.1 Boyé, A., Legendre, P., Grall, J., & Gauthier, O. (2017). Constancy despite variability: Local and

Accepted Article

This article is protected by copyright. All rights reserved regional macrofaunal diversity in intertidal seagrass beds. Journal of Sea Research, 130, 107-122. doi: 10.1016Research, 130, 107-122. doi: 10. /j.seares.2017.06.004 .06.004 Brun, F. G., Vergara, J. J., Peralta, G., García-Sánchez, M., Hernández, I., & Pérez-Lloréns, J. L. (2006). Clonal building, simple growth rules and phylloclimate as key steps to develop functional-structural seagrass models. Marine Ecology Progress Series, 323,[133][134][135][136][137][138][139][140][141][142][143][144][145][146][147][148]/meps323133 Cabaço, S., Apostolaki, E. T., García-Marín, P., Gruber, R., Hernández, I., Martínez-Crego, B., Mascaró, O., Pérez, M., Prathep, A., Robinson, C., Romero, J., Schmidt, A. L., Short, F. T., van Tussenbroek, B. I., & Santos, R. (2013). Effects of nutrient enrichment on seagrass population dynamics: evidence and synthesis from the biomass-density relationships. Journal of Ecology, 101(6), 1552-1562. doi: 10.1111/1365-2745.12134 Cabaço, S., Machás, R., & Santos, R. (2007). Biomass-density relationships of the seagrass Zostera noltii: A tool for monitoring anthropogenic nutrient disturbance. Estuarine, Coastal and Shelf Science, 74, 557-564. doi: 10.1016/j.ecss.2007.05.029 Cabaço, S., & Santos, R. (2010). Reproduction of the eelgrass Zostera marina at the species southern distributional limit in the Eastern Atlantic. Marine Ecology, 31(2), 300-308. doi: 10.1111/j.1439-0485.2009.00331.x Caliński, T., & Harabasz, J. (1974). A dendrite method for cluster analysis. Communications in Statistics-Theory and Methods, 3(1), 1-27. Callaway, R. M., & Walker, L. R. (1997). Competition and facilitation: a synthetic approach to interactions in plant communities. Ecology, 78, 1958-1965. doi: 10.1890/0012-9658(1997)078[1958:CAFASA]2.0.CO;2 Collier, C. J., Villacorta-Rath, C., van Dijk, K., Takahashi, M., & Waycott, M. (2014). Seagrass Proliferation Precedes Mortality during Hypo-Salinity Events: A Stress-Induced Morphometric Response. PLoS ONE, 9(4), e94014. doi: 10.1371/journal.pone.0094014 de Kroon, H. (1993). Competition between Shoots in Stands of Clonal Plants. Plant Species

Accepted Article

This article is protected by copyright. All rights reserved Biology, 8 (2-3), 85-94. doi: 10.1111/j.1442-1984.1993.tb00060.x de Los Santos, C. B., Brun, F. G., Bouma, T. J., Vergara, J. J., & Pérez-Lloréns, J. L. (2010).

Acclimation of seagrass Zostera noltii to co-occurring hydrodynamic and light stresses. Series, 398,[127][128][129][130][131][132][133][134][135].3354/meps08343 de los Santos, C. B., Krause-Jensen, D., Alcoverro, T., Marbà, N., Duarte, C. M., van Katwijk, M. M., Pérez, M., Romero, J., Sánchez-Lizaso, J. L., Roca, G., Jankowska, E., Pérez-Lloréns, J. L., Fournier, J., Montefalcone, M., Pergent, G., Ruiz, J. M., Cabaço, S., … Santos, R. ( 2019). Recent trend reversal for declining European seagrass meadows. Nature Communications, 10(1). doi: 10.1038/s41467-019-11340-4 Duarte, C. M., Marbà, N., Krause-Jensen, D., & Sánchez-Camacho, M. (2007). Testing the predictive power of seagrass depth limit models. Estuaries and Coasts, 30(4), 652-656. doi: 10.1007/BF02841962 Duffy, J. E. (2006). Biodiversity and the functioning of seagrass ecosystems. Marine Ecology Progress Series, 311, 233-250. doi: 10.3354/meps311233 Ellison, A. M. (2019). Foundation Species, Non-trophic Interactions, and the Value of Being Common. IScience, 13, 254-268. doi: 10.1016/j.isci.2019.02.020 Enríquez, S., Olivé, I., Cayabyab, N., & Hedley, J. D. (2019). Structural complexity governs seagrass acclimatization to depth with relevant consequences for meadow production, macrophyte diversity and habitat carbon storage capacity. Scientific Reports, 9(1), 1-14. doi: 10.1038/s41598-019-51248-z Eriksson, O. (1993). Dynamics of genets in clonal plants. Trends in Ecology and Evolution, 8(9), 313-316. doi: 10.1016/0169-5347(93)90237-J Ferguson, A. J. P., Gruber, R. K., Orr, M., & Scanes, P. (2016). Morphological plasticity in Zostera muelleri across light, sediment, and nutrient gradients in Australian temperate coastal lakes. Series, 556, 91-104. doi: 10.3354/meps11830 Firth, L. B., Harris, D., Blaze, J. A., Marzloff, M. P., Boyé, A., Miller, P. I., Curd, A., Vasquez, M.,

Marine Ecology Progress

Marine Ecology Progress

Accepted Article

This article is protected by copyright. All rights reserved Nunn, J. D., O'Connor, N. E., Power, A. M., Mieszkowska, N., O'Riordan, R. M., Burrows, M.

Accepted Article

This article is protected by copyright. All rights reserved Rao, R. (1964). The Use and Interpretation of Principal Component Analysis in Applied Research.

The Indian Journal of Statistics, 26(4), 329-358. Rasheed, M. A., & Unsworth, R. K. F. (2011). Long-term climate-associated dynamics of a tropical seagrass meadow: implications for the future. Marine Ecology Progress Series, 422, 93-103. doi: 10.3354/meps08925 Reusch, T. B. H., Ehlers, A., Hämmerli, A., & Worm, B. (2005). Ecosystem recovery after climatic extremes enhanced by genotypic diversity. Proceedings of the National Academy of Sciences, 102(8), 2826-2831. doi: 10.1073/pnas.0500008102 Roca, G., Alcoverro, T., Krause-Jensen, D., Balsby, T. J. S., van Katwijk, M. M., Marba, N., Santos, R., Arthur, R., Mascara, O., Fernandez-Torquemada, Y., Perez, M., Duarte, C. M., & Romero, J. (2016). Response of seagrass indicators to shifts in environmental stressors: A global review and management synthesis. Ecological Indicators, 63, 310-323. doi: 10.1016/j.ecolind.2015.12.007 Schubert, N., Colombo-Pallota, M. F., & Enríquez, S. (2015). Leaf and canopy scale characterization of the photoprotective response to high-light stress of the seagrass Thalassia testudinum. Limnology and Oceanography. doi: 10.1002/lno.10024 SHOM. (2014). SHOMAR -Édition 2014-2015, disponible sur http://www.shom.fr/lesproduits/bases-de-donnees-numeriques/maree-et-courant/shomar/. Short, F., Carruthers, T., Dennison, W., & Waycott, M. (2007). Global seagrass distribution and diversity: A bioregional model. Journal of Experimental Marine Biology and Ecology, 350(1-2), 3-20. doi: 10.1016/j.jembe.2007.06.012 Stockbridge, J., Jones, A. R., & Gillanders, B. M. (2020). A meta-analysis of multiple stressors on seagrasses in the context of marine spatial cumulative impacts assessment. Scientific Reports, 10(1), 1-11. doi: 10.1038/s41598-020-68801-w Stoeckel, S., & Masson, J. P. (2014). The exact distributions of FIS under partial asexuality in small finite populations with mutation. PLoS ONE, 9(1). doi: 10.1371/journal.pone.0085228

Accepted Article

This article is protected by copyright. All rights reserved Templ, M., Alfons, A., Kowarik, A., & Prantner, B. (2015). VIM: Visualization and Imputation of Missing Values. R package version 4.4.1. http://CRAN.R-project.org/package=VIM. Thom, R., Southard, S., & Borde, A. (2014). Climate-linked Mechanisms Driving Spatial and Temporal Variation in Eelgrass (Zostera marina L.) Growth and Assemblage Structure in Pacific Northwest Estuaries, U.S.A. Journal of Coastal Research, 68(sp1), 1-11. doi: 10.2112/si68-001.1 Tréguer, P., Goberville, E., Barrier, N., L'Helguen, S., Morin, P., Bozec, Y., Rimmelin-Maury, P., Czamanski, M., Grossteffan, E., Cariou, T., Répécaud, M., & Quéméner, L. (2014). Large and local-scale influences on physical and chemical characteristics of coastal waters of Western Europe during winter. Journal of Marine Systems, 139(November), 79-90. doi: 10.1016/j.jmarsys.2014.05.019 Unsworth, R. K. F., Collier, C. J., Waycott, M., Mckenzie, L. J., & Cullen-Unsworth, L. C. (2015). A framework for the resilience of seagrass ecosystems. Marine Pollution Bulletin, 100(1), 34-46. doi: 10.1016/j.marpolbul.2015.08.016 Valdemarsen, T., Canal-Verges, P., Kristensen, E., Holmer, M., Kristiansen, M. D., & Flindt, M. R. (2010). Vulnerability of Zostera marina seedlings to physical stress. Marine Ecology Progress Series, 418, 119-130. doi: 10.3354/meps08828 Van Lent, F., & Verschuure, J. M. (1994a). Intraspecific variability of Zostera marina L. (eelgrass) in the estuaries and lagoons of the southwestern Netherlands. I. Population dynamics. Aquatic Botany, 48(1), 31-58. doi: 10.1016/0304-3770(94)90072-8 Van Lent, F., & Verschuure, J. M. (1994b). Intraspecific variability of Zostera marina L. (eelgrass) in the estuaries and lagoons of the southwestern Netherlands. II. Relation with environmental factors. Aquatic Botany, 48(1), 59-75. doi: 10.1016/0304-3770(94)90072-8 Verheyen, K., De Frenne, P., Baeten, L., Waller, D. M., Hédl, R., Perring, M. P., Blondeel, H., Brunet, J., Chudomelová, M., Decocq, G., De Lombaerde, E., Depauw, L., Dirnböck, T., Durak, T., Eriksson, O., Gilliam, F. S., Heinken, T., … Bernhardt-Römermann, M. (2017). Combining biodiversity resurveys across regions to advance global change research. BioScience, 67(1),

Accepted Article

Accepted Article

This article is protected by copyright. All rights reserved Yoda, K., Kira, T., Ogawa, H., & Hozumi, K. (1963). Self-thinning in overcrowded pure stands under cultivated and natural conditions (Intraspecific competition among higher plants XI).

Journal of the

Institute of Polythecnics, Osaka City University, Series D, 14, 107-129. Zabarte-Maeztu, I., Matheson, F. E., Manley-Harris, M., Davies-Colley, R. J., & Hawes, I. (2021). Interaction of substrate muddiness and low irradiance on seagrass: A mesocosm study of Zostera muelleri. Aquatic Botany, 175. doi: 10.1016/j.aquabot.2021.103435 Zimmerman, R., Hill, V., Jinuntuya, M., Celebi, B., Ruble, D., Smith, M., Cedeno, T., & Swingle, W. (2017). Experimental impacts of climate warming and ocean carbonation on eelgrass Zostera marina. Marine Ecology Progress Series, 566, 1-15. doi: 10.3354/meps12051 Accepted Article

synergistic impacts on seagrass (Stockbridge et al., 2020). Moreover, occurrence of a given phenotypic type is likely to affect (positively or negatively) the resilience of shallow seagrass populations in the face of extreme events (e.g. [START_REF] Maxwell | Phenotypic plasticity promotes persistence following severe events: Physiological and morphological responses of seagrass to flooding[END_REF]. As these are becoming more frequent and intense [START_REF] Oliver | Longer and more frequent marine heatwaves over the past century[END_REF], future seagrass monitoring programmes will require finer-scale environmental data in order to better apprehend the effect of these drivers of seagrass changes. Increasing hydrodynamics and climatic constraints [START_REF] Harley | The impacts of climate change in coastal marine systems[END_REF] are already leading to strong negative interactive effects in conjugation with water quality degradation [START_REF] Krause-Jensen | Century-long records reveal shifting challenges to seagrass recovery[END_REF][START_REF] Lefcheck | Multiple stressors threaten the imperiled coastal foundation species eelgrass (Zostera marina) in Chesapeake Bay, USA[END_REF]. Water quality policies are therefore key pressure points that need to be used to ensure the resilience of seagrass meadows (Unsworth et al., 2015), especially for intertidal and shallow populations as our results suggest they face high risks of strong cumulative effects under global change stressors.

ACKNOWLEDGMENTS

We thank Camille Mazoyer and Mickaël Vasquez for their help with the retrieval of the MARS3D hydrological and EuSeaMap PAR data. This study was supported by the REBENT programme coordinated by Sandrine Derrien (MNHN) and its funding partners (Agence de l'eau Loire-Bretagne, Région Bretagne, DREAL Bretagne); the "Laboratoire d'Excellence" LabexMER (ANR-10-LABX-19) co-funded by grants from the French government under the "Investissements d'Avenir" programme and the Regional Council of the Région Bretagne; and a CNRS-UBO grant to O.G.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest. All authors gave final approval for publication.

AUTHORS' CONTRIBUTIONS

DATA AVAILABILITY

All the data used in the study, in addition to R scripts to reproduce the analyses and the figures, can be found online at http://doi.org/10.5281/zenodo.5018108. The complete I] being the interaction terms between sites and years.

[I|S+T] represents the pure effect of interaction, once the effect of spatial and temporal contrasts was removed, and likewise for the pure spatial and temporal fractions. Spatial and temporal structures were first quantified and tested on the standardized Zostera marina traits (Raw data). A second analysis was then performed on the residuals of the RDA modelling Z. marina traits with the selected environmental variables (Fig. 4 and 5) to test for remaining spatial and temporal signals once the explanatory variables effects were removed. Unadjusted R 2 are presented as percentage of total variance and p-values were obtained after 9999 permutations with *** < 0.001,** < 0.01 and * < 0.05. Residual unadjusted R 2 were also expressed as percentage of total variance by multiplying the R 2 with the unexplained variance of the RDA model (1-R 2 adj = 59.8%)