Rare earth elements in oyster shells: provenance discrimination and potential vital effects

Rare earth elements (REEs) and yttrium in seawater originate from atmospheric fallout, continental weathering, and transport from rivers, as well as hydrothermal activity. Previous studies have reported the use of REE and Y measurements in biogenic carbonates as a means to reconstruct these surface processes in ancient times. As coastal seawater REE and Y concentrations partially reflect those of nearby rivers, it may be possible to obtain a regional fingerprint of these concentrations from bivalve shells for seafood traceability and environmental monitoring studies. Here, we present a dataset of 297 measurements of REE and Y abundances by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) from two species (Crassostrea gigas and Ostrea edulis). We measured a total of 49 oyster specimens from six locations in France (Atlantic Ocean and Mediterranean Sea). Our study reports that there is no significant difference in concentrations from shell parts corresponding to winter and summer periods for both species. Moreover, interspecific vital effects are reported from specimens from both species and from the same locality. REE and Y profiles as well as t-distributed stochastic neighbour embedding processing (t-SNE; a discriminant statistical method) indicate that REE and Y measurements from C. gigas shells can be discriminated from one locality to another, but this is not the case for O. edulis, which presents very similar concentrations in all studied localities. Therefore, provenance studies using bivalve shells based on REEs and Y have to first be tested for the species. Other methods have to be investigated to be able to find the provenance of some species, such as O. edulis.

1 Introduction

Rare earth elements (REEs) are a group of 15 elements (La to Lu) with a similar electronic configuration of atoms, similar properties, and similar chemical behaviour (Elderfield, 1988), with the exceptions of Ce and Eu, which present multiple oxidation states. The main sources of REEs in seawater are atmospheric fallout (Elderfield and Greaves, 1982; De Baar et al., 1983) and riverine input through continental weathering (Goldstein et al., 1984; Frost et al., 1986), as well as hydrothermal activity (Olivarez and Owen, 1991). In addition to these various sources, the concentrations of REEs in seawater are impacted by adsorption processes of REEs to mineral surfaces and complexation with organic ligands, especially near the ocean surface (Sholkovitz et al., 1994; Schijf et al., 2015).

Reconstruction of the REE compositions of seawater is generally used to provide information on past continental weathering, tectonic activity, and water mass circulation (Greaves et al., 1991; Censi et al., 2004; Haley et al., 2005; Piper and Bau, 2013). Moreover, specific elemental ratios, such as Y ∕ Ho, have been investigated as potential provenance proxies. Indeed, although the average Y ∕ Ho value is equivalent to that of chondritic meteorites and mid-ocean-ridge basalt (MORB; Jochum et al., 1986; Taylor and McLennan, 1988), it has been shown that Y ∕ Ho fractionates in sediment particles, not only in seawater with depth, but also probably in waterbodies from watersheds (rivers and estuaries), whose composition has been modified depending on the weathered continental rocks (Bau et al., 1995; Nozaki et al., 1997; Prajith et al., 2015). The Y ∕ Ho ratios in estuaries could therefore exhibit different values according to the regional inputs related to the mineralogical variability in continental covers.

It has been demonstrated that the seawater composition in REEs and yttrium is at least partially recorded in carbonate materials, such as ooids (Li et al., 2019), brachiopod shells (Zaky et al., 2015, 2016), foraminifera tests (Osborne et al., 2017), and coral skeletons (Sholkovitz and Shen, 1995). In addition, the anthropogenic REE contamination of the Rhine (Germany) has been demonstrated by the shell composition of freshwater mussels (Merschel and Bau, 2015). The reconstruction of the REE and Y fingerprints in a coastal environment from mollusc shells could therefore be useful, not only to monitor potential contaminations from anthropic activities (Le Goff et al., 2019), but also as a provenance proxy for quality control of cultured organisms prior to human consumption (Bennion et al., 2019; Morrison et al., 2019). Indeed, public interest is rising regarding the origin of food for various reasons including decreased confidence in the quality and safety of the remote food supply (Kelly et al., 2005; Gopi et al., 2019). The geographic origin of seashells can also be of interest in archaeology, since mollusc shells can sometimes be unearthed very far from the nearest shoreline (Bardot-Cambot, 2014), in order to rebuild historic trade routes. Seafood traceability is therefore important for both modern and archaeological contexts. However, “vital effects” (compositional shifts between inorganic and biogenic carbonate due to metabolic activity; Urey et al., 1951) have been reported to alter this regional fingerprint of REEs between species (Akagi et al., 2004); therefore, a feasibility study needs to be conducted on molluscs before any extended seafood traceability perspective is considered.

To this end, REE and Y measurements have been performed by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) on modern and archaeological oyster shells from various French localities. The aim of this study is to assess interspecific effects using specimens from two edible species: the flat oyster Ostrea edulis (Linnaeus, 1758), commonly found at antique sites but still found on modern shores, and the cupped oyster Crassostrea gigas (Thunberg, 1793), which first appeared on the French coastline during the 20th century. Multiple measurements were performed on each of the specimens to evaluate the intraspecific variations, in particular regarding the potential seasonal fluctuations of rainfall intensity, weathering, and rate of river flows. A recent statistical method, t-SNE (t-distributed stochastic neighbour embedding; van der Maaten and Hinton, 2008), is introduced as an attempt to discriminate the nine oyster groups of the study in order to highlight significant interspecific and inter-regional differences between REE incorporations.

2 Material and methods

2.1 Modern-day settings and specimens

2.1.1 Baie des Veys (Normandy)

The Géfosse area, in Baie des Veys (Normandy, France), is currently used as a commercial oyster farm location. This open-sea area is characterized by a semidiurnal tidal range of 8 m from the English Channel and an overall siltation due to weaker ebb than flow, inducing a poor resuspension of sediment particles (Le Gall, 1970). The Baie des Veys is recharged with fresh water by the Isigny Channel, formed by the Vire and Aure rivers, and the Carentan Channel, constituted by the Douve and the Taute rivers, which drain a large part of the Bessin plain and the Cotentin. The watershed comprises limestone as well as basalt, acidic, and alkaline metavolcanic rock and diorite (Baize et al., 1997). The respective flows of the Isigny and Carentan channels are 19 and 33 m³ s⁻¹, with these rates having no significant impact on the salinity of the shoreline (Sylvand, 1995).

Oyster specimens from this locality were gathered during a previous rearing experiment conducted between 2005 and 2006 (Lartaud et al., 2010a, b; Mouchi et al., 2013). Although these specimens have been transplanted in several localities during their lives (for details, see Lartaud et al., 2010a), the part of the shells analysed in the present study is restricted to that corresponding to their Baie des Veys stay. This period is recognized in the shells owing to in vivo chemical labelling performed during the rearing experiment (Huyghe et al., 2019). Both Crassostrea gigas (n=5) and Ostrea edulis (n=5) specimens from this locality having shared the same breeding location were considered (Fig. 1, Table 1, Fig. S1 in the Supplement). Crassostrea gigas was recently renamed Magallana gigas by Salvi and Mariottini (2016); however, as this genus change is still debated (Bayne et al., 2017) and C. gigas is the most commonly found occurrence in the literature, this paper refers to this species by its original genus name. As modern specimens, these are referred to as Mod_BDV_Cgig and Mod_BDV_Oedu groups in this paper for C. gigas and O. edulis, respectively.

Figure 1 Map of the localities of modern (squares) and archaeological (stars) specimens. The coordinates of the modern sites are indicated in Table 1.

2.1.2 Marennes–Oléron bay (Charente-Maritime)

The Marennes–Oléron bay is located on the Atlantic coast of France. It is bordered by the Oléron island on the western side and opens to the Atlantic on its northern and southern borders. Its surface area is approximately 180 km². Atlantic seawater runs through the bay from north to south (Dechambenoy et al., 1977). In the area, the tidal range is 5 m with a semidiurnal rhythm. The large watershed covers 10 000 km² of land and comprises Cenozoic river deposits, limestones, and clays (Bourgueil et al., 1968, 1972; Platel et al., 1977, 1978; Bambier et al., 1982; Hanztpergue et al., 1984; Mourier et al., 1989). The bay is recharged with fresh water by two rivers, the Charente (north side), with a 36 m³ s⁻¹ flow, and the Seudre (south side) whose flow is 30 times less significant (Soletchnik et al., 1998).

The bay hosts 30 km² of aquaculture domains, including oysters. This site was also used during a rearing experiment (Lartaud et al., 2010a, b): both the C. gigas (n=3) and O. edulis (n=3) shell parts analysed in the present study refer to the shell portions corresponding to their stay at this site, highlighted by in vivo labels. These two groups are referred to as Mod_MO_Cgig and Mod_MO_Oedu for C. gigas and O. edulis specimens, respectively.

Table 1 Specimen groups and information on their respective localities. Temperature and salinity ranges at Marennes–Oléron, Tès, and Baie des Veys are from Lartaud et al. (2010b), and those at Leucate are from Andrisoa (2019).

2.1.3 Tès (Arcachon basin)

The Arcachon basin is a lagoon of 156 km² on the French Atlantic coastline. The area is subdivided into a subtidal zone and an intertidal zone with a semidiurnal tidal range of 3 m, where the studied oysters grew. Fresh water is provided from a watershed of 4138 km² by three main channels, the Eyre, the Porge, and the Landes, as well as 26 other streams and local groundwater, for a total supply of $\mathrm{1}\times \mathrm{34}\times {\mathrm{10}}^{\mathrm{6}}$ m³ of fresh water per year (Lamour and Balades, 1979; Auby et al., 1994). The largest part of this watershed consists of Cenozoic river deposits with, to a lesser extent, limestone, clay (Dubreuilh and Bouchet, 1992), and some iron oxide deposits historically used as building material (Gourdon-Platel and Maurin, 2004).

Crassostrea gigas specimens (n=8; Fig. 1, Table 1) from this locality originate from the same rearing experiment as those that were placed at Baie des Veys, and measurements were restricted to the parts of the shells corresponding to the period spent at this locality (Lartaud et al., 2010a). These specimens are referred to as Mod_TES_Cgig.

2.1.4 Leucate (Aude)

The Salses–Leucate lagoon is located on the southwestern French Mediterranean coast. It corresponds to a shallow coastal basin 14 km long and 5 km wide, separated from the Mediterranean Sea by a sandy barrier interrupted by three narrow marine inlets. The average water depth is 1.7 m, and the hydrology is a balance between the entrance of marine waters from the Mediterranean Sea, a supply of groundwater discharge from two main karstic springs with flows of 3×10⁵ and 2×10⁵ m³ d⁻¹, respectively (Fleury et al., 2007), and rainfall of approximately 500 mm yr⁻¹ restricted to the fall and spring periods. The superficial watershed covers 162 km², but the total area including the karstic waters is not yet known accurately, but it is likely extended to 60 km inland far from the shore (Salvayre, 1989; Ladagnous and Le Bec, 1997), with karstic waters penetrating Jurassic and Cretaceaous limestone and dolomite. While the tidal range is restricted, the seawater level changes in the lagoon are controlled by strong northwesterly winds regularly exceeding 10 m s⁻¹ (Rodellas et al., 2018).

Crassostrea gigas oysters (n=5; Fig. 1, Table 1) originate from a wild brood stock in the vicinity of the local oyster farming area. We were not able to collect reliable O. edulis specimens from the Mediterranean Sea shoreline for comparison, as only aquaculture specimens of C. gigas are now available here. Specimens from this group are referred to as Mod_LEU_Cgig.

2.2 Archaeological sites and specimens

2.2.1 Lyon, Auvergne–Rhône–Alpes

In the area of the Fourvière hill in Lyon, where the remains of a building were tentatively identified as a sanctuary of the goddess Cybèle, a pit was filled with food waste, which included around 200 valves of the flat oyster O. edulis (Bardot-Cambot, 2013). Absolute dating of this pit is currently being re-evaluated and is approximated to the beginning of the current era or during the 1st century CE. The provenance of these oysters is debated (Bardot-Cambot, 2013). Two groups of animals, one originating from the Mediterranean Sea coastline and the other from the Atlantic coastline, were identified based on morphometric measurements and associated mollusc shells (Bardot-Cambot, 2013). Six O. edulis specimens were selected from the first group (later referred to as the Anc_CYB1_Oedu group) and seven more from the second group (Anc_CYB2_Oedu group) for the preservation quality of their umbo (Fig. 1, Table 1).

2.2.2 La Malène, Occitanie

This medieval site (circa the 6th century CE) is located on top of a cliff and is comprised of the remains of a fortified construction (Schneider and Clément, 2012). This castrum corresponds to one of the last antique sites where oysters were consumed by the elite, with supposedly spatially restricted commercial travel (Bardot-Cambot and Forest, 2014). Still, the shells found in a dump, along with other evidence of the high social status of the occupants (such as golden currency and silver nails; Schneider and Clément, 2012), had been transported over 120 km from the Mediterranean Sea. This origin is certified because some valves are fixed on valves of Flexopecten glaber, which is endemic to the Mediterranean Sea. Six O. edulis specimens were selected for the preservation quality of their umbo (Fig. 1, Table 1). As ancient specimens, those are referred to as Anc_MAL_Oedu.

2.3 Sample preparation

All specimens were mechanically cleaned of any epibiont and were selected according to the preservation state of their umbo region. The umbo was cut from the rest of the shell and embedded in Huntsman Araldite 2020 epoxy resin. Longitudinal thick sections (approx. 750 µm thick) were manufactured to expose the preserved internal structures (Fig. 2) in order to perform geochemical analyses on this protected region away from shell surface organic or chemical contaminants. An extensive chemical cleaning of the section surfaces, as advised by Zaky et al. (2015), was not performed as the analytical surface was preserved from any external contaminants over the history of the shell. However, the influence of organic matter occluded in the crystal lattice cannot be discarded, as in any LA-ICP-MS work on biominerals.

Figure 2 Typical archaeological Ostrea edulis (a) and modern Crassostrea gigas (b) specimens. The umbo region (c) is cut following the dashed white line; laser ablation craters (200 µm in diameter) are indicated by circles. Multiple measurements have been performed on each shell for both winter and summer parts. Red and blue circles represent measurements corresponding to summer and winter periods, respectively, based on cathodoluminescence.

All sections were observed under cathodoluminescence using a cathodyne OPEA cold cathode at ISTeP, Sorbonne Université (Paris, France). Observation settings were 15–20 kV and 200–400 µA mm⁻² at a pressure of 0.05 Torr. Areas potentially affected by diagenesis or damaged were identified in order to avoid any analysis of these regions by LA-ICP-MS. In addition, cathodoluminescence (CL) observations were used to define the seasonal calibration of the umbo, according to Langlet et al. (2006) and Lartaud et al. (2010a). Only Leucate specimens did not have seasonal calibration as the CL signal was uniform and nearly absent for these specimens. A second seasonal calibration method from Kirby et al. (1998), based on a sclerochronological record on the ligamental area in the form of external convex and concave bands, was attempted. Unfortunately, umbos from Leucate specimens did not exhibit the necessary curved surface to conduct such a study. Consequently, measurement data from Leucate specimens were removed from the dataset used for the study of seasonal contrasts in the REE and Y fingerprints.

2.4 Geochemical analyses

Chemical analyses were carried out by LA-ICP-MS at ISTE (University of Lausanne, Switzerland). Measurements were performed using an Element XR (Thermo Scientific) ICP-MS coupled with a RESOlution 193 nm ArF excimer ablation system equipped with an S155 two-volume ablation cell (Australian Scientific Instruments). A pulse repetition rate of 20 Hz and an on-sample energy density of 4 J cm⁻² were used. Pre-ablation of spots was first conducted in order to clean the surface of potential contaminants that could possibly be introduced during the sanding and polishing of the samples. The analytical spots were 200 µm in diameter. Ablation was performed on the areas of each sample section corresponding to winter and summer periods (according to the cathodoluminescence seasonal calibration). This protocol allows the REE and Y incorporation to be compared for different seasons throughout the life of the oysters (Fig. 2c), with the exception of Leucate specimens, which did not exhibit seasonal cathodoluminescence signals. Sections from Leucate specimens were analysed at random positions over the umbo region instead. Multiple measurements were performed on each section to avoid bias from potential internal variability. Repeated measurements of NIST SRM 612, prior to and after each 15-sample analytical series, were used for external standardization. Accuracy was checked against measurements of the BCR-2 basalt reference material from the USGS according to the GeoReM preferred values (Jochum et al., 2005). The relative standard deviation from 22 measurements of BCR-2 was always better than 2.8 % for all REEs and Y. Measured and expected values are indicated in Table S1. Measured elements were La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu, Hf, Y, and Ca as the internal standard. Data reduction was performed using the LAMTRACE software (Jackson, 2008). A total of 297 measurements were made.

2.5 Data processing

Data processing was conducted using the MATLAB software (MathWorks; https://www.mathworks.com, last access: 17 April 2020; v. R2017a). None of the measured elements exhibit a normal distribution (Kolmogorov–Smirnov test): the distribution is right-skewed at larger element abundances. To facilitate further statistical treatment of such data, they need to be transformed to normality, for which several mathematical transforms can be used. Here, we use the cubic root transform (Chen and Deo, 2004). Seasonal differences (from the seasonal age models from cathodoluminescence) were estimated by hierarchical cluster analyses of 30 measurements from C. gigas specimens (n=5) and 41 measurements from O. edulis specimens (n=5) from Baie des Veys, as well as 24 measurements from C. gigas specimens (n=3) and 13 measurements from O. edulis specimens (n=3) from Marennes–Oléron (Fig. S2). For both species, two methods for calculating cluster distances were tested: (i) unweighted average distance and (ii) the Ward inner squared distance. A cophenetic correlation coefficient was calculated, as it is the linear correlation coefficient between the distances obtained from the cluster tree and the original distances (in the multivariate space). This coefficient is an indicator of the accuracy of the distances (estimated on the tree) to faithfully represent the dissimilarities among the observations. Multivariate analysis of variance (MANOVA) was used to compare the REE and Y fingerprints obtained from measurements of C. gigas and O. edulis from Baie des Veys against the null hypothesis that both datasets belong to the same population. Kruskal–Wallis tests were performed to compare the Y ∕ Ho ratios of multiple groups against the null hypothesis that all groups belong to the same population. A recent statistical method, t-SNE (t-distributed stochastic neighbour embedding; van der Maaten and Hinton, 2008), was used to compare and classify the multivariate dataset (exact Euclidean method). The idea of the t-SNE method is to embed high-dimensional data points in low dimensions in a way that respects similarities between points. Nearby data points in the high-dimensional space correspond to nearby embedded low-dimensional points, and distant points in high-dimensional space correspond to distant embedded low-dimensional points (MathWorks; http://www.mathworks.com, last access: 17 April 2020; v. R2017a).

3 Results

3.1 Comparisons of the interspecific and intraspecific seasonal record

Firstly, heavy REEs (Tm, Yb, Lu) and Hf were usually not detected (below 0.1 ng g⁻¹) in all the specimen groups and were therefore removed from the dataset. Secondly, measurements from the Baie des Veys and Marennes–Oléron specimens allow for interspecific and intraspecific comparisons, as they were performed on specimens from both C. gigas and O. edulis species. Data collected from the different seasons for each species did not show any significant difference in the incorporation of REEs and Y between winter and summer. For Mod_BDV_Oedu, records from both winter (n=25) and summer (n=16) in O. edulis shell samples are mixed together without clustering samples with respect to seasons (Fig. S2); in addition, the cophenetic correlation coefficients are 0.92 and 0.76 for the average and Ward methods, respectively, which emphasizes the quality of the classification. Results are similar for Mod_BDV_Cgig shells (n=19 and n=11 for winter and summer, respectively; Fig. S2), with cophenetic correlation coefficients of 0.76 and 0.69 for the average and Ward methods, respectively. Equivalent results are found for Mod_MO_Cgig (cophenetic correlation coefficients of 0.75 and 0.74 for the average and Ward methods, respectively) and Mod_MO_Oedu (cophenetic correlation coefficients of 0.94 and 0.87 for the average and Ward methods, respectively; Fig. S2). However, a comparison between C. gigas and O. edulis (interspecific comparison) clearly exhibits significant differences for both the winter (MANOVA; p value $=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$) and summer (MANOVA; p value $=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$) periods between the two species. Therefore, C. gigas and O. edulis seem to differently record their respective seasonal signals.

3.2 The Y ∕ Ho ratio as a provenance proxy

The Y ∕ Ho ratio is commonly used as a provenance proxy (Bau et al., 1995; Prajith et al., 2015). The Y ∕ Ho ratios in the present study (Fig. 3) do not display significant differences between the four localities for modern C. gigas specimens (Kruskal–Wallis; p values =0.70 between groups Mod_TES_Cgig and Mod_LEU_Cgig, 0.95 between groups Mod_TES_Cgig and Mod_BDV_Cgig, 0.98 between groups Mod_TES_Cgig and Mod_MO_Cgig, 1.00 between groups Mod_LEU_Cgig and Mod_BDV_Cgig, 0.31 between groups Mod_LEU_Cgig and Mod_MO_Cgig, and 0.57 between groups Mod_BDV_Cgig and Mod_MO_Cgig). Also, all O. edulis modern and archaeological specimens from the other localities, except Anc_CYB2_Oedu, are similar to each other (Kruskal–Wallis; p values =0.98 between groups Anc_CYB1_Oedu and Anc_MAL_Oedu, 0.99 between groups Anc_CYB1_Oedu and Mod_BDV_Oedu, 1.00 between groups Anc_CYB1_Oedu and Mod_MO_Oedu, 1.00 between groups Anc_MAL_Oedu and Mod_BDV_Oedu, 0.97 between groups Anc_MAL_Oedu and Mod_MO_Oedu, and 0.97 between groups Mod_BDV_Oedu and Mod_MO_Oedu) but are different from C. gigas shells (Table S2). Moreover, a significant difference between modern C. gigas and O. edulis shells from the same localities (Baie des Veys and Marennes–Oléron) is also reported (Table S2). However, the Anc_CYB2_Oedu group does not share the homogeneity of other modern and ancient O. edulis populations. On the contrary, the ratios measured in O. edulis specimens from this group are not significantly different from the ones obtained in modern C. gigas specimens (Kruskal–Wallis; p values =0.97, 0.26, 0.48, and 1.00 when compared with specimens from Mod_TES_Cgig, Mod_LEU_Cgig, Mod_BDV_Cgig, and Mod_MO_Cgig, respectively).

Figure 3 Box plots of Y ∕ Ho ratios for all groups. The letters on top of the boxes (a and b) identify the significant differences between groups from Kruskal–Wallis tests. Note that for the same locality at Baie des Veys (BDV) and Marennes–Oléron (MO), the Y ∕ Ho ratios are significantly different depending on the species considered (C. gigas and O. edulis groups). Grey bars represent median values; the lower and higher large black bars represent the 25th and 75th percentiles, respectively, and the lower and higher small black bars represent the minimum and maximum values not considered outliers, respectively. Outliers are represented by grey crosses.

3.3 REE incorporation in shells

For all the specimens, a gradual decrease in REE dispersion (i.e. amplitude of values relative to Y) is generally observed with increasing atomic number (Fig. 4). Indeed, several groups can be identified with light REEs (e.g. Pr and Sm), such as the C. gigas groups. In contrast, middle REE (e.g. Dy and Er) distributions appear similar for all groups, with the exception of the Mod_LEU_Cgig group, whose measurements are isolated from those of the other groups. The REE median profiles (normalized to Post-Archean Australian Shale; McLennan, 1989) also present similar trends of middle REEs (from Gd to Er) for most groups (Fig. 5), except for the Mod_LEU_Cgig group, which exhibits lower abundances than the other groups for all REEs. However, light REEs are generally substantially depleted in C. gigas specimens compared to modern and ancient O. edulis groups (approx. 1 order of magnitude difference). The only exceptions are La and Ce in C. gigas shells from Baie des Veys (Mod_BDV_Cgig), which present values in the range of those from modern and ancient O. edulis specimens. This pattern of enriched Ce (and to a lesser extent, La) is not shared by O. edulis specimens from this same locality (Mod_BDV_Oedu). Although these two elements have similar abundances for both species in this locality, all the other REE abundances are different.

Figure 4 Gradual decrease in REE abundances in oyster shells according to the REE atomic number, presented against Y. Values are expressed in the cubic root of abundances (µg g⁻¹) to approach normality. Measurements from C. gigas and O. edulis are indicated by crosses and filled circles, respectively.

Figure 5 REE median profiles for all groups of oyster specimens. Crassostrea gigas groups are symbolized by dashed lines, and continuous lines are for O. edulis groups. Seawater profiles of the Atlantic Ocean (van der Flierdt et al., 2012) and the Mediterranean Sea (Censi et al., 2004) are indicated for comparison. Values are normalized to Post-Archean Australian Shale (PAAS) according to McLennan (1989).

Results from the entire dataset (i.e. 297 measurements and 12 elements per measurement) are compared in Fig. 6 using t-SNE. Some groups are well identified by this method, such as the Mod_LEU_Cgig and Mod_BDV_Cgig groups (Mediterranean Sea and British Channel coastlines). Both Mod_MO_Cgig and Mod_TES_Cgig groups, originating from the Atlantic Ocean coastline, appear as one unique group by t-SNE, with a restricted dispersion. The modern and ancient O. edulis groups are, however, not discriminated by t-SNE and present a substantially larger dispersion than the Atlantic Ocean C. gigas groups. In this sample set, the Anc_CYB1_Oedu, Mod_MO_Oedu, and Mod_BDV_Oedu groups are relatively similar in terms of the range of distribution, and Anc_CYB2_Oedu and Anc_MAL_Oedu specimens share similarities, but the five groups remain poorly differentiated.

Figure 6 Visualization of shell group partitioning using t-SNE applied to all REE and Y measurements as variables. Crosses and large dots refer to C. gigas and O. edulis, respectively.

4 Discussion

In both studied species, the decrease in the range of variation of REE abundances with increasing atomic number (except for Tm, Yb, Lu, and Hf, which were not quantified in the shells) can be explained by the increased affinity of heavy REEs for complexation in seawater, as has been demonstrated in previous studies (Cantrell and Byrne, 1987; Byrne and Kim, 1990; De Baar et al., 1991). As these elements are trapped in complexed forms or ligands, their bioavailability in seawater is strongly reduced, limiting their insertion in the oyster ionic pumps leading to the mineralization locus. These bioavailability restrictions of REEs have already been demonstrated in the freshwater mussel Corbicula fluminea for Gd (Merschel and Bau, 2015) and heavy REEs in Mytilus edulis (Ponnurangam et al., 2016). Another explanation can be advanced regarding the technique used. An LA-ICP-MS device analyses both the mineral and organic phases ablated from the biomineral without the possibility to assess their relative proportions. Although the mean proportion of organic compounds in oyster shells is limited (<0.5 % for C. gigas; Mouchi et al., 2016), it is known that organic REE abundances are depleted in heavy REEs (Freslon et al., 2014). The decreasing abundance with increasing atomic number may then be caused by protein and polysaccharide contents. Only extensive cleaning for solution-based ICP-MS analyses would be able to entirely remove the organic molecules before measurements, but this would not fit the fast REE assessment of a large number of specimens we aimed to conduct in this study.

In this study, the Y ∕ Ho ratios, which are usually proposed as a provenance proxy (Bau et al., 1995; Prajith et al., 2015), are affected by strong vital effects for both oyster species, potentially due to the decrease in REE abundance with increasing atomic number. In addition, significant differences between species from the same locality are also reported. Hence, the Y ∕ Ho ratio generally does not depend on the original location (Fig. 3). Consequently, Y ∕ Ho should not be used directly as a provenance proxy (at least from LA-ICP-MS data collected on biogenic carbonates), or it should be used with extreme caution after having discarded any potential vital effect. As the decreased REE abundance with increasing atomic number discussed above prevents a locality-specific variation of Ho, other Y ∕ REE ratios have been tested as alternative provenance proxies using lighter REEs (Fig. 4). Y ∕ La, Y ∕ Ce, Y ∕ Pr, and Y ∕ Nd were all unsuccessful in providing the identification of locality groups and also present similar values for all modern and ancient O. edulis. For Y ∕ La and Y ∕ Ce ratios, C. gigas specimens from Baie des Veys were identical to all modern and ancient O. edulis specimens, while for Y ∕ Pr and Y ∕ Nd ratios, C. gigas specimens from Baie des Veys were identical to all C. gigas specimens from the other localities. Overall, Y ∕ REE ratios appear unsuccessful for provenance discrimination.

Measurements performed on Baie des Veys and Marennes–Oléron specimens (Mod_BDV_Cgig and Mod_BDV_Oedu, as well as Mod_MO_Cgig and Mod_MO_Oedu) have been used to study the incorporation of REEs inside a single species in order to evaluate its intraspecific variation. Contrary to the model experiments performed for different temperatures on the mussel Mytilus edulis (Ponnurangam et al., 2016), the seasonal conditions do not have any impact on the REE incorporation for O. edulis or C. gigas shells (in the range of 5–20 ^∘C; Table 1). Parameters other than temperature and pH, used by Ponnurangam et al. (2016), are probably in effect, which lowers the impact of temperature on REE incorporation. This observation implies that any part of a shell can be sampled without necessarily having to define a temporal calibration of the umbo. However, REE abundances fluctuate widely within a single specimen, and the “local fingerprint” of these elements needs to be clarified and based on multiple measurements performed on each of several specimens. These intra-individual fluctuations cannot be due to seasonally controlled environmental factors, such as temperature, precipitation, or plankton blooms. However, a potential source of REEs for oysters can be pore water or resuspended sediment (Haley et al., 2004; Crocket et al., 2018), and therefore the relative abundances may fluctuate without precise temporal cyclicity.

The reasons for these intraspecific and interspecific vital effects remain unknown. Indeed, as far as we know, no study has ever shown evidence or suspicion of the use of REEs in metabolic processes that could induce an effective filter of these elements between seawater and the extrapallial cavity where shell mineralization occurs. Nevertheless, it has been reported that REEs, or other unsuspectedly useful elements, are indeed used by organisms in specific environmental settings; an example is provided by diatoms in Zn-depleted conditions, under which Zn is used as a co-factor of carbonic anhydrase (Lee et al., 1995). Another example is given by methanotrophic archaea, which use Cd, a toxic element, as a co-factor of methanol dehydrogenase (Pol et al., 2014).

The incorporation of REEs differs between the two studied species. C. gigas shells not only present different REE profiles between groups (unlike those of O. edulis), but the positive Gd anomaly, a characteristic of modern coasts under pressure from anthropic activities (Bau and Dulski, 1996; Nozaki et al., 2000; Le Goff et al., 2019), is also observed solely for “modern” C. gigas specimens from the Mod_BDV_Cgig, Mod_MO_Cgig, and Mod_TES_Cgig groups (Fig. 5), which correspond to the only watersheds with major cities. The Gd anomaly is visible neither in modern C. gigas from Mod_LEU_Cgig nor in modern O. edulis. The systematically low abundances of REEs in the Leucate shells can also be explained by the regional geology (Ladagnous and Le Bec, 1997), as watersheds in the other localities of C. gigas specimens present substratum types with higher REE contents than karsts (e.g. basalts; Baize et al., 1997). The dispersion of measurements using t-SNE seems to be ineffective in discriminating between the Mod_MO_Cgig and Mod_TES_Cgig groups. Both these groups are located on the Atlantic coastline of France, and both respective watersheds comprise the same rock types (i.e. mainly limestones and some clays and river deposits). One could therefore expect a similar riverine water REE content. Overall, these species-specific characteristics indicate that C. gigas can be used as a sentinel species regarding REE pollution of coastal waters. However, we cannot confirm that O. edulis is a proper candidate for such studies.

Several reasons for these interspecific differences can be advanced. It is known that oysters can be selective in their diet, composed mainly of diatoms (Yonge, 1928; Paulmier, 1971) of a specific size range, and preferentially digest specific species of diatoms over others (Shumway et al., 1985; Cognie et al., 2001). If food is a source of REEs, it may be possible that each oyster species does not feed on the same prey, which can present different abundances of these elements. Aquarium experiments have reported different ingestion rates in the 5–15 µm algal size range between these oyster species (Nielsen et al., 2017), but there is no indication of the REE content of the food. Alternatively, O. edulis, which exhibits generally higher abundances of REEs (nearly 1 order of magnitude higher for light REEs, except for La and Ce for Mod_BDV_Cgig specimens; Fig. 5), could present a higher bioaccumulation of these elements in its soft tissues (and eventually shell) compared to C. gigas. Ong et al. (2013) presented trace element measurements from the soft tissues of both species from the Baie de Quiberon (Brittany, France), indicating that the soft tissues of O. edulis contain generally less Cu and Zn but more Cd and Pb than those of C. gigas. Such species-specific bioaccumulation and incorporation differences could also be in effect for REEs. Finally, it is possible to explain the higher abundance of light REEs in O. edulis shells compared to C. gigas by suggesting that O. edulis ingests more clay particles. As heavy REEs are trapped in complexed form in seawater, mainly light REEs must be available for adsorption on clay particles and are eventually integrated into the forming carbonate shell.

In any case, this study shows that t-SNE can be used on REE and Y measurements from C. gigas shells to identify the regions of origin of specimens from this species. However, it appears that intraspecific vital effects prevent its efficiency in O. edulis, whose specimens exhibit the same fingerprint for several localities of origin. For this reason, we cannot confirm or refute a different origin of the two populations of Cybèle archaeological specimens (Anc_CYB1_Oedu and Anc_CYB2_Oedu) with these elements.

5 Conclusions

Multiple types of vital effects and other factors regarding REE incorporation in C. gigas and O. edulis oyster shells have been highlighted in this study. Intraspecific variations in REE abundances are significant but not related to seasonal fluctuations. A gradual decrease in REE incorporations with increasing atomic numbers has been observed, and it appears that heavy REEs are less discriminant than light REEs in identifying the various studied groups. The Y ∕ Ho ratio, previously reported as a proxy for provenance studies, is ineffective in oyster shells. Finally, interspecific variations underline the ability of the t-SNE procedure to correctly separate C. gigas specimens of various regions of origin but not O. edulis specimens, which also implies that only C. gigas can be used as a monitor species of light REE pollution. Reconstruction of the provenance of oyster specimens will therefore have to be performed separately for each studied species, as the regional geochemical fingerprints of the shells appear to be species-dependant. In order to be able to identify the regions of origin of species affected by strong vital effects (such as O. edulis), it is necessary to investigate other chemical elements as potential provenance proxies. Moreover, increased efficiency in identifying the locality (and not only the region) of origin of C. gigas specimens could be achieved by also measuring elements other than REEs.