Bivalve monitoring over French coasts: multi-decadal records of carbon and nitrogen elemental and isotopic ratios as ecological indicators of global change

Recent changes in climate and environment, influenced by both global and local factors, have had profound impacts on coastal ecosystem functioning and trajectories. By examining archived samples from ongoing ecological monitoring efforts, particularly focusing on bivalves like mussels and oysters, we gain a valuable long-term perspective on how ecosystems are responding at various scales. We conducted analyses on carbon and nitrogen content (C %, N %) as well as elemental and isotopic ratios (C:N, δ¹³C, δ¹⁵N) of mussel and oyster soft tissues collected annually at 33 sites along the French coast from 1981 to 2021. This extensive dataset (https://doi.org/10.17882/100583, Liénart et al., 2024a) offers a comprehensive view spanning multiple decades and ecosystems, allowing us to track how coastal ecosystems and marine species record changing climate, physical-chemical environments, and organic matter cycles. Additionally, these data are crucial for establishing isotope baselines for studying food webs. Ultimately, this dataset provides valuable information for more effective ecosystem conservation and management strategies in our rapidly changing world.

1 Introduction

Over the past 2 decades, rapid and sometimes abrupt changes in climate and environmental conditions have significantly impacted ecosystem functioning with human activities widely acknowledged as a primary factor driving and modulating these changes (Cloern et al., 2016). At the land–ocean interface, coastal and estuarine environments are particularly vulnerable to climate change alongside local human-induced pressures such as pollution and habitat destruction (Harley et al., 2006; Lotze et al., 2006). In such dynamic environments, pressures vary simultaneously, with different directions and magnitudes, over different time (short- vs. long-term) and spatial (global vs. local) scales (Cabral et al., 2019). Hence, a fundamental challenge persists in elucidating how the various facets of global change affect coastal ecosystem responses at local and regional scales and consequently the services they provide to human societies. Recurrent ecological monitoring programs conducted at large spatial scales are valuable in assessing ecosystem status and the pressures affecting important properties (e.g., biodiversity, carbon, and nutrient cycling) and detecting changes over multiple decades (Hofmann et al., 2013; Sukhotin and Berger, 2013). Yet only a few ecological indicators are both integrative (i.e., measured in long-living species' tissues) and sensitive enough (e.g., measured at a molecular level, such as metallothionein, used as a biomarker of metal exposure; Amiard et al., 2006) to various disturbances while also exhibiting predictable responses with a low variability in its response (Dale and Beyeler, 2001; Niemi and McDonald, 2004).

Filter-feeding bivalves include important reef habitat forming species, which promote benthic–pelagic coupling and nutrient recycling (Ray and Fulweiler, 2020). Some of them, such as mussels (Mytilus) and oysters (Magallana gigas, formerly Crassostrea gigas), represent an economical value through aquaculture and are widely used as bioindicators for climate and environmental change as well as contaminant monitoring (Kanduč et al., 2018; Karlson and Faxneld, 2021; Mazaleyrat et al., 2022; Briand et al., 2023; Chahouri et al., 2023; Liénart et al., 2024b). The use of ecological indicators is the basis of ecosystem monitoring to detect early-warning signals in ecosystem changes or disturbances (Dale and Beyeler, 2001). Indeed, bivalve tissues record the environmental conditions of their sampling area over time. Since the end of the 1970s, the French monitoring network ROCCH (Réseau d'Observation de la Contamination CHimique, coordinated by Ifremer) has been monitoring chemical contamination in the environment along French coastlines using oysters and mussels as bioindicators (Briant et al., 2018; Chahouri et al., 2023). The ROCCH network focuses on chemical contaminants, but its archived tissues represent a large sample bank that can be useful for other ecological purposes, such as understanding long-term changes in ecosystem functioning and responses to global or local pressures.

The C:N ratio is mostly an indicator of bivalve condition and physiology, reflecting the balance between organisms' requirements and elemental availability in the environment; i.e., ecological stoichiometry sensu Elser et al., 2003; N content increases (thus the C:N ratio decreases) as the protein content increases whereas C content increases (thus the C:N ratio increases) as the lipid or carbohydrate content increases. When measured long-term in bivalve tissues, such a dataset allows scientists to understand ecosystem responses to the changing physical–chemical environment and organic matter cycling and to track trends in climate changes and its effect on coastal ecosystems (Liénart et al., 2020, 2024b). We took advantage of the large-scale multi-decadal sample bank of the ROCCH monitoring network to analyze C and N content as well as C:N, δ¹³C, and δ¹⁵N ratios in mussels (Mytilus edulis, Mytilus galloprovincialis) and oysters (Magallana gigas) tissues from a set of 33 stations distributed along the French coasts. In this paper, we present a unique dataset of multi-decadal and multi-ecosystem carbon and nitrogen content and elemental and isotopic ratios from three widespread bivalve species. Ultimately, it could provide valuable input for developing predictive models of bivalve physiology (Emmery et al., 2011) or trophic ecology (Marín Leal et al., 2008), explaining the ecosystem response to future environmental changes and possibly forecasting potential impacts of climate change and human activities on coastal ecosystems. Overall, this long-term dataset provides insights into ecosystem dynamics and is essential for advancing scientific understanding in the face of ongoing environmental challenges.

2 Methods and data

2.1 ROCCH: network and sampling procedure

Over the last 4 decades, the French national monitoring network for chemical contaminant, ROCCH (Réseau d'Observation de la Contamination CHimique), coordinated by Ifremer (Institut français de recherche pour l'exploitation de la mer, https://littoral.ifremer.fr/Reseaux-de-surveillance/Environnement/ROCCH-Reseau-d-Observation-de-la-Contamination, last access: 19 February 2025), has been annually sampling bivalves as bioindicators of chemical contamination. The number of stations has varied since the initiation of the network, so has the sampling frequency. Approximately 150 stations are nowadays monitored along the French coastlines, with the historical station having been monitored since 1979. Sampling takes place once a year during winter (currently mid-February), with a tolerance of one tidal cycle before and after the target date, meaning a 6-week amplitude spreading from the end of January to the beginning of March on all sites of the network. Three different species are targeted – namely, the Pacific oyster (Magallana gigas) and the blue mussel (Mytilus edulis) in the English Channel and Atlantic coastal areas, and Mytilus galloprovincialis in the Mediterranean coastal area.

The sampling protocol has been identical since the start of the monitoring and is designed to acquire bivalve samples with consistent and homogeneous characteristics for measuring chemical concentrations of contaminants. Bivalves are collected alive at fixed points (maximum tolerance of 180 m around the selected point) chosen so they are away from known anthropogenic discharges. Bivalves are sourced from either wild beds or dedicated rearing facilities, ensuring that they remain on site for at least 6 months before being sampled. The selected individuals are adults of the same species and of uniform size (30 to 60 mm long for mussels, 90 to 140 mm long for oysters, i.e., 2 to 3 years old). A minimum of 50 mussels or 10 oysters is required to constitute a representative pooled sample accounting for inter-individual variability and to get enough material for chemical analysis and long-term storage. Bivalves are first depurated for 18 to 26 h in decanted seawater collected near the collection site. Next, once extracted from the shell, the whole bodies (i.e., total soft tissues) of each individual are pooled and drained all together for 30 min. The resulting pooled bivalve samples are placed in clean glass containers (washed and baked for 8 h at 450 °C) covered with calcined aluminum foil and sealed with a plastic lid and frozen (−20 °C). Frozen samples are sent to the central laboratory in Nantes (France), where the samples are thawed, grounded, and homogenized in a stainless-steel-bladed blender and then freeze-dried. After the chemical analysis of a sample aliquot, the containers are stored indefinitely at the Nantes Ifremer center at room temperature, protected from light, and under moisture-regulation control to prevent moisture pick-up. The original protocol (Grouhel, 2023, in French) is available at https://doi.org/10.13155/97878.

Figure 1 Map of the 33 stations from the three sea coastal areas of France analyzed for bivalves' elemental and isotopic ratios. Bivalve species are the mussels Mytilus edulis and Mytilus galloprovincialis and the oyster Magallana gigas. The main rivers are presented in blue. Basemap source: Esri, HERE, © OpenStreetMap contributors, and the GIS user community.

2.2 Study sites and bivalve dataset

In the present dataset, bivalve samples were selected from 33 stations of the ROCCH network distributed along the three French sea coastal areas: the English Channel (EC, 10 stations), the Atlantic Ocean (AO, 17 stations), and the Mediterranean Sea (MS, 6 stations; Fig. 1, Table 1). The stations were selected to span over a wide diversity of geomorphological and environmental conditions encountered by bivalves along the French coastlines. Ecosystems vary between open-sea coasts (rocky or sandy littoral shores), open or semi-enclosed bays and rias, lagoons (shallow waterbodies almost closed by narrow landforms), and estuary mouths, ranging from eutrophic to oligotrophic systems, under temperate oceanic or Mediterranean climates and along gradients of influence from the main French rivers (Seine, Loire, Garonne, Dordogne, and Rhône; Table 2). The complete dataset consists of 1141 bivalve winter samples for time series spanning 11 to 40 years (depending on stations) over the period 1981–2021 (Table 1; https://doi.org/10.17882/100583, Liénart et al., 2024a).

2.3 Sample analysis and data quality

Aliquots of 400–700 µg of each dry and grounded archived sample of pooled bivalve tissues (not acidified) were analyzed for carbon and nitrogen elemental and isotopic composition (C, N, δ¹³C, δ¹⁵N) at the Center for Physical Science and Technology (Vilnius, Lithuania) with a FlashEA 1112 Series Elemental Analyzer (Thermo Finnigan) connected to a Delta V Advantage isotope ratio mass spectrometer (Thermo Fisher). Carbon and nitrogen elemental contents are expressed as the percentage (%) of mass ratio (g g⁻¹) dry weight, while the C:N ratio is express in mol mol⁻¹. Isotope ratios are expressed using the conventional delta notation: δ¹³C_sample or δ¹⁵N${}_{\mathrm{sample}}=\left[\right(R_{\mathrm{sample}}/R_{\mathrm{standard}})-\mathrm{1}]$, where $R={}^{\mathrm{13}}\mathrm{C}/{}^{\mathrm{12}}\mathrm{C}$ or ${}^{\mathrm{15}}\mathrm{N}/{}^{\mathrm{14}}\mathrm{N}$, in per mill deviation (‰) from international reference, Vienna Pee Dee Belemnite for δ¹³C and atmospheric N² for δ¹⁵N. The mass spectrometer was calibrated against external certified standards (IAEA600, caffeine: δ¹³C of −27.77 ‰, δ¹⁵N of 1.00 ‰; USGS24, graphite: δ¹³C of −16.05 ‰) and internal standards (casein: δ¹³C of −23.30 ‰, δ¹⁵N of 6.30 ‰, glycine: δ¹³C of −45.20 ‰, δ¹⁵N of 3.00 ‰), controlled using IAEA-CH-6 (sucrose); IAEA-N-1 and IAEA-N-2 (ammonium sulfate) and USGS24 (graphite) were added at the beginning and after every 10 samples within each batch of samples to control analytical performance.

Figure 2 Box plots of δ¹³C, δ¹⁵N, C, N and C:N values of bivalves for the 33 stations of the three French sea coastal areas over the period from 1981–2021 (Magallana gigas, Mytilus edulis, Mytilus galloprovincialis). The lower and upper edges of the box indicate the 25th and 75th percentiles and the middle line the median. The whiskers indicate the maximum/minimum values and the open circles the extreme values. The black arrows and associated values represent extreme values out of graphical bounds.

Analytical precision was always better than 0.2 ‰ for δ¹³C and δ¹⁵N and 0.1 mol mol⁻¹ for the C:N ratio (median of the standard deviations of each standard used for all batches is 0.12 ‰ for δ¹³C and δ¹⁵N and 0.04 mol mol⁻¹ for C:N ratio). Each sample was usually analyzed once, and for some peculiar values that seemed too high/low or outside of the average trend, samples were reanalyzed to verify the values. If values were similar (<0.5 ‰), we chose the average value of both analyses; otherwise, we chose the value associated with the best analytical precision. The dataset presented here consists of a single value for a given year at a given station (the standard deviation associated with repeated measures is not presented). One data point was considered an outlier in regard of its carbon isotope value (reanalysis confirmed the value): the station 14-PALL in 2014 (Fig. 2) exhibited unexpected low values of δ¹³C (−26.96 ‰) and high δ¹⁵N (10.26 ‰) compared to the average values at this station (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=-\mathrm{18.71}\pm \mathrm{0.75}$ ‰, ${\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=\mathrm{8.41}\pm \mathrm{0.93}$ ‰), suggesting potential issues in the sampling and/or storage. This value was therefore removed from average calculations and statistical tests. Within the C and N % datasets, there were a few exceptionally high/low values (Fig. 2) for the stations 1-AMBL in 1990 (low C % only), 5-CLHV in 1994 (high C and N %), 18-BOUR in 1990 (low C and N %), 20-RIVE in 1989 (high C and N %), and 30-STMM in 2016 (high C and N %). These values must be considered with caution when interpreting or using C and N % data and were excluded from statistical calculations (see Sect. 2.4). Nevertheless, these data pairs all led to a consistent C:N ratio at these sites and years and were kept for C:N statistical analyses. All the other values presented in this dataset are considered analytically valid and were scientifically validated taking into consideration expert knowledge of each ecosystem.

Table 1 Metadata description for the 33 stations analyzed for bivalve carbon and nitrogen elemental and isotopic ratios. The species studied are the mussels Mytilus edulis and Mytilus galloprovincialis and the oyster Magallana gigas. The station 31-CART represents two stations merged over time (at a distance of 720 m). The geographic coordinate system is WGS 84.

Table 2 Environmental characteristics of the 33 stations studied for bivalve elemental and isotopic ratios. Average values for salinity, water temperature, and chlorophyll a were calculated, when available, over the studied period from monitoring stations located nearby bivalves sampling stations, with data retrieved from Surval (https://surval.ifremer.fr, last access: November 2022, Ifremer) and Somlit (https://www.somlit.fr/, last access: November 2022, INSU) databases. Average annual flow rates of the main rivers were calculated over the study period using the HydroPortail database (https://hydro.eaufrance.fr/, last access: November 2022). The average trophic status for each station is based on the literature (Liénart et al., 2017; Lheureux et al., 2023) and the knowledge of local experts. n/a: not available.

2.4 Statistical analyses

Statistical analyses were performed with the R software (R Core Team, 2022, version 4.3.1).

Due to the nature of the dataset (i.e., only one of the three species was sampled for each sea coastal area except for stations 15-PEBE and 16-CHEM from the Atlantic coastal area), it was not possible to test whether the observed geographical differences in elemental and isotopic ratios were due to differences between species or sea coastal area, nor was it possible to ascertain if there was any interaction between these two factors. Consequently, we have tested the effects of both factors (species and sea coastal area) independently. Therefore, particular attention should be paid to interpreting the results of these tests. The effects of species (three levels) and sea coastal area (three levels) were tested independently for each elemental and isotopic variable over the entire dataset (all stations and time series) with non-parametric Kruskal–Wallis tests followed by Dunn post hoc tests (R package PMCMRplus version 1.9.10, functions kruskalTest( ) and kwAllPairsDunnTest( ); Pohlert, 2023) since normality and homogeneity of variance (inspected using Shapiro–Wilk and Levene tests, respectively) were not met, precluding the use of analyses of variance (ANOVA, parametric test).

Monotonous temporal trends in bivalve elemental and isotopic times series were examined for each station using Mann–Kendall tests corrected for autocorrelation (R package modifiedmk version 1.6, function mmkh( ), Patakamuri and O'Brien, 2021). Linear models (R package stats version 4.2.1, function lm( ), R Core Team, 2022) were applied to calculate the value of the slope (in ‰ per decade) for each time series (note that all time series do not have the same number of years). Complete time series were tested for shifts, allowing for the detection of one unique shift per time series (R package cpm version 2.3, function detectChangePoint( ), Gordon, 2015).

3 Data characteristics and general overview

This dataset illustrates the temporal and spatial variability in carbon and nitrogen content (C and N %) and elemental and isotopic ratios (C:N, δ¹³C, δ¹⁵N) of three species of bivalve along the French coasts from 1981 to 2021. Overall, the average δ¹³C and δ¹⁵N values of bivalves were significantly lower for the stations of the Mediterranean Sea compared with the stations of the English Channel and the Atlantic Ocean. We observed that stations closest to river mouths displayed the lowest average δ¹³C values compared to those more distant from river influence. Over the last 10 to 40 years, nearly all stations exhibited a significant decrease in δ¹³C and δ¹⁵N. There was no spatial nor temporal pattern in C:N, but values differed between bivalve genera (mussel vs. oyster). It is important to note that, due to the nature of the dataset (one species present for each sea coastal area), it was not possible, across all our study sites, to determine whether the differences observed in elemental and isotopic ratios between the three bivalve species and across the three sea coastal areas were due to species/genera or spatial differences. However, this could be assessed at the scale of the Atlantic coastal area (see Sects. 2.4 and 3.1).

3.1 Taxonomic patterns

Different species may exhibit different C and N elemental and isotopic ratios under the same environmental and growth conditions (Mele et al., 2023), and within a single species, varying environmental conditions can lead to differences in these values (e.g., Magni et al., 2013; Briant et al., 2018). Overall, bivalve δ¹³C ranged from −23.29 ‰ to −16.98 ‰, δ¹⁵N from 2.47 ‰ to 13.15 ‰, C:N from 4.13 to 8.12 mol mol⁻¹, C from 29.81 % to 48.64 %, and N from 5.81 % to 10.57 % (out of extreme data, see Sect. 2.3). There was a significant difference in δ¹³C and C content (p<0.001) between Mytilus galloprovincialis (δ¹³C: $-\mathrm{20.05}\pm \mathrm{0.98}$ ‰, C: 36.10±1.53 %) and the other two species, Mytilus edulis (δ¹³C: $-\mathrm{19.38}\pm \mathrm{1.04}$ ‰, C: 39.80±1.93 %) and Magallana gigas (δ¹³C: $-\mathrm{19.25}\pm \mathrm{1.06}$ ‰, C: 39.30±1.90 %), which were not significantly different from each other (Fig. 2). There was a significant difference in N content (p<0.001) between Magallana gigas (8.02±0.81 %) and the two other species, Mytilus edulis (9.16±0.46 %) and Mytilus galloprovincialis (9.06±0.44 %; Fig. 2). Finally, there was a significant difference in δ¹⁵N (p<0.001) between the three species of bivalves (Magallana gigas: 8.42±1.09 ‰, Mytilus edulis 9.41±1.19 ‰, Mytilus galloprovincialis 5.70±1.49 ‰) as well as in C:N (p<0.001, Magallana gigas: 5.77±0.62 mol mol⁻¹, Mytilus edulis 5.07±0.31 mol mol⁻¹, Mytilus galloprovincialis 4.65±0.23 mol mol⁻¹; Fig. 2). Since the mussel Mytilus galloprovincialis is exclusively present in the Mediterranean Sea and the two other species, Mytilus edulis and Magallana gigas are mostly present in the English Channel and Atlantic Ocean, respectively, it is difficult to separate the species factor from the spatial factor (see Sect. 3.2) in elemental and isotopic ratios (especially in δ¹³C, where Mytilus galloprovincialis differs significantly). However, within the Atlantic coastal area, where mostly oysters are sampled ($\mathrm{15}/\mathrm{17}$ stations), two stations are sampled for mussels (stations 15-PEBE and 16-CHEM), hence species effect can be tested for this coastal area. There were no significant differences in δ¹³C between the two species (Magallana gigas: $-\mathrm{19.24}\pm \mathrm{1.07}$ ‰, Mytilus edulis: $-\mathrm{19.36}\pm \mathrm{1.02}$ ‰), but δ¹⁵N, C and N content, and C:N ratio differed significantly (p<0.001) between Magallana gigas (δ¹⁵N: 8.43±1.12 ‰, C: 39.40±1.93 %, N: 8.06±0.82 %, C:N: 5.76±0.62 mol mol⁻¹) and Mytilus edulis (δ¹⁵N: 10.54±0.88 ‰, C: 38.20±1.33 %, N: 9.34±0.39 %, C:N: 4.78±0.19 mol mol⁻¹). Note that the unequal sample sizes (15 stations for oysters, 2 stations for mussels) can reduce the power of the test and the ability to detect real differences.

3.2 Spatial patterns

3.2.1 Among sea coastal areas

The three sea coastal areas showed significant differences in bivalve δ¹³C (p<0.001). The Mediterranean Sea had lower values ($-\mathrm{20.05}\pm \mathrm{0.98}$ ‰) than the other two sea coastal areas, the English Channel ($-\mathrm{19.39}\pm \mathrm{1.04}$ ‰) and the Atlantic Ocean ($-\mathrm{19.26}\pm \mathrm{1.06}$ ‰), which were not significantly different from each other (Fig. 2). There was also a significant difference in bivalve δ¹⁵N among the three sea coastal areas (p<0.01), showing a latitudinal gradient with high δ¹⁵N values in the English Channel (9.02±1.07 ‰), intermediate values in the Atlantic Ocean (8.71±1.30 ‰), and low values in the Mediterranean Sea (5.70±1.49 ‰; Fig. 2). Bivalve C content and the C:N ratio were also significantly different among the three sea coastal areas (p<0.001), with high values in the English Channel (C: 40.0±1.89 %, C:N: 5.24±0.42 mol mol⁻¹) and Atlantic Ocean (C: 39.20±1.90 %, C:N: 5.63±0.67 mol mol⁻¹) and low values in the Mediterranean Sea (C: 36.10±1.53 %, C:N: 4.65±0.23 mol mol⁻¹; Fig. 2). Finally, N content differed significantly (p<0.001) between the Atlantic Ocean (8.23±0.89 %) and the two other sea coastal areas (English Channel: 8.94±0.69 %, Mediterranean Sea: 9.06±0.44 %), which is likely dependent on the genera/species (oysters vs. mussels; see Sect. 3.1).

3.2.2 Among stations

The lowest δ¹³C values were mostly observed in bivalves from the stations inside or close to the main river mouths: 5-CLHV ($-\mathrm{20.13}\pm \mathrm{0.70}$ ‰) and 6-VILL ($-\mathrm{20.60}\pm \mathrm{1.03}$ ‰) for the Seine River, 16-CHEM ($-\mathrm{20.15}\pm \mathrm{0.74}$ ‰) and 17-JOAL ($-\mathrm{19.90}\pm \mathrm{0.60}$ ‰) for the Loire river, 24-FOSS ($-\mathrm{20.71}\pm \mathrm{1.04}$ ‰) for the Gironde Estuary, and 30-STMM ($-\mathrm{20.34}\pm \mathrm{0.81}$ ‰) and 31-CART ($-\mathrm{20.58}\pm \mathrm{1.14}$ ‰) for the Rhône River (Fig. 2). The lowest δ¹⁵N values were mostly observed for bivalves from the Mediterranean Sea stations (all below <6 ‰), with minimum average values for station 28-BANY (4.27±0.94 ‰) and maximum for station 30-STMM (5.91±0.75 ‰). The highest δ¹⁵N values were observed for stations under the Loire and Seine river plumes (>9 ‰), e.g., 16-CHEM (10.90±0.90 ‰), 17-JOAL (9.79±0.43 ‰), 5-CLHV (9.41±1.08 ‰), and 6-VILL (9.78±1.26 ‰), and locally under the influence of smaller rivers (Fig. 2), e.g., 9-BDVG (10.16±0.52 ‰), 13-AULN (9.81±0.42 ‰), and 15-PEBE (10.19±0.71 ‰). This pattern is less clear for the Garonne and Dordogne rivers (Gironde Estuary), where δ¹⁵N shows intermediate values 23-PONT (7.74±0.60 ‰), 24-FOSS (8.04±0.73 ‰). Bivalve C and N content and C:N ratio did not show any clear spatial pattern between stations and seemed mostly linked with genera/species (differing between coastal areas; Fig. 2; see also Sect. 3.1).

Interestingly, for stations 16-CHEM and 17-JOAL, both located at the mouth of the Loire river and sampled for mussels and oysters, respectively, the δ¹³C signal did not differ significantly (mussels: $-\mathrm{19.95}\pm \mathrm{0.74}$ ‰, oysters: $-\mathrm{19.90}\pm \mathrm{0.60}$ ‰) over the same time period (2010–2021), despite being different species and genus. As for the comparison among the Atlantic coastal area (see Sect. 3.1), δ¹⁵N values significantly differed between the two species (mussels: 10.30±0.44 ‰, oysters: 9.79±0.43 ‰); however, the C:N ratio was not significantly different between the two species for this specific area (mussels: 4.74±0.19 ‰, oysters: 5.16±0.64 ‰). This suggests that spatial location could be the main driver of bivalve carbon isotope signal but that nitrogen is likely influenced by both global and local processes, including bivalve physiology, that compensate for species difference.

3.2.3 Key messages

Bivalves' elemental and isotopic ratios were generally lower for the stations of the Mediterranean Sea compared to those located in the Atlantic Ocean and English Channel. We found significant differences in carbon and nitrogen content and elemental and isotopic ratios between species (see Sect. 3.1), each of which is present in only one sea coastal area (with a few exceptions). Nevertheless, this difference is also linked with spatial variability. Indeed, C and N baselines (for isotope signal) and trophic status (i.e., nutrient and chlorophyll amounts) vary between each sea coastal area (oligotrophic Mediterranean Sea vs. meso- to eutrophic Atlantic Ocean and English Channel; Table 2) but also locally (i.e., proximity to river mouth). To overcome the lack of nutrients available in the Mediterranean Sea, some specific groups of phytoplankton (i.e., diazotrophs) uses ¹⁵N-depleted atmospheric nitrogen as the source of N, which decreases both the δ¹⁵N of the overall particulate organic matter (Kerhervé et al., 2001; Wannicke et al., 2010; Landrum et al., 2011; Liénart et al., 2017), and, by trophic propagation, the low δ¹⁵N signal is reflected in bivalves (Liénart et al., 2022, 2023). There is a clear difference in the δ¹⁵N baseline between the Mediterranean Sea and the other French sea coastal areas (Liénart et al., 2017). Generally, there was a clear difference in the δ¹³C signal close to river mouths, decreasing along river plumes, and in δ¹⁵N for some of the main (mostly eutrophicated) rivers. The enriched ¹⁵N signal of nutrients coming from rivers with extensive watershed agricultural activities and urban outlets is reflected in the high δ¹⁵N values of bivalves sampled close to river mouths (Fukumori et al., 2008; Thibault et al., 2020). Similarly, low δ¹³C values for bivalves sampled close to river mouths are mostly linked to the inputs of continental particulate material bearing this specific signal (Liénart et al., 2017). δ¹³C values of dissolved inorganic carbon (DIC) are lower in riverine/freshwater ecosystems compared to marine environments (Mook and Rozanski, 2000), resulting in more negative δ¹³C values in phytoplankton, which may be ingested locally by bivalves. Additionally, δ¹³C of particulate organic carbon (POC) from terrestrial sources tends to be more negative (Liénart et al., 2017) and can be consumed by bivalves. However, whether terrestrial POC is a significant food source for bivalves remains debated (Malet et al., 2008; Marín Leal et al., 2008; Dubois et al., 2014). Similarly, the δ¹⁵N signal of freshwater and marine phytoplankton differs due to the distinct nitrogen sources in these two environments. Therefore, it is not surprising that the δ¹⁵N signal in bivalves varies along a gradient related to proximity to river mouths.

3.3 Temporal patterns

3.3.1 Pluri-decadal δ¹³C dataset

Bivalve δ¹³C showed a significant decrease over the period from 1981–2021 for 82 % of the stations (27/33 stations; Fig. 3). The Mediterranean Sea exhibited a significant decrease for all stations (6/6 stations), the Atlantic Ocean for 82 % of the stations (14/17 stations), and the English Channel for 70 % of the stations (7/10 stations). The average decrease over the total period was consistent across all coastal areas, averaging $-\mathrm{0.58}\pm \mathrm{0.28}$ ‰ per decade, ranging from −0.22 ‰ to −1.21 ‰ per decade (i.e., −4.85 ‰ to −0.89 ‰ over the last 40 years) when the slopes were significant. This decrease in bivalve δ¹³C was more pronounced for the shorter time series starting from the late 1990s (17-JOAL, 2010–2021; 11-PASS and 12-PERS, 2001–2021; 33-POME, 1998–2021; and 7-OUIS, 1994–2021; Fig. 3). Overall, the slopes were more pronounced for the Mediterranean Sea and for the stations in the western Bay of Seine, i.e., under the influence of the Seine plume (stations 7 to 10). Some of the stations showed strong interannual variability in bivalve δ¹³C (e.g., 1-AMBL, 6-VILL, and 24-FOSS; Fig. 3).

Figure 3 Bivalve δ¹³C time series over the period from 1981–2021 for the 33 stations of the three French sea coastal areas (Magallana gigas, Mytilus edulis, Mytilus galloprovincialis). Black lines correspond to significant trends (Mann–Kendall tests; p value < 0.05), and the slope value (‰ per decade) appears in the upper-right corner. Shifts are represented by the vertical dashed line, with the year mentioned to its bottom left. One outlier is shown in grey (21-PALL, −27 ‰) but was not considered for statistical analyses.

Shifts in bivalve δ¹³C time series were detected for 79 % of the stations (26/33 stations), equally represented in each sea coastal area. The shifts occurred mostly around the year 1999–2000 (6/26 stations) or after, around the years 2006±1 year and 2012±1 year, with no specific pattern per coastal area. In the Mediterranean Sea, most of the shifts occurred before 2000. The earliest shift was in 1989 (31-CART), there were two in 1999, and the latest was in 2014 (15-PEBE). Usually, no shift was detected at the stations close to river mouths, i.e., the Seine (5-CLHV, 6-VILL), Loire (16-CHEM, 17-JOAL), and Gironde (24-FOSS) estuaries, and at some of the mouths of smaller rivers (e.g., 15-PEBE, Vilaine River, and 29-HRLT, Hérault River).

Figure 4 Bivalve δ¹⁵N time series over the period from 1981–2021 for the 33 stations of the three French sea coastal areas (Magallana gigas, Mytilus edulis, Mytilus galloprovincialis). Black lines correspond to significant trends (Mann–Kendall tests; p value < 0.05), and the slope value (‰ per decade) appears in the upper-right corner. Shifts are represented by the vertical dashed line, with the year mentioned to its bottom left.

Figure 5 Bivalve C:N time series over the period from 1981–2021 for the 33 stations of the three French sea coastal areas (Magallana gigas, Mytilus edulis, Mytilus galloprovincialis). Black lines correspond to significant trends (Mann–Kendall tests; p value < 0.05), and the slope value (‰ per decade) appears in the upper-right corner. Shifts are represented by the vertical dashed line, with the year mentioned to its bottom left.

3.3.2 Pluri-decadal δ¹⁵N dataset

Bivalve δ¹⁵N also showed a significant decrease over the period from 1981–2021 for 64 % of the stations (22/33 stations; Fig. 4). The Mediterranean Sea exhibited a significant decrease for 83 % of the stations (5/6 stations), the English Channel for 80 % of the stations (8/10 stations), and the Atlantic Ocean for only 53 % of the stations (9/17 stations). Only one station showed a significant increase (25-FERR, 0.34 ‰ over 39 years). The average decrease was $-\mathrm{0.44}\pm \mathrm{0.20}$ ‰ per decade, ranging from −0.18 ‰ to −0.84 ‰ per decade (i.e., −0.72 ‰ to −3.35 ‰ over the last 40 years; Fig. 4). However, this decrease was more pronounced in the English Channel ($-\mathrm{0.50}\pm \mathrm{0.21}$ ‰) in the eastern part and close to the Seine river mouth (e.g., 3-VARE, 4-ANTI, and 6-VILL) and in the Mediterranean Sea ($-\mathrm{0.54}\pm \mathrm{0.20}$ ‰) compared with the Atlantic Ocean ($-\mathrm{0.32}\pm \mathrm{0.08}$ ‰), where the lowest slopes were near the Loire river mouth and the inner two stations of the Arcachon Lagoon. The interannual variability was relatively low for δ¹⁵N, with the exception of some stations (e.g., 1-AMBL, 6-VILL, 16-CHEM, and 33-POME).

Shifts in bivalve δ¹⁵N time series were detected for 70 % of the stations (23/33 stations), mostly in the English Channel and Mediterranean Sea. The shifts occurred mostly around the year 1999–2000 (5/23 stations) and then around the years 1995–1996 (4/23 stations) and 2004 ± 1 year (4/23 stations), with shifts mostly occurring before 2000 in the Mediterranean Sea. For each sea coastal area, there were often shifts occurring the same year or in two close years for nearby stations (e.g., 15-PEBE, 2003, and 16-CHEM, 2004). The earliest shift was in 1991 (29-HRLT) and the latest in 2013 (14-GUIL).

Figure 6 Bivalve C content (%) time series over the period from 1981–2021 for the 33 stations of the three French sea coastal areas (Magallana gigas, Mytilus edulis, Mytilus galloprovincialis). Black lines correspond to significant trends (Mann–Kendall tests; p value < 0.05, calculated excluding extreme values mentioned in Sect. 2.3), and the slope value (‰ per decade) appears in the upper-right corner. Shifts are represented by the vertical dashed line with the year mentioned to its bottom left. Extreme values mentioned in Sect. 2.3 are not represented and were excluded from the calculations of the slopes and shifts.

Figure 7 Bivalve N content (%) time series over the period from 1981–2021 for the 33 stations of the three French sea coastal areas (Magallana gigas, Mytilus edulis, Mytilus galloprovincialis). Black lines correspond to significant trends (Mann–Kendall tests; p value < 0.05, calculated excluding extreme values mentioned in Sect. 2.3), and the slope value (‰ per decade) appears in the upper-right corner. Shifts are represented by the vertical dashed line with the year mentioned to its bottom left. Extreme values mentioned in Sect. 2.3 are not represented and were excluded from the calculations of the slopes and shifts.

3.3.3 Pluri-decadal C and N datasets

Bivalve C and N contents and C:N ratio showed significant trends for less than half of the stations (between 30 % and 36 %) over the period from 1981–2021 (C: 12/33 stations, N: 11/33 stations, C:N: 10/33 stations; Figs. 5–7). For C content, there was the same number of increasing and decreasing trends (both 18 %, 6/33 stations), with no spatial pattern, while the N content showed only decreasing trends (33 %, 11/33 stations). However, most of the significant trends were increasing for the C:N ratio (24 %, 8/33 stations), exclusively in the Atlantic Ocean (6/33 stations) and the English Channel (2/33 stations), and only 6 % (2/33 stations) of the stations showed decreasing trends, all in the Mediterranean Sea. The average increase in the C:N ratio was 0.27±0.37 mol mol⁻¹ per decade, ranging from 0.05 to 0.38 mol mol⁻¹ per decade (i.e., 0.20 to 1.53 mol mol⁻¹ over the last 40 years, excluding 17-JOAL with a 1.45 mol mol⁻¹ increase for 11 years; Fig. 5). The average C:N decrease in the Mediterranean Sea was $-\mathrm{0.05}\pm \mathrm{0.02}$ mol mol⁻¹ per decade. Interannual variability in bivalve C:N was either very low for the two mussel species (e.g., 7-OUIS, 16-CHEM, and 33-POME) or very large for the oysters (e.g., 12-PERS, 19-NOIR, and 22-BOYV). The decreasing trends in the C:N ratio were mostly linked to decreases in bivalve N % (6/8 stations), whereas only 3/8 stations showed a significant increase in C %. Shifts in bivalve C, N, and C:N time series were detected only for 18 % to 33 % of the stations (C: 6/33 stations, N: 9/33 stations, C:N: 11/33 stations), without any peculiar spatial pattern. Most of the shifts occurred after the year 2000 (C: 4/33 stations, N: 5/33 stations, C:N: 7/11 stations). The earliest C:N shift was in 1989 (27-COMP) and the latest in 2014 (6-VILL, 22-BOYV). There was no apparent spatial pattern in the C, N, and C:N shifts.

3.3.4 Key messages

In summary, over the last 30 to 40 years, bivalves have shown a general decrease in both δ¹³C and δ¹⁵N for most stations of the three French sea coastal areas. These trends are most likely linked with global processes occurring at large spatial scales, such as climate change. Liénart et al. (2024b) explained the decrease in δ¹³C is primarily associated with the increase in atmospheric CO₂ concentrations due to fossil-fuel burning, which generates a decrease in atmospheric CO₂ δ¹³C values occurring since the industrial revolution in 1850 (the Suess effect; Keeling, 1979; Gruber et al., 1999). This decrease in δ¹³C values propagates along marine food webs and is visible in organisms' tissues (Schloesser et al., 2009; Liénart et al., 2022, 2024b). However, other global factors, such as rising temperatures or climate indices, may also contribute to this overall decrease in δ¹³C through passive (e.g., Suess effect) or active processes related to bivalve physiology (e.g., change in fractionation; Liénart et al., 2024b). Decreasing δ¹⁵N trends are potentially linked to global factors, indirectly, through a general decrease in nutrients and particles inputs from rivers (Milliman et al., 2008; Bauer et al., 2013). Since organic particles generated in areas with high human activities bear a high δ¹⁵N signal, the decrease in such ¹⁵N-enriched inputs is reflected in coastal water and marine fauna (Connolly et al., 2013). However, such changes are strongly influenced by regional and local factors, such as regional climate, watershed activities, and damming (Milliman et al., 2008), resulting in different slope values in isotope time series close to river mouths (more pronounced for both C and N isotopes). The observed decrease in both isotopic ratios likely results from the cumulative, synergistic, or antagonistic effects of global and regional/local influences (Cabral et al., 2019). Such interacting pressures are more likely to occur in regions strongly affected by climate change, such as the Mediterranean Sea (Tuel and Eltahir, 2020). In this dataset, isotopes reveal both global (climate) effects, mostly through δ¹³C, and some local variability in the organic matter used by bivalves, as shown by the different δ¹³C (low) and δ¹⁵N (high) signals along the river-mouth-to-sea gradient. The decrease in bivalve δ¹³C and in δ¹⁵N over the past decades may indicate a decrease in nutrient and particle inputs from the rivers over the same period. Finally, as an indicator of bivalve condition and physiology (Elser et al., 2003), the increase in C:N ratio observed for most stations in the Atlantic and the English Channel can be attributed to either an increase in carbon content (i.e., increase in lipid or carbohydrate content), which is unlikely given rising temperatures and appears significant for only few of our stations, or a decrease in nitrogen content (i.e., proteins and amino acids), as observed in our data, as is more likely due to a decrease in N nutrient availability in the environment, which is consistent with the observed decrease in δ¹⁵N. The absence of trends at most of the stations and the large interannual variability in the C:N ratio may result from the complex interplay between global and local environmental effects and bivalve physiology related to this parameter.

4 Data availability

The dataset presented in this article and the metadata associated can be freely accessed on the SEANOE open-access repository (https://www.seanoe.org/, last access: 19 February 2025) using the following DOI: https://doi.org/10.17882/100583 (Liénart et al., 2024a). The data give the δ¹³C and δ¹⁵N (‰), C and N (%), and C:N ratio (mol mol⁻¹), for each station and year, with the sampling date and species name. For each station, it contains the sea coastal area, station ID, station name, ROCCH codes, latitude and longitude (geographic coordinate system WGS 84) as presented in Table 1. Carbon isotopic time series have already been published (Liénart et al., 2024b, except from station 17-JOAL) in a study focused on global effect on coastal ecosystems by large-scale anthropogenic and natural climate proxies, including the Suess effect, over the period from 1981–2021. Bivalve δ¹³C data and δ¹³C corrected from the Suess effect are available via the permanent identifier https://doi.org/10.6084/m9.figshare.24884871.v1 (Liénart et al., 2024c). The companion data from ROCCH (contaminants) are available in the Surval database (https://surval.ifremer.fr, Ifremer, 2022).

5 Conclusion and recommendation for use

By providing this dataset to the scientific community, including caveats for interpreting such data spatially and over time, we expect it to be useful for numerous ecological studies. Such data are relevant to trace nutrients origin and to set accurate baselines to study organism diets and food web structure and can be and indicator of water quality over space and time through bivalve physiology. It could also provide valuable input for developing predictive models of bivalve physiology (Emmery et al., 2011) or trophic ecology (Marín Leal et al., 2008). Multi-decadal time series allow scientists to understand coastal ecosystem responses to global change through biological and biogeochemical processes. The carbon isotope dataset is already part of a study where trends are thoughtfully interpreted taking into consideration global proxies for climate and anthropogenic changes and corrected for the Suess effect (Liénart et al., 2024b). Part of the observed δ¹³C trends was linked with the Suess effect, leading to shifts in the isotope baseline over recent decades (i.e., Suess effect for δ¹³C) rather than specific changes in ecosystem functioning. Hence, we would like to draw users' attention to the need to correct for the Suess effect before comparing biological samples collected 1 or more decades apart as recommended in literature (Dombrosky, 2020; Clark et al., 2021; Liénart et al., 2024b). Additionally, in order to take into account the low δ¹³C and high δ¹⁵N signal observed near river mouths, we advise users to compare slope values rather than absolute values when assessing temporal changes between stations. Overall, the complex interplay between global and regional/local effects needs to be considered when interpreting time series. This long-term dataset of elemental and isotopic values in suspension-feeder tissues provides insights into ecosystem dynamics and holds broader significance for advancing scientific understanding in the face of ongoing environmental challenges.