Cormy Gaetan

Le Mao Patrick

Desroy Nicolas

Characterizing spatio-temporal changes in benthic communities: Taxonomic and functional trajectories of intertidal assemblages in the bay of Saint-Brieuc (English Channel)

Keywords: community trajectory analysis, ecological trajectories, benthic macrofauna, stability, variability, marine protected area

 and 2019. The taxonomic and functional trajectories of intertidal benthic assemblages were studied at local and assemblage scales. Whatever the perennial contribution of dominant prevalent species at the site scale, some species were highly variable over time and space, influencing the structuration of assemblages. Taxonomic changes were not systematically associated with significant changes in the functional properties at local and community scale. The consideration of local changes and variations at the species level explained the differences in composition observed at the scale of assemblages. Observed changes mainly reflected random population dynamics of structuring prevalent species in habitat under strong natural pressures, rather than specific changes of the community. Few species reveal the long-term effect of green tides in some assemblages. Overall, this study confirms the complexity to disentangle factors of changes in the absence of a monitoring strategy specifically devoted the measurement of environmental factors and anthropogenic pressures. The study also highlights the need for the development of benthic monitoring strategies at multiple spatial and temporal scales.

Introduction

Costal ecosystems face severe natural and anthropogenic pressures that affect, sometimes strongly, both the structure and functioning of benthic communities [START_REF] Harris | Anthropogenic Threats to Benthic Habitats[END_REF][START_REF] Cimon | Impact of multiple disturbances and stress on the temporal trajectories and resilience of benthic intertidal communities[END_REF]. Although assessment of stability/variability of benthic communities, often in relation with management considerations, has received a large attention by scientists, it remains a central point of ecological research and conservation [START_REF] Ellis | Broad-scale disturbance of intertidal and shallow sublittoral soft-sediment habitats; effects on the benthic macrofauna[END_REF]. Studying stability/variability in ecological systems requires the consideration of disturbances defined as "any discrete event in time that disrupt ecosystems, community or populations structure" [START_REF] Picket | The Ecology of Natural Disturbance and Patch Dynamics[END_REF]. Ecological theory suggests that small scale disturbance events at intermediate frequencies may have positive effects on the system at larger scales whereas disturbance on larger scales are much more likely to reduce heterogeneity in communities and induce slower recovery rates [START_REF] Bazzaz | Characteristics of Populations in Relation to Disturbance in Natural and Man-Modified Ecosystems[END_REF][START_REF] Caswell | Disturbance, interspecific interaction and diversity in metapopulations[END_REF][START_REF] Connell | Diversity in tropical rainforests and coral reefs[END_REF][START_REF] Ellis | Broad-scale disturbance of intertidal and shallow sublittoral soft-sediment habitats; effects on the benthic macrofauna[END_REF]. Understanding the effects of environmental changes on natural communities requires to distinguish natural and anthropogenic variability [START_REF] Osenberg | Detecting Ecological Impacts Caused by Human Activities[END_REF][START_REF] Stewart-Oaten | Goals in Environmental Monitoring[END_REF]. However, the lack of reference, non-impacted sites in modern ecosystems, makes this task a real challenge [START_REF] Boldina | Effect of long-term mechanical perturbation on intertidal soft-bottom meiofaunal community spatial structure[END_REF]. Temporal surveys assessing the evolution of benthic communities are therefore essential to establish baselines and detect ecological changes [START_REF] Desroy | Macrobenthic resources of the shallow soft-bottom sediments in the eastern English Channel and southern North Sea[END_REF][START_REF] Hardman-Mountford | Diagnostic monitoring of a changing environment: An alternative UK perspective[END_REF][START_REF] Schückel | Macrofauna communities of tidal channels in Jade Bay (German Wadden Sea): spatial patterns, relationships with environmental characteristics, and comparative aspects[END_REF][START_REF] Veiga | Relationship between structure of macrobenthic assemblages and environmental variables in shallow sublittoral soft bottoms[END_REF].

Benthic communities are commonly sampled via a network of stations and a sampling periodicity varying from seasonal to decadal. Such an approach allows detecting major changes in benthic macrofauna without being able to infer on the causes of changes [START_REF] Bacouillard | Long-term spatio-temporal changes of the muddy fine sand benthic community of the Bay of Seine (eastern English Channel)[END_REF]. The simultaneous sampling of several stations distributed in the area of interest makes possible to investigate temporal changes occurring at different inter-connected scales: (1) at the scale of a station, where local processes may drive the composition of communities (sedimentary changes, local anthropogenic pressures), and (2) at the scale of cluster of stations with similar composition (hereafter called assemblages), where processes operating at larger scale may drive the temporal variability (climate, diffuse and chronic anthropogenic disturbances).

Although there is a large consensus that diversity losses threaten ecosystem functioning, man-induced changes do not systematically impact diversity. Studying the temporal variability of benthic assemblages through analysis of composition and structure gains or losses seems consequently limited regarding functional considerations. For example, substantial changes in species composition and structure in space and time may not systematically results in proportional α diversity [START_REF] Whittaker | EVOLUTION AND MEASUREMENT OF SPECIES DIVERSITY[END_REF][START_REF] Dornelas | Assemblage Time Series Reveal Biodiversity Change but Not Systematic Loss[END_REF] and trait-based functional changes [START_REF] Mclean | Fish communities diverge in species but converge in traits over three decades of warming[END_REF][START_REF] De Sousa Gomes-Gonçalves | Functional stability despite anthropogenic influences on the ichthyofauna of a tropical bay[END_REF]. In this context, trait-based approaches offer a complementary integrative framework to link functional consequences to biodiversity changes [START_REF] Suding | Scaling environmental change through the community-level: a trait-based response-and-effect framework for plants[END_REF][START_REF] Boyé | Trait-based approach to monitoring marine benthic data along 500 km of coastline[END_REF].

Studies combining composition, structure and functioning assessments of assemblages in a multidimensional scale allow the detection of multi-faceted responses to environmental variability in faunal communities [START_REF] Mclean | Fish communities diverge in species but converge in traits over three decades of warming[END_REF].

The detection of future potential changes in species distribution and assemblage compositions, which can be helpful for monitoring and conservation issues [START_REF] Desroy | Macrobenthic resources of the shallow soft-bottom sediments in the eastern English Channel and southern North Sea[END_REF][START_REF] Schückel | Spatial variability in structural and functional aspects of macrofauna communities and their environmental parameters in the Jade Bay (Wadden Sea Lower Saxony, southern North Sea)[END_REF][START_REF] Veiga | Relationship between structure of macrobenthic assemblages and environmental variables in shallow sublittoral soft bottoms[END_REF], requires the establishment of baseline knowledge. In the English Channel, scientific studies have largely focused on benthic communities in the main coastal bays and estuaries [START_REF] Desroy | Macrobenthic resources of the shallow soft-bottom sediments in the eastern English Channel and southern North Sea[END_REF].

In the bay of Saint-Brieuc (western English Channel), where intertidal habitats are partially protected since 1998 by a National Nature Reserve, the analysis of diachronic changes of benthic assemblages has not been updated for 20 years [START_REF] Gros | Benthic macrofauna data set of intertidal and subtidal softbottom sediments in the Bay of Saint-Brieuc[END_REF]. Previous studies underlined the close link between benthic resources and waders which prey on macro-invertebrates [START_REF] Sturbois | Exploitation of intertidal feeding resources by the Red Knot Calidris canutus under megatidal conditions[END_REF][START_REF] Ponsero | How do macrobenthic resources concentrate foraging waders in large megatidal sandflats? Estuarine[END_REF], highlighting the need to study and conserve benthic habitats.

After two decades, an analysis of benthic assemblages using the same protocols as those used for the definition of historical baselines was needed to assess temporal changes of intertidal habitats in this bay representative of coastal areas characterised by range of impacts due to the complex co-occurrence of management's mosaic and anthropogenic uses [START_REF] Ellis | Broad-scale disturbance of intertidal and shallow sublittoral soft-sediment habitats; effects on the benthic macrofauna[END_REF]. Natural and anthropogenic contexts (Supplementary material, Appendix A) make the study site suitable to contribute to the debate about the definition of impact assessment at spatial and temporal scales and the feasibility of such research. It is also essential to provide up-to-date knowledge to policy makers and managers about the evolution of the composition, structure and functioning of benthic assemblages.

This study, based on a long-term survey in the intertidal soft-bottom sediments of the bay of Saint-Brieuc, aims to analyse the taxonomic and functional variability of benthic communities over the last 30 years at different spatial scales using 42 stations sampled in 1987, 2001 and 2019. The approach specifically focused on the detection of trajectory patterns by applying the recent Community Trajectory Analysis Framework [START_REF] De Cáceres | Trajectory analysis in community ecology[END_REF]Sturbois et al., 2021b). Results were discussed regarding potential environmental and anthropogenic driving factors and some propositions were made for such long-term study approaches.

Material and methods

Study area

Fieldwork was conducted in the bay of Saint-Brieuc, France (Figure 1). The study area includes about 2900 ha of tidal flats divided in two coves (Yffiniac and Morieux) mainly dominated by fine to medium sands. The most sheltered parts are dominated by salt marshes. The bay is under the influence of a semi-diurnal megatidal regime.

Tidal range varies between 4 m at neap tides and nearly 13 m during spring tide (Figure 2). Divagation of channels is particularly dynamic, with increasing magnitude in low levels (Supplementary material, Appendix B). A large part of the bay (1140 ha) is included in the National Nature Reserve of the bay of Saint-Brieuc, created in 1998 to protect wader populations (Figure 1). The study area undergoes a variety of anthropogenic pressures (Supplementary material, Appendix A): sediment extraction conducted until 2004 (1000m 3 .y -1 , around stations 3 and 5) and muddy sediment deposits (near station 11) dredged from the Saint-Brieuc harbour since 2004 (120 000 m 3 .y -1 ). The bay also suffers from eutrophication, resulting in strong ephemeral macroalgae proliferation and cyclic green tides episodes. Professional (Cerastoderma edule) and recreational (Cerastoderma edule, Pecten maximus, Arenicola marina) fishing activities has been an important activity for many decades, especially in Yffiniac cove. In the lower level of Morieux cove, mussel farming on wooden poles, called "bouchots" (Figure 1) has increased over the study period, from 277 ha in 1987 to 386 ha in 2001 and 312 ha in 2019.

Field sampling

Macrofauna was sampled during winter (March) in 1987(March) in , 2001(March) in and 2019 at 42 at 42 stations, following the sampling grid (Figure 1) and the protocol defined in 1987 by [START_REF] Gros | Benthic macrofauna data set of intertidal and subtidal softbottom sediments in the Bay of Saint-Brieuc[END_REF]. At each site, four replicates were collected for macrobenthos analysis using a 1/32m² hand corer (diameter: 20cm), at a depth of 25 cm. The content of the cores was gently sieved on site through 1-mm square mesh sieve. Retained material was preserved for analysis in 5% buffered formaldehyde.

Density per m² of Arenicola marina, which is difficult to sample with a handcorer, was estimated visually by counting worm casts at each station within 10 quadrats of 1 m².

A visual description of the sediment type was recorded and two samples of sediment were collected from additional cores (diameter 5 cm, depth 5 cm) and subsequently analysed for grain size distribution, and water and organic matter contents. All replicates were collected at the four corners of a 1m² square, using GPS positionfixing (GPS Etrex Garmin).

Laboratory analysis

Macrofauna was identified to the lowest possible level (usually species level), counted, and abundances were expressed as densities per 1 m 2 . Taxon were named following the World Register of Marine Species. For granulometric analysis, sediment was desalted with freshwater during 48 h before being dried at 60°C for 48 h, sieved through AFNOR standard sieves (from 40 to 25000µm), and weighed. Water and organic matter contents were respectively determined after drying (60°C for 48 h) and combustion of sediment at 550 °C for 4 h [START_REF] Hedges | Carbon and nitrogen determinations of carbonate-containing solids [In sediments, sediment trap materials and plankton[END_REF][START_REF] Salonen | A versatile method for the rapid and accurate determination of carbon by high temperature combustion1[END_REF].

Traits collection

The functional dataset was composed of six traits divided in 23 categories (Supplementary material, Appendix C). These traits characterized the morphology (body size, flexibility, fragility) and behavioural traits [feeding behaviour, living habit, tolerance, [START_REF] Degen | The Arctic Traits Database -a repository of Arctic benthic invertebrate traits[END_REF]]. This set of traits is related to the vulnerability of species to mechanical disturbances (associated to recreational and professional fishing activity and the circulation of vehicles) and organic enrichment (tolerance). They were chosen to analyse functional changes in the community, in relation to the main pressures identified in the study area. Information was retrieved from primary literature on specific taxa, expert knowledge and three publicly available databases: polytraits (http://polytraits.lifewatchgreece.eu), the World Register of Marine Species (WoRMS Editorial Board, 2020), and Biological Trait Information Catalogue (BIOTIC, http://www.marlin.ac.uk/biotic/). Scores were attributed to species according to their affinity to each category of traits, using a fuzzy coding from 0 (no affinity) to 3 [high affinity [START_REF] Chevene | A fuzzy coding approach for the analysis of long-term ecological data[END_REF]]. Community-weighted trait values (CWMs) were calculated to characterize the functional structure [START_REF] Ricotta | CWM and Rao's quadratic diversity: a unified framework for functional ecology[END_REF].

Data analysis

Data collected in 1987, 2001 and 2019 were analysed to track taxonomic and functional changes. Analysis of sediments (sub-section 2.5.1) and α-diversity (2.5.2) were followed by spatial taxonomic clustering (2.5.3) and taxonomic and functional trajectory analysis (2.5.4). All the statistical analysis was performed with R 4.0.

Sedimentary dataset

Statistical description of sediments, based on grain-size distributions, was computed with the R package 'G2Sd' [START_REF] Fournier | G2Sd: a new package for the statistical analysis of unconsolidated sediments[END_REF]. Ten sedimentary parameters were selected to describe sediment characteristics: water content (H20), organic matter content (OMC), arithmetic mean, median (d50), gravel (> 2mm), coarse sand (1mm to 250µm), fine sand (250 to 125 µm), very fine sand (125 to 63µm), and mud (<63µm) contents.

Taxonomic α-diversity

In order to assess the different components of α-diversity, the mean number of individuals per m² (N), taxa richness (S), Shannon-Weaver index (H'), Simpson index (D) and Pielou's species evenness (J) were calculated for each station prior to any transformations of the dataset.

Spatial taxonomic clustering

Hierarchical Cluster Analysis (HCA) was performed to distinguish benthic assemblages (i.e. site groups) at each campaign by using the Hellinger distance between each pair of samples [START_REF] Legendre | Ecologically meaningful transformations for ordination of species data[END_REF] and by applying the Ward's clustering method. The indicator value method (Indval) was performed to characterize assemblages [START_REF] Dufrene | Species assemblages and indicator species: the need for a flexible asymmetrical approach[END_REF].

Detecting trajectory patterns

PCAs were performed on Hellinger transformed data and followed by Community Trajectory Analysis (CTA). CTA is based on geometrical properties of ecological trajectories [START_REF] De Cáceres | Trajectory analysis in community ecology[END_REF]. The original framework and its new extension (Sturbois et al., 2021b) 2019) developed a geometrically-based approach to trajectory resemblance which includes the shape, size, direction and position of trajectories with respect to the resemblance between all observations (state) belonging to a same trajectory. Dissimilarities between ecological trajectories (1987-2001-2019) were calculated to test their resemblance [START_REF] De Cáceres | Trajectory analysis in community ecology[END_REF] to underline potential spatial pattern with HCAs. Indicator value method [START_REF] Dufrene | Species assemblages and indicator species: the need for a flexible asymmetrical approach[END_REF] was performed on taxonomic trajectory clusters to test differences between 1987, 2001 and 2019. The relative contribution of all categories of functional traits between 1987, 2001 and 2019 were compared for each functional trajectory clusters and at the scale of the whole study area.

Results

General description of sediment

The sampling area was dominated by fine sediment at the three sampling dates:

(arithmetic means= 144.4µm ± 52.5, 135.9µm ± 44.9 and 159.5µm ± 70.6 in 1987, 2001 and 2019, respectively). Sand was the most represented grain size class over the last three decades: 97.0% ± 1.9, 98.5% ± 4.1 and 95.6% ± 5.6. Mud and gravel were less represented and respectively ranged between 1.2% ± 4.0 to 3% ± 5.5, and 0.3% ± 0.8 to 1.4% ± 5.0. Overall, surficial sediment of the study area can be summarized as fine sands organised along a continuum of fine to very fine sands characterised by local variations in mud or gravel contributions.

General description of the benthos

The number of taxa reported in the intertidal flat was variable across the sampling period (Table I): 50 in 198750 in , 53 in 200150 in and 85 in 201950 in . In comparison to 1987, 11 , 11 taxa disappeared and 14 taxa appeared in 2001 and 14 taxa disappeared and 46 taxa appeared in 2019. 36 taxa were common to the three datasets and represented 95.5% ± 2.9 of the total abundance. Polychaetes, crustaceans, and molluscs were dominant. The number of taxa of these three dominant groups increased over the study period. Other categories contributed to less than 5% of the taxa richness, except for echinoderms, which represented 8.26% of the total abundance in 2019.

The high total abundance observed in 2001 was due to a strong increase of the mollusc Donax vittatus in the low shore (Figure3), responsible for 51,8% of the total abundance.

In 1987, macrofauna was dominated by molluscs excepted in the low shore, where annelids were the most numerous group (Figure 3). In 2019, distribution between taxonomic groups was less contrasted, especially in the low shore, mostly due to an increase in mollusc abundances. Abundance of crustaceans also increased in most of the study area and a clear extension of the echinoderms distribution was observed in the western part of the study area, mostly due to the ophiuroid Acrocnida spatulispina.

Taxa richness increases from high to low shore whereas abundances follow an inverse pattern (Table I). Density, taxa richness and diversity (H, D) recorded from all samples were spatially and temporally variable. Compared to 1987, abundance (815 ind.m -² ± 1617 vs 1206 ind.m -² ± 1530), taxa richness (10.42 ± 3.21 vs 16.80 ± 6.69) and species diversity slightly increased in 2019.

Benthic assemblages

The hierarchical classification separated four clusters in 1987 and 2001 and five clusters in 2019. Mainly influenced by bathymetry (Figure 2), the distribution and structuration of assemblages was stable from upper to lower flat, despite the southwestward extension of assemblage IV, and the apparition of a new assemblage in 2019 (Figure 4). II) are detailed in Supplementary material, Appendix D.

Detecting trajectory patterns

Taxonomic trajectories

The first period was characterized by less change (total S1 lengths = 453.05, mean= 10.79 ± 3.46) than the second [total S2 lengths= 582.34, mean= 13.87 ± 5.56 (Figure 5 III) are detailed in Supplementary material, Appendix D.

Functional trajectories

Functional trajectory paths were, for all sites, similar for both periods, revealing similar level of functional variability between 1987 and 2001 (254.55, 6.06 ± 2.88) and 2001 and 2019 (223.41, 5.32 ± 2.61). Seven stations concentrated 30 % (143.67) of the overall trajectory path, illustrating a high functional variability in these areas over the last 30 years. RDT revealed recovering trajectories for 25 stations (60 %), evenly distributed in the study area (Figure 6,A). The low net change ratio (46.78% ± 0.11) was influenced by recovering trajectories.

The HCA performed on CTA of trajectory similarities separated four main clusters according to the shape of functional trajectories (Figure 6, B andC).

Cluster E grouped six stations located in sheltered and upper parts of the shore and characterised by the highest net changes (7.50 ± 1.27) and moderate trajectory path (10.80 ± 1.09). Cluster E trajectories were mainly defined by an increase in the contribution of very sensitive, free living and small species (Figure 7 The main functional shift in cluster H concerned a positive shift of suspensionfeeders, free living, medium fragile, and rigid species, notably attributable to strong variation of D. vittatus observed at these stations.

Despite some local functional shift depending on functional trajectory clusters, and some slight significant shifts for some categories of traits, the overall functional properties remained stable over time at the site scale (Figure 8).

Discussion

Local ecological dynamics influence benthic assemblages

As expected in such coastal environment [START_REF] Mclachlan | The Ecology of Sandy Shores[END_REF], macrofauna was dominated by molluscs, polychaetes, and crustaceans, and 95.5% of the total abundance was due to 36 species common to the three sampling periods. This revealed a strong stability in the pool of species influencing community abundance and composition. Some gains and losses of occasional and rare species influenced variations of the relative contribution of taxonomic group to overall richness. Changes in assemblages mainly involved some of the prevalent species, suggesting a strong impact of fluctuating population dynamics rather than species replacement in the community. Different indicators and numerous species were involved in the definition of assemblages and taxonomic trajectory clusters but the contribution of some species to these descriptions was constant over time, notably P.ulvae, L. balthica, U. [START_REF] Beukema | Some like it cold: populations of the tellinid bivalve Macoma balthica (L.) suffer in various ways from a warming climate[END_REF][START_REF] Genelt-Yanovskiy | Trends and drivers of Macoma balthica L. dynamics in Kandalaksha Bay, the White Sea[END_REF][START_REF] Jansen | Macoma balthica in Spain, a few decades back in climate history[END_REF][START_REF] Ong | Physiological responses to ocean acidification and warming synergistically reduce condition of the common cockle Cerastoderma edule[END_REF][START_REF] Philippart | Climate-related changes in recruitment of the bivalve Macoma balthica[END_REF][START_REF] Van Colen | Biogeographic vulnerability to ocean acidification and warming in a marine bivalve[END_REF]. The bay of Saint-Brieuc being located in the southern part of L. balthica's distribution, populations are likely to be affected by climate warming. At some stations, local taxonomic variability was of higher magnitude than the withinassemblages dispersion observed in 1987, which resulted in changes in the affectation of assemblage over time. Depending on stations, differences between assemblages were driven by differences in species compositions and/or dynamic of densities (see differences in table II). Some stations were also located at the border of some assemblages implying that they can be influenced by one assemblage or another depending on the year. Because changes in species composition and/or abundance (maybe favoured by hydrodynamics and/or sediments), some stations move from one assemblage to another. For instance, the station 15 located at the border of two assemblages in an area of high tidal channels dynamic was characterised by successive changes over time: some species were sampled at least at two dates (e.g.

Cerastoderma edule, Sigalion mathildae, Macomangulus tenuis or Bathyporeia sarsi)

whereas others rised high abundance at only one date (e.g. Scolelepis squamata in 1987; Macomangulus tenuis and Polycirrus medusa in 2001; Bathyporeia elegans, Cumopsis longipes in 2019).

At the opposite, others trajectories buffered local functional changes in some stations at the assemblage scale, especially in the low part of the shore. Taken separately, stations were characterized by different functional trajectories, and the concentration of local recovering trajectories in these assemblages contributed to explain the overall functional persistence. Macrobenthic communities are often associated to a mosaic of micro-habitats where natural disturbance events result in the creation of patches characterized by different recovering stages [START_REF] Defeo | Threats to sandy beach ecosystems: A review[END_REF][START_REF] Ellis | Broad-scale disturbance of intertidal and shallow sublittoral soft-sediment habitats; effects on the benthic macrofauna[END_REF].

According to [START_REF] Grimm | Babel, or the ecological stability discussions: an inventory and analysis of terminology and a guide for avoiding confusion[END_REF], persistence is defined as a stability property of ecological systems which nevertheless integrate moderate temporal variations. It differs of the concept of constancy, which describes unchanged systems. Despites changes observed at the scale of assemblages coupled with fluctuating population dynamics of structuring species, we assume an overall persistence of the whole benthic community in the bay of Saint-Brieuc. Such assumptions are in accordance with [START_REF] Bacouillard | Long-term spatio-temporal changes of the muddy fine sand benthic community of the Bay of Seine (eastern English Channel)[END_REF], which suggest that local turnover can be very high, especially in response of hydrologic and sedimentary dynamics [START_REF] Desroy | The intraannual variability of soft-bottom macrobenthos abundance patterns in the North Channel of the Seine estuary[END_REF], without adversely affecting the long-term persistence of the community.

Low shore concentrates taxonomic changes

The CTA approach revealed that the taxonomic variability increased following an inshore-offshore gradient. In the upper flat, sheltered conditions, strong influence of salinity variations and emersion time limit the number of species and the variability due to species gain. Inversely, the taxa richness increasing towards the lower flat contributes to the magnitude of change in low levels. As observed in 2019 for Aponuphis bilineata and Nucula nitidosa, the local presence of subtidal species on the intertidal flat constitutes a source of variation.

The northwest part of the bay, characterised by professional and recreational fishing activities known to impact benthic community properties [START_REF] Boldina | Effect of long-term mechanical perturbation on intertidal soft-bottom meiofaunal community spatial structure[END_REF][START_REF] Mosbahi | Immediate Effect of Clam Harvesting on Intertidal Benthic Communities in the Mudflat Zones of Kneiss Islands (Central Mediterranean Sea)[END_REF], concentrated the highest net changes. Tides associated to hydrodynamics, particularly strong in our study area, is recognized as a dominant factor shaping the variability of coastal systems [START_REF] Monbet | Control of Phytoplankton Biomass in Estuaries: A Comparative Analysis of Microtidal and Macrotidal Estuaries[END_REF][START_REF] Wildish | Factors controlling marine and estuarine sublittoral macrofauna[END_REF][START_REF] Wildish | Tidal Energy and Sublittoral Macrobenthic Animals in Estuaries[END_REF]. This factor is responsible for surficial sediment coverage and, consequently, for the dynamics of the distribution of macrofauna as demonstrated by [START_REF] Fujii | Spatial patterns of benthic macrofauna in relation to environmental variables in an intertidal habitat in the Humber estuary, UK: Developing a tool for estuarine shoreline management[END_REF] in the Humber estuary (UK), [START_REF] Compton | Distinctly variable mudscapes: Distribution gradients of intertidal macrofauna across the Dutch Wadden Sea[END_REF] in the Wadden sea, or [START_REF] Veiga | Relationship between structure of macrobenthic assemblages and environmental variables in shallow sublittoral soft bottoms[END_REF] on North Portuguese coasts. The influence of intertidal channels divagation (within a 1420 ha area since 1988, Supplementary material, Appendix A) must also be considered as a driver of sediments resuspension and habitat resetting, particularly important in the low western shore. In the eastern part, low net changes associated with high trajectory path and recovering trajectories could be explained by sand banks displacement from north to south. Such impact of sediment transport has already been highlighted in different coastal areas in New Zealand [START_REF] Thrush | Habitat change in estuaries: predicting broad-scale responses of intertidal macrofauna to sediment mud content[END_REF], France [START_REF] Gouleau | Short-and long-term sedimentation on Montportail-Brouage intertidal mudflat, Marennes-Oléron Bay (France)[END_REF] or experimentally in the UK [START_REF] Whomersley | Response of intertidal macrofauna to multiple disturbance types and intensities -An experimental approach[END_REF]. In this way, we expect any anthropogenic activity potentially modifying sedimentary parameters (bio deposition, sediment deposition, beach nourishment, dredging…) to, in turn, modify intensely the composition, structure and functioning of soft bottom assemblages as suggested by [START_REF] Veiga | Relationship between structure of macrobenthic assemblages and environmental variables in shallow sublittoral soft bottoms[END_REF]. [START_REF] Whomersley | Response of intertidal macrofauna to multiple disturbance types and intensities -An experimental approach[END_REF] suggested that benthic communities frequently disturbed by sediment movements or other natural environmental conditions would show greater resilience against further physical disturbance. Such resilient responses to disturbance, due to the selection of adapted species, depend on the history of prior disturbance and on the inherent ecological plasticity exhibited by species [START_REF] Davic | Linking Keystone Species and Functional Groups: A New Operational Definition of the Keystone Species Concept[END_REF]. In our study, 36 prevalent species contribute to the persistence of the community over time. Most of these species are adapted to living in mobile sediments [START_REF] Mclachlan | The Ecology of Sandy Shores[END_REF] and facing erosion and sedimentation.

They could then persist because of their tolerance to the harshness of the environment as suggested by [START_REF] Callaway | Historical Data Reveal 30-Year Persistence of Benthic Fauna Associations in Heaviy Modified Waterbody[END_REF].

In many coastal ecosystems, multiple disturbance sources, both natural and anthropogenic, are often confounded [START_REF] Whomersley | Response of intertidal macrofauna to multiple disturbance types and intensities -An experimental approach[END_REF], which makes difficult the definition of their respective contribution [START_REF] Harris | Anthropogenic Threats to Benthic Habitats[END_REF]. Our study confirms the complexity to disentangle factors of changes in dynamic areas, where strong natural variability results in regular community resets [START_REF] Defeo | Global patterns in sandy beach macrofauna: Species richness, abundance, biomass and body size[END_REF], which may overshadow anthropogenic impact. On the opposite, areas characterized by low hydrodynamics where benthic communities are not naturally adapted to disturbances would therefore be very sensitive to additional physical disturbance as showed by [START_REF] Jac | Detecting adverse effect on seabed integrity. Part 2: How much of seabed habitats are left in good environmental status by fisheries?[END_REF]. In the bay of Saint-Brieuc characterised by severe natural constrains, further studies coupling experimental and high frequency field samplings should deal with the respective role of hydrodynamics and anthropogenic activities in resetting benthic communities (shear stress measurement in area concerned and not by the circulation of vehicles from mussels farming, analysis of the benthic macrofauna variability coupled with different gradients of fishing activities and or local dynamics…).

Functional properties of the community persist over time

The CTA resulted in contrasted taxonomic and functional trajectory patterns (i.e. net changes, departing and recovering trajectory). [START_REF] Clare | Community variability and ecological functioning: 40 years of change in the North Sea benthos[END_REF] showed that density compensation by redundant taxa buffers ecological functioning changes, but that functional stability is subjected to aperiodic disruption due to substitutions of dissimilar taxa or uncompensated population fluctuations. Accordingly, high trajectory paths observed in the low shore reflected an important functional variability, but recovering trajectories notably influenced by D. vittatus fluctuating population dynamics led to low net changes. Functional properties persisted over time at the scale of the whole community. Low sensitivity of traits such as body size, flexibility, fragility to genus and, in a lesser extent, species changes also contribute to limit the influence of taxonomic variability on the functional stability. Lower influence of traits relative to feeding behaviour finally suggested no major changes at the scale of intertidal sedimentary habitats, in accordance with (Sturbois et al., 2021b) who showed that high sedimentary variability was limited to few sites.

To choose the set of traits in relation to main pressures, we hypothesised that fishing activities in the lower flat may influence benthic assemblages over time according to changes in relative values of body size, flexibility, fragility. Low net changes in the lower part of the shore for the widespread functional trajectory cluster G and the influence of D. vittatus in the functional trajectory of cluster H does not support such a functional response to human pressures.

A similar interest was devoted to the tolerance of species to organic matter enrichment according to the eutrophication process, mainly due to green tides. No bidentata were stimulated in soft bottom sediments impacted by green tides, contrary to the crustaceans B. elegans and Cumopsis spp., which were impacted. In our study, the abundance and occurrence of these four taxa increased over the study period, revealing a more complex gradient of sensitivity to green tides and/or nutrient enrichment and the implication of other driving factors. As suggested by [START_REF] Bolam | Effects of macroalgal mats on intertidal sandflats: an experimental study[END_REF], [START_REF] Everett | Macroalgae in marine soft-sediment communities : effects on benthic faunal assemblages[END_REF][START_REF] Raffaelli | Ecological impact of green macroalgal blooms[END_REF], generalisations on the effect of algae mat are not straightforward. In the study area, we believe that the important green-tides associated to nutrient uptake buffer the influence of nutrient enrichment on macrofauna.

The crucial need of data at multiple spatial and temporal scales

The CTA applied to marine communities offers new trajectory metrics and tools [START_REF] De Cáceres | Trajectory analysis in community ecology[END_REF]Sturbois et al., 2021b), useful to the accurate quantification and description of changes as suggested by [START_REF] Cimon | Impact of multiple disturbances and stress on the temporal trajectories and resilience of benthic intertidal communities[END_REF], which claimed for the development of new metrics to document effects on community dynamics.

The potential and the interest of such framework focused on the analysis of ecological dynamics increases significantly with the size of the time series.

There is virtually no appropriate benchmark data describing natural standards and alterations of most environmental systems concerned by impacts accumulation [START_REF] Ellis | Broad-scale disturbance of intertidal and shallow sublittoral soft-sediment habitats; effects on the benthic macrofauna[END_REF]. Like many estuarine and coastal areas, the bay of Saint-Brieuc has been exposed to strong human influence over the last decades. The ecological state identified in 1987 and considered as the initial state in this study, was defined under a pressure context already well established (notably nutrient enrichment, macroalgae proliferation, mussel farming, fishing activities). The absence of a pristine reference state constitutes a strong limitation in the analysis of long-term trends [START_REF] Hardman-Mountford | Diagnostic monitoring of a changing environment: An alternative UK perspective[END_REF][START_REF] Rubal | Effects of subtle pollution at different levels of biological organisation on species-rich assemblages[END_REF][START_REF] Veiga | Relationship between structure of macrobenthic assemblages and environmental variables in shallow sublittoral soft bottoms[END_REF]. Furthermore, our study suffers from the absence of long-term measurements of environmental factors and anthropogenic pressures, especially to disentangle natural vs anthropogenic changes.

Our results illustrate the value of historical data for marine and coastal management, as well as the importance of maintaining consistent protocols, as underlined by [START_REF] Callaway | Historical Data Reveal 30-Year Persistence of Benthic Fauna Associations in Heaviy Modified Waterbody[END_REF]. Overall, the complementary methods and scales used to analyse the data set contributed to the research frameworks focused on the definition of changes and impact assessments at spatial and temporal scales in marine ecosystems.

However, it pointed the limitations of sampling strategies that only imply long-term networks. The patterns of changes described were based only on three sampling occasions, which limits the understanding of processes occurring between each surveys, and potentially overshadowing important intermediate changes (e.g. occasional dominance by very large numbers from very few species as observed for D. vittatus). As suggested by [START_REF] Bacouillard | Long-term spatio-temporal changes of the muddy fine sand benthic community of the Bay of Seine (eastern English Channel)[END_REF], the coupling of different spatial and temporal scales in sampling strategy [few stations with high frequency sampling [START_REF] Hewitt | Multiple stressors, nonlinear effects and the implications of climate change impacts on marine coastal ecosystems[END_REF] vs low frequency larger sampling networks [START_REF] Kröncke | Changes in North Sea macrofauna communities and species distribution between 1986 and 2000[END_REF]] could help to track changes between long-term networks. [START_REF] Bacouillard | Long-term spatio-temporal changes of the muddy fine sand benthic community of the Bay of Seine (eastern English Channel)[END_REF] propose in the bay of Seine 60 stations sampled every five years to assess long-term changes in the beta diversity (i.e. assemblages distributions) and five stations sampled An overarching question in conservation biology is the assessment of marine protected area effect on marine habitats [START_REF] Appolloni | Does full protection count for the maintenance of β-diversity patterns in marine communities? Evidence from Mediterranean fish assemblages[END_REF][START_REF] Bevilacqua | Mitigating human disturbance: can protection influence trajectories of recovery in benthic assemblages?: MPAs and benthic assemblage recovery[END_REF][START_REF] Day | Categories to Marine Protected Areas (Best Practice Protected Area Guidelines Series No. 19[END_REF][START_REF] Klein | Shortfalls in the global protected area network at representing marine biodiversity[END_REF][START_REF] Thompson | Multi-Scale Sampling to Evaluate Assemblage Dynamics in an Oceanic Marine Reserve[END_REF]. In this study, it was not possible to assess accurately any conservation effect of the National Nature Reserve as 1) species and assemblages remained locally very contrasted between unprotected and protected areas and 2) such assessment ideally requires a global network of protected and unprotected areas.

However, this study provides an up-to-date assessment to policy makers and managers about the evolution of the composition, structure and functioning of benthic assemblages. This knowledge will integrate future governance rounds, in a multiscale approach involving policy makers and stakeholders, and will make marine habitats, a central part of the conservation process [START_REF] Greathead | A generic framework to assess the representation and protection of benthic ecosystems in European marine protected areas[END_REF], which aim notably to extend the current natural reserve perimeter. Such approaches should not only focus on aggregating local present day practices, but also on exploring future visions and potential scenarios (Pereira L., 2021).

Conclusion

The CTA framework brings new quantitative and qualitative insights to characterize and illustrate changes in marine habitats. The consideration of both local changes and variations at population levels helps to explain the differences observed at the scale of assemblages. Observed changes mainly reflected random population dynamics of structuring prevalent species in habitat under strong natural constrains, rather than strong specific changes of the community. Few species may underline the long-term effect of green tides in some assemblages. This study also claims for ambitious monitoring strategy in order to integrate multiple spatial and temporal scales, and efficient measurements of environmental and anthropogenic driven factors.

Tables:

Table I: General description of macrofauna

  were used to study and represent temporal changes at the station scale taking into account all the dimensions of the multivariate analysis through trajectories of surveys in the multivariate space. A set of distance-based metrics were calculated. Segment-length represents, for a station, the distance between two consecutive surveys. S1 informs about changes between 1987 and 2001, and S2 between 2001 and 2019. The trajectory path represents the sum of segment length (S1+S2). Net changes were calculated to analyse trajectories between 1987 and 2019, considering the 1987 data set as an initial state (i.e. the first ecological state of the time series). Net changes informs about the overall change occurring during the study. The net change ratio represents the ratio between the net changes and the trajectory path. It informs about the straightness of recovering and or departing processes with respect to the initial state. Recovering and departing consecutive trajectory segments (RDT) were identified by subtracting Net change n-1 to Net change n (Sturbois et al., 2021b). De Cáceres et al. (

  Assemblage I, located in muddy sands on the upper part of the flat, was characterized Pygospio elegans in 1987, and Peringia ulvae in 2001 and 2019 (first rank indval, Table II). P. ulvae was dominant in 1987 and 2019, whereas Corophium arenarium was the most abundant species in 2001. Assemblage II, associated to very fine to fine sand of the upper mid-flat, was characterized by Macomangulus tenuis in 1987 and 2019 and Urothoe poseidonis in 2019. M. tenuis remained the most abundant species over time. Assemblage III, associated to very fine to fine sand of the lower mid-flat, was characterized by U. poseidonis in 1987 and Acrocnida spatulispina in 2001, and 2019. This assemblage was dominated by M. tenuis in 1987, U. poseidonis in 2001, and A. spatulispina in 2019. Assemblage IV, associated to sands mainly located in low tide levels, was characterised by Scoloplos armiger in 1987, Donax vittatus in 2001, and Magelona mirabilis in 2019. It was dominated by P. medusa in 1987 and D. vittatus in 2001 and 2019. A new assemblage was identified in 2019, in the upper part of the shore, on bare sands associated to major sand intakes. Eurydice affinis characterized this assemblage and Bathyporeia pilosa was the most abundant species. Distribution maps of species used to describe assemblages (Table

  , A)]. Highest net changes occurred in the low shore where eight stations concentrated 31% of changes (179.24). A positive correlation (0.529, p=0.0002) was found between net change and latitude. RDT was positive for 31 stations (71.43%) and negative for 11 stations, showing an overall departure dynamics from the initial taxonomic state. Total net change (582.34) represented only 56 % of the total trajectory path (1035.39). This low net change ratio (53.14% ± 0.12) indicated significant direction changes (i.e. driving species groups) in the multivariate space between 1987-2001 and 2001-2019 for departing trajectories, and included recovering ones. The HCA performed on CTA of trajectory similarities separated four main clusters of taxonomic trajectories (Figure 5, B and C). Cluster A was composed of 8 stations from the sheltered part of the bay, and describing the evolution from an assemblage characterized P. elegans in 1987 and B. pilosa in 2019 and dominated by P. ulvae in 1987 and 2019 and C. arenarium in 2001. Cluster B, the most common group (28 stations), was characterised by M. tenuis in 1987, U. poseidonis in 2001 and E. dolfusi in 2019. M. tenuis remained the most numerous species over time. Cluster C (NC=25.06 ± 0.38; TP=35.33 ± 0.83) and D (NC=20.78 ± 1.91; TP=40.65 ± 5.28) were characterised by longer net changes and trajectory path compared to cluster A (NC=10.49 ± 0.83; TP=21.20 ± 1.43) and B (NC=12.74 ± 0.81; TP=23.94 ± 1.11), but represented a limited number of trajectories (14 % of stations). Cluster C grouped 4 stations located in the north-west, characterized by C. gibber at the three campaigns and dominated by P. medusa in 1987 and D. vittatus in 2001 and 2019. Cluster D which regrouped only two stations located in the eastern low shore, was characterised by Bathyporeia guilliamsoniana in 1987, P. medusa in 2001 and Nucula nitodosa in 2019. Magelona mirabilis was the most numerous species in 1987 and 2019 and D. vittatus remained dominant in 2001. Distribution maps of species used to describe taxonomic trajectory clusters (Table

  ). Cluster F, which concerned only one station, was characterised by the highest trajectory path (27.29) but one of the smallest net change (2.63) in relation with an effective recovering trajectory between 1987-2001 and 2001-2019. Cluster G was the most common pattern of functional trajectories (30 stations) characterized by moderate net changes (4.66 ± 0.28) and trajectory path (9.53 ± 0.48). This cluster was characterised by moderate functional changes, except for an increase in the contribution of flexible and small size species. Cluster H grouped five stations of the low shore characterised by high trajectory path (19.97 ± 1.72) and moderate net changes (5.71 ± 0.94), in relation with recovering trajectories occurring at this sites.

  poseidonis, C. edule, M. tenuis, C. gibber and M. mirabilis. At the scale of stations, despite an overall departing taxonomic trend, we observed mainly recovering functional trajectories. Two main mechanisms are involved in the long-term variability of benthic assemblages: the dynamics of structuring species and the existence of strong local variability. 1. Coupled with random gain or losses of species (Legendre, 2019), the dynamics of structuring species results in moderate specific taxonomic changes in the multivariate structure of assemblages. However, several species can induce strong trends, such as D. vittatus, whose abundance varied from 13.4 ± 15 ind.m² (20 stations) in 1987 to 3107 ± 6347 ind.m² (16 stations) in 2001 and 72.8 ± 101 ind.m² (19 stations) in 2019 or A. spatulispina and Kurtiella bidentata, whose occurrence and abundance increased. In assemblage I, the decline in L. balthica populations is potentially attributable to climate change (temperature increase) as suggested by different works in European regions

2.

  The stronger local variability is able to induce changes in assemblage definition, according to the initial state defined in 1987. The identification of a new assemblage in 2019, in relation to sedimentary variations, and the re-affectation of some stations to different assemblages illustrate such spatial changes. The distribution of A. spatulispina and K. bidentata in particular significantly extended in the sampling area over time.

  major change was observed at the scale of the whole study area in the contribution of the different categories of tolerance. In accordance, changes were moderate in the different functional clusters except in cluster E characterized by an important positive shift in the contribution of very sensitive species. Despite a slow long-term water quality improvement trend and a progressive decrease of green tides associated to local water development and management plan, these events still strongly affect the study area. Quillien et al. (2015a) discussed the potential influence of Ulva sp. mats on abundance peak of D. vittatus (improved recruitments and negative effect on predators of larvae), which could explain high densities reaching 21 487 ind.m² in 2001 and 14685 ind.m² in 2011 (unpublished data). Quillien et al. (2015b) also highlighted large-scale effects of green tides on macrotidal sandy beaches in Brittany, showing notably that the echinoderm A. spatulispina and the mollusc K.

  at least once a year to analyse the heterogeneity in the responses of local community to multiple stressors. Sturbois et al. (2021a) have also demonstrated such interests of coupling different spatial and temporal scales by the use of two complementary data sets (38 stations in 1987 and 2019, and one station from 2005 to 2019). In this study, the 15 years time series suggested that changes observed in a subtidal benthic community over the last 30 years occurred probably from 2015. It claims for the local development of an ambitious monitoring strategy in the intertidal area based on the 42 stations network sampled every 10 years coupled with at least one station sampled once a year in each assemblage.

Figure 1 :

 1 Figure 1: Location of the study area and sampling strategy. Black dots corresponds to the 42 stations sampled in 1987, 2001 and 2019. The North limit of the Nature reserve is indicated by a black line and mussel bouchots by polygons.

Figure 2 :

 2 Figure 2: Bathymetry of the study area. Isolines are coloured according to altitude from the top (red) to the low tide level (blue). Mussel bouchots are represented by polygones.

Figure 3 :

 3 Figure 3: Variations in relative abundance of the main numerous taxonomic group between 1987 (A), 2001 (B) and 2019 (C). Size of pies corresponds to total abundance and sections to taxonomic group contributions.

Figure 4 :

 4 Figure 4: Spatial distribution of benthic assemblages identified with Hierarchical Cluster Analysis performed on 42 stations in 1987, 2001 and 2019.

Figure 5 :

 5 Figure 5: Taxonomic trajectories. A.: Maps of local trajectories. Distance based metrics are represented with circular and triangle symbols. Net changes are represented with blue circles between 1987 and 2019. Bottom triangles represent the trajectory segment S1 (1987 to 2001) and top ones S2 (2001 to 2019). Size of symbols corresponds to lengths. For triangles, colors are used to distinguish recovering (black) from departing trajectories (grey). B: Trajectory ordination diagram. The three ecological state of each station were connected with arrows in order to represent the whole trajectory path. Each trajectory path was coloured according to taxonomic trajectory clusters defined with Hierarchical Cluster Analysis based on CTA trajectory similarity analysis. Despite trajectory distance-based metrics and dissimilarity were calculated considering all of the dimensions of the PCA, the first two dimensions show only 14% of the total variation, which may potentially hide some patterns. C.: Localisation of trajectory clusters within the study area.

Figure 6 :

 6 Figure 6: Functional trajectories. A.: Maps of local trajectories. Distance based metrics are represented with circular and triangle symbols. Net changes are represented with green circles between 1987 and 2019. Bottom triangles represent the trajectory segment S1 (1987 to 2001) and top ones S2 (2001 to 2019). Size of symbols corresponds to lengths. For triangles, colors are used to distinguish recovering (black) from departing trajectories (grey). B: Trajectory ordination diagram. The three ecological state of each station were connected with arrows in order to represent the whole trajectory path. Each trajectory path was coloured according to functional trajectory clusters defined with Hierarchical Cluster Analysis based on CTA trajectory similarity analysis. C.: Localisation of trajectory clusters within the study area.

Figure 7 :

 7 Figure 7: Barplots of Community-weighted trait values (CWM) values shift for each categories of traits (mean 2019mean 1987). Bars are coloured according to functional trajectory clusters.

Figure 8 :

 8 Figure 8: Box plots of Community-weighted trait values (CWM) traits values in levels of respective categories according to year (1987: red; 2001: green; 2019: blue). Box plot indicates the median (line in the box), the first and third quartiles (hinges) and extreme data points (black points). Significance of student permutation paired tests are indicated for each category traits (ns: non significant; *: pval<0.05; **: pval<0.01; *** pval:<0.001).

  Table II: Characteristics of assemblages identified with the Hierarchical Cluster Analysis in 1987, 2001 and 2019: Assemblages (Ass), number of station (N st), Taxa richness (S), Mean density per m² (M Dens), Shannon-Weaver Index (H), Simpson Index (D) and Pielou's evenness (J) (± = standard error). Distribution of species used to describe assemblages are detailed in Supplementary material, Appendix D. Table III: Characteristics of ecological states of trajectory clusters in 1987, 2001 and 2019. Taxonomic trajectory are coloured according to trajectory clusters and described with distance-based metrics mean ± se (NC: net change, TP: Trajectory path). Distribution of species used to describe taxonomic trajectory clusters (Ind.sp.: Indicative species, 1st. Num.: two first numerous species) are detailed in Supplementary material, Appendix D.

Acknowledgments

We thank the three anonymous reviewers of previous versions of this manuscript for their constructive comments. We acknowledge all people involved in field sampling and analysis with a special mention for Alexandre Le Moal and Charles Veilllard, which significantly contribute to laboratory analysis.

Funding sources

This work was supported by Agence de l'Eau Loire-Bretagne (grant number 180212501), the Région Bretagne (OSIRIS PFEA621219CR0530023), the Europe for the European maritime and fisheries fund (grant number FEAMP 621-B) and the Ministère de la Transition Ecologique et Solidaire (grant number EJ N°2102930123) who fund the ResTroph Baie de Saint-Brieuc research program. This work was carried out as part of the PhD thesis of A. Sturbois for Université de Bretagne Occidentale.

Supplementary material:

Appendix A: Time line of sampling dates with respect to main anthropogenic uses,