Using Mg/Ca Ratios from the Limpet Patella depressa Pennant, 1777 Measured by Laser-Induced Breakdown Spectroscopy (LIBS) to Reconstruct Paleoclimate

: Measurement of the elemental composition of shells is increasingly emerging as an avenue for obtaining high-resolution insights into paleoclimate and past seasonality. Several studies have shown signiﬁcant correlations between Mg/Ca ratios measured on shell carbonate and the sea surface temperature (SST) within which this carbonate was precipitated. However, other investigations have reported large variability in this relationship between species. Therefore, further studies, including taxa previously not considered are still required in order to validate these new species as suitable climate proxies. Here, we measured Mg/Ca ratios for limpet Patella depressa Pennant, 1777 samples live-collected in northern Spain for the ﬁrst time. The elemental ratio was measured using laser-induced breakdown spectroscopy (LIBS), a technique that signiﬁcantly decreases the time required for sample preparation and increases the number of shells that can be analyzed. In this study, calibration-free LIBS (CF-LIBS) methods were applied to estimate molar concentrations of chemical elements on biogenic calcium carbonate. The Mg/Ca ratio evolution along the shell growth axis was compared with stable oxygen isotope ( δ 18 O) proﬁles obtained from these same limpets and the SST at the place where the mollusk grew to determine if the sequences obtained correctly reﬂected environmental conditions during the life-span of the mollusk. The results showed a signiﬁcant correlation between Mg/Ca ratio series and both δ 18 O proﬁles and SST, highlighting the paleoenvironmental and archaeological potential of LIBS analyses on this mollusk species that is frequently found in archaeological contexts in the western Europe.


Introduction
There has long been a keen theoretical and methodological interest in archaeology at elucidating the impacts of climatic and environmental change on hominin behavior, settlement, and culture in different parts of the world at different points of time [1][2][3][4][5], and how this might be mobilized to inform challenges facing our species in the 21st century [6,7]. In western Europe, climate changes have been variously argued to have had deep implications for Neanderthal behavior and eventual extinction [8][9][10], the tempo and nature of human hunter-gatherer settlement across the Last Glacial Maximum and Early-Middle Holocene [11][12][13], and the arrival and adaptations of the first farmers [14,15]. Crucial to such debates, however, is the ability to develop a set of proxies that can reconstruct the impacts of broader climatic systems on local environments used by past human populations, both marine and terrestrial [16][17][18]. Furthermore, from a behavioral perspective, some of the most critical climatic changes impact intra-annual seasonality of temperature and rainfall, something that has been notoriously difficult to reconstruct for the past [19,20].
The shells of marine mollusks are routinely found in many archaeological deposits due to littoral resources that were exploited by past hominins at least from the Middle Pleistocene [21][22][23][24], but especially during Upper Paleolithic and Mesolithic periods [25][26][27]. Previous investigations have demonstrated that shell remains can offer very valuable information on human settlement and mobility patterns, as well as on subsistence strategies of these past populations [28][29][30][31]. At the same time, shells can also be considered as a particularly relevant paleoenvironmental proxy in the context of exploring past marine conditions [16,32]. Marine mollusk species have very specific environmental requirements and are very sensitive to climate change [33,34]. Warmer conditions during the Early Holocene provoked significant changes in species representation along the Atlantic coasts of Europe compared to the Pleistocene assemblages [35,36]. Meanwhile, the current rising of seawater temperatures is provoking a rapid expansion of the northern limit of warm-adapted marine species [37,38]. Beyond taxonomic studies, during the last few decades, the isotopic measurement of biogenic carbonates within mollusk shells has been found to provide a highly effective paleothermometer over the course of a mollusk's lifespan [16,20,39,40]. Stable oxygen isotope ratios derived from marine mollusk shells (δ 18 O shell ) are powerful recorders of the seasonal sea surface temperature (SST) variations experienced during the precipitation of carbonate year-round in the past [39,[41][42][43]. Not only that, but isotopic series can be used to define the season(s) when the mollusks were harvested by past hominin groups, thus deciphering seasonal exploitation patterns of this food supply during prehistory [44][45][46][47].
Although previous studies on modern marine mollusk shells have shown a high correlation between δ 18 O shell values and SST [48][49][50][51], oxygen isotope fractionation is also influenced by oxygen isotope composition of the seawater (δ 18 O sw ) [20,52]. This represents a significant limitation for paleoreconstructions where δ 18 O sw values are unknown, while the cost-and time-consuming nature of serial δ 18 O shell sampling can also be prohibitive. This has led to efforts to develop independent paleotemperature archives. Concentration ratios of some trace elements to calcium have been successfully tested on corals and foraminifera [53][54][55]. However, despite extensive research on mollusk shells, relatively few studies have shown a clear relationship between the element/Ca ratios obtained and SST [39,[56][57][58], with many reporting a low (or no) correlation between element/Ca ratios and SST in different mollusk species [59][60][61][62][63][64]. This indicates significant species-specific variability in this regard, as well as divergence in how trace elements are incorporated into carbonate in different habitats [65]. Therefore, more analyses involving new shell taxa and species collected from coastal areas previously not considered are required to better understand trace element precipitation mechanisms and the wide-scale effectiveness of this approach as a past environmental proxy.
In this study, we measured, for the first time, Mg/Ca ratios of limpet Patella depressa Pennant, 1777 samples live-collected in northern Spain. This mollusk species was consumed in high numbers by Mesolithic humans during the Early Holocene period, along the Atlantic coast of Europe, representing an important percentage of all shell taxa recovered from shell middens excavated in northern Spain and France [66,67]. Determining if Mg/Ca ratios derived from P. depressa can be used as a suitable paleothermometer to reconstruct paleoclimate conditions and seasonality would have significant implications for future archaeological and paleoenvironmental investigations. Mg/Ca ratios were measured using laser-induced breakdown spectroscopy (LIBS), one of the most novel techniques for measuring element/Ca ratios on shells due to significantly easier sample preparation and lower measurement times [57,58,[68][69][70]. A calibration-free LIBS (CF-LIBS) technique [71] was applied in order to estimate the molar concentrations of both chemical elements. This procedure was applied for the first time to marine mollusk shell biogenic carbonates. Although a recent isotopic investigation carried out on this species has shown that δ 18 O shell derived from this species correctly reflects seasonal temperature variations year-round in northern Spain [72], LIBS-derived Mg/Ca ratio could offer significant improvements in relation to δ 18 O shell as a SST recorder, allowing an increase in the number of shells measured for each stratigraphic unit, as well as the wider accuracy of archaeological conclusions. In order to determine if elemental series measured on P. depressa are a suitable paleothermometer, Mg/Ca ratio series obtained were compared with stable oxygen isotope ratio profiles previously obtained from these same four marine mollusk shells [72], as well as with measured SST information at the place where mollusks grew.

Modern Shells and Sea Surface Temperature (SST)
Four P. depressa live-collected limpets were collected at Langre Beach (Cantabria, northern Spain, Figure 1) (43 • 28 37" N, 3 • 41 44" W) on 1 October 2012, corresponding to the end of the warmer season (i.e., summer). The specimens were gathered from the rocky intertidal shore and they were sacrificed immediately after collection by immersion in boiling water for one minute, thus avoiding further deposition of calcium carbonate. To determine if the LIBS-derived Mg/Ca series from these shells are dependent on environmental conditions, they were compared with measured year-round SST variations during the life-span of the mollusks. Daily instrumental SST data were provided by the Spanish Institute of Oceanography (Santander, Cantabria) (Figure 1), which is located close to Langre beach (<10 km). The salinity conditions of the sea (35.6 ± 0.7 PSU), as well as the oxygen isotope composition of the seawater (δ 18 O sw = 0.9 ± 0.2‰ [VSMOW]) are similar in both areas, with no influence of continental runoff or sea currents with different salinity levels [73,74]. Appl. Sci. 2021, 11, 2959 3 of 17 significant implications for future archaeological and paleoenvironmental investigations. Mg/Ca ratios were measured using laser-induced breakdown spectroscopy (LIBS), one of the most novel techniques for measuring element/Ca ratios on shells due to significantly easier sample preparation and lower measurement times [57,58,[68][69][70]. A calibration-free LIBS (CF-LIBS) technique [71] was applied in order to estimate the molar concentrations of both chemical elements. This procedure was applied for the first time to marine mollusk shell biogenic carbonates. Although a recent isotopic investigation carried out on this species has shown that δ 18 Oshell derived from this species correctly reflects seasonal temperature variations year-round in northern Spain [72], LIBS-derived Mg/Ca ratio could offer significant improvements in relation to δ 18 Oshell as a SST recorder, allowing an increase in the number of shells measured for each stratigraphic unit, as well as the wider accuracy of archaeological conclusions. In order to determine if elemental series measured on P. depressa are a suitable paleothermometer, Mg/Ca ratio series obtained were compared with stable oxygen isotope ratio profiles previously obtained from these same four marine mollusk shells [72], as well as with measured SST information at the place where mollusks grew.

Modern Shells and Sea Surface Temperature (SST)
Four P. depressa live-collected limpets were collected at Langre Beach (Cantabria, northern Spain, Figure 1) (43°28'37" N, 3°41'44" W) on 1 October 2012, corresponding to the end of the warmer season (i.e., summer). The specimens were gathered from the rocky intertidal shore and they were sacrificed immediately after collection by immersion in boiling water for one minute, thus avoiding further deposition of calcium carbonate. To determine if the LIBS-derived Mg/Ca series from these shells are dependent on environmental conditions, they were compared with measured year-round SST variations during the life-span of the mollusks. Daily instrumental SST data were provided by the Spanish Institute of Oceanography (Santander, Cantabria) ( Figure 1), which is located close to Langre beach (<10 km). The salinity conditions of the sea (35.6 ± 0.7 PSU), as well as the oxygen isotope composition of the seawater (δ 18 Osw = 0.9 ± 0.2‰ [VSMOW]) are similar in both areas, with no influence of continental runoff or sea currents with different salinity levels [73,74].

LIBS Experimental Set-Up
LIBS is an atomic emission spectroscopic technique that ablates a small volume of material with a high-energy pulsed laser, creating a plasma. The light emitted from the plasma contains emission lines at different wavelengths depending on the sample composition at the elemental and molecular levels [75]. The Mg/Ca ratios reported in this work were obtained from the processing of Mg and Ca lines' emission intensities. The experimental setup (Figure 2) consisted of a Nd:YAG laser (Lotis LS-2134D) operating at a wavelength of 1064 nm, with 16 ns pulse width and 10 Hz repetition rate. The pulse energy was set to 40 mJ. The laser spot was focused with a lens of 75 mm focal length, producing a crater of 200 µm diameter and 0.2 µm depth (average) for each laser shot. The light emission from the plasma was collected with an optical fiber (800 µm diameter) coupled to a bundle of 8 fibers for each channel of an Avantes ULS2048-USB2-RM CCD spectrometer with a total wavelength range from 178 to 889 nm, and a resolution from 0.015 to 0.06 nm. The capture window to collect the plasma-emitted light was set to start 5 microseconds after the laser shot, with a length of 1 millisecond (the minimum capture time for this non-gated spectrometer). A trough-lens vision system with a CCD camera allowed the inspection of the sample surface. The sample position was controlled by a motorized XYZ positioner, programed to follow a measurement path along the shell surface, which is described in the next section.

Shell Sample Preparation and Mg/Ca Ratio Analysis
Following the methodological procedures previously utilized in investigations conducted on Patella limpet species, shells were sectioned along the growth axis before being analyzed by LIBS [57,69]. The selected specimens were partially coated with a metal epoxy resin along maximum growth axis (i.e., from the anterior to the posterior margin) to avoid the shell breaking when sectioned ( Figure 3a). Sectioning was performed using a Buehler Isomet low-speed saw and a diamond wheel at the IIIPC-University of Cantabria (Spain). Firstly, two thick sections (~3 mm each) were obtained from each limpet, which were previously used for conducting the first sclerochronological analysis carried out on P. depressa species [72]. One out of the two remaining sections was used for the elemental analysis conducted here (Figure 3a). The remaining right hand limpet sections were ground on glass plates using 600 and 800 SiC grit powder and polished with a 1 µm diamond suspension grit until the internal growth lines and increments were clearly visible. According to the methodology previously applied to several limpet species [39,57,65,69], Mg/Ca ratios were sequentially measured every 0.1 mm on the calcite layer equivalent to the m + 2 layer [72,76] from the shell edge (i.e., the last portion of the shell growth) to the shell apex (i.e., the first portion of the shell growth) (Figure 3b). The measurement path from the edge was perpendicular to the growth lines and the diameter of the LIBS spot was always lower than 0.2 mm. shell apex (i.e., the first portion of the shell growth) (Figure 3b). The measurement path from the edge was perpendicular to the growth lines and the diameter of the LIBS spot was always lower than 0.2 mm.

Elemental Spectra Processing and Mg/Ca Ratio (mmol/mol) Calculation
In the recorded spectra, due to the high resolution and wide spectral range of the eight-channel spectrometer, many emission lines can be seen. Most of them are CaI and CaII lines, with several MgI and MgII lines and a few from SrI emission, in agreement with the expected chemical composition of mollusk shell carbonate. Although the presence of lithium in the carbonate of other species has been reported, no Li lines were found in our measurements. Due to the absence of specific matrix-matched calibration standards for marine mollusk shell biogenic carbonates, previous investigations carried out on marine mollusk shell biogenic carbonates obtained the Mg/Ca values from the raw ratio of two elemental lines, expressed in arbitrary units [57,58,68,69]. In this study, a calibrationfree LIBS (CF-LIBS) technique was applied for the first time on biogenic carbonates in order to estimate molar concentrations of magnesium and calcium chemical elements. The algorithm used is an implementation of the method originally proposed by Ciucci et al. [71]. It obtains the concentrations from the Boltzmann and Saha-Boltzmann plots of many emission lines of the species of interest (in our case, CaI, CaII, MgI, and MgII) (Figure 4), and takes into account the electronic temperature and electron density of the plasma, thus reducing the variability of the signal due to shot-to-shot plasma fluctuations. To improve the accuracy of the results, the emission lines used by the CF algorithm have been carefully selected based on their signal to noise ratio in the spectra, energy levels, and selfabsorption properties. One common issue of the CF-LIBS approach, in applications that use many measurements at different spatial points to create a sequence of concentrations, is that decisions to consider or discard specific emission lines are made individually for each spectrum. These decisions are dynamically made based on the particular signal to noise ratio of each emission peak, the quality of fit to a given peak function, or whether peak parameters are outliers in the Saha-Boltzmann plots or in the calculation of electron density. To reduce the variability of the concentrations determined across all the measurements of the sequences due to changes in the used emission lines, only the set of lines that were valid for all spatial points were considered here. This results in 16 lines of CaI, 5 of CaII and 6 of MgI being considered for the Saha-Boltzmann plots. As all MgII lines were discarded due to their high self-absorption, the concentration of ionized species was extrapolated using the Saha equilibrium equation [71]. The Matlab implementation of the CF-LIBS algorithm used in this work is openly available at github.com/acobo/CFLIBS. This procedure allowed us to obtain the relative concentration in mmol/mol units and provided significant advantages to the line ratio approach, decreasing LIBS-derived

Elemental Spectra Processing and Mg/Ca Ratio (mmol/mol) Calculation
In the recorded spectra, due to the high resolution and wide spectral range of the eight-channel spectrometer, many emission lines can be seen. Most of them are CaI and CaII lines, with several MgI and MgII lines and a few from SrI emission, in agreement with the expected chemical composition of mollusk shell carbonate. Although the presence of lithium in the carbonate of other species has been reported, no Li lines were found in our measurements. Due to the absence of specific matrix-matched calibration standards for marine mollusk shell biogenic carbonates, previous investigations carried out on marine mollusk shell biogenic carbonates obtained the Mg/Ca values from the raw ratio of two elemental lines, expressed in arbitrary units [57,58,68,69]. In this study, a calibrationfree LIBS (CF-LIBS) technique was applied for the first time on biogenic carbonates in order to estimate molar concentrations of magnesium and calcium chemical elements. The algorithm used is an implementation of the method originally proposed by Ciucci et al. [71]. It obtains the concentrations from the Boltzmann and Saha-Boltzmann plots of many emission lines of the species of interest (in our case, CaI, CaII, MgI, and MgII) (Figure 4), and takes into account the electronic temperature and electron density of the plasma, thus reducing the variability of the signal due to shot-to-shot plasma fluctuations. To improve the accuracy of the results, the emission lines used by the CF algorithm have been carefully selected based on their signal to noise ratio in the spectra, energy levels, and self-absorption properties. One common issue of the CF-LIBS approach, in applications that use many measurements at different spatial points to create a sequence of concentrations, is that decisions to consider or discard specific emission lines are made individually for each spectrum. These decisions are dynamically made based on the particular signal to noise ratio of each emission peak, the quality of fit to a given peak function, or whether peak parameters are outliers in the Saha-Boltzmann plots or in the calculation of electron density. To reduce the variability of the concentrations determined across all the measurements of the sequences due to changes in the used emission lines, only the set of lines that were valid for all spatial points were considered here. This results in 16 lines of CaI, 5 of CaII and 6 of MgI being considered for the Saha-Boltzmann plots. As all MgII lines were discarded due to their high self-absorption, the concentration of ionized species was extrapolated using the Saha equilibrium equation [71]. The Matlab implementation of the CF-LIBS algorithm used in this work is openly available at github.com/acobo/CFLIBS. ment, or spectrometer configuration) between specimens, but also the variability found along a single cross-section of a limpet. In addition, to better assess the long-term (temperature related) ratio variations, the well-known singular spectrum analysis (SSA) was applied to the sequences [77,78], a methodology that has already been applied on LIBSderived series from marine mollusk shells elsewhere [57]. Each sequence was decomposed into 2 signal components and 6 noise components; the latter then being discarded. The effect of SSA is shown in Figure 5.   This procedure allowed us to obtain the relative concentration in mmol/mol units and provided significant advantages to the line ratio approach, decreasing LIBS-derived variability caused by diverse experimental conditions (room temperature, optical alignment, or spectrometer configuration) between specimens, but also the variability found along a single cross-section of a limpet. In addition, to better assess the long-term (temperature related) ratio variations, the well-known singular spectrum analysis (SSA) was applied to the sequences [77,78], a methodology that has already been applied on LIBS-derived series from marine mollusk shells elsewhere [57]. Each sequence was decomposed into 2 signal components and 6 noise components; the latter then being discarded. The effect of SSA is shown in Figure 5. This procedure allowed us to obtain the relative concentration in mmol/mol units and provided significant advantages to the line ratio approach, decreasing LIBS-derived variability caused by diverse experimental conditions (room temperature, optical alignment, or spectrometer configuration) between specimens, but also the variability found along a single cross-section of a limpet. In addition, to better assess the long-term (temperature related) ratio variations, the well-known singular spectrum analysis (SSA) was applied to the sequences [77,78], a methodology that has already been applied on LIBSderived series from marine mollusk shells elsewhere [57]. Each sequence was decomposed into 2 signal components and 6 noise components; the latter then being discarded. The effect of SSA is shown in Figure 5.

Estimation of Mg/Ca Ratio by CF-LIBS
The elemental Mg/Ca ratio was obtained on a predefined measurement path of the calcite layer for the four shells analyzed (Figure 3b). These ratio sequences are shown in Figure 6a-d. Values on the left-hand side correspond to the shell edge (date of collection) and values on the right-hand side correspond to the area close to the apex, i.e., earlier dates in the mollusk's lifespan.

Estimation of Mg/Ca Ratio by CF-LIBS
The elemental Mg/Ca ratio was obtained on a predefined measurement path of the calcite layer for the four shells analyzed (Figure 3b). These ratio sequences are shown in Figure 6a-d. Values on the left-hand side correspond to the shell edge (date of collection) and values on the right-hand side correspond to the area close to the apex, i.e., earlier dates in the mollusk's lifespan.  /mol) were relatively consistent between the four shells measured, although LAN554 and LAN559 reported a maximum and a minimum value lower than those derived from the other shells, respectively. Maximum (59.2 mmol/mol) and minimum (31.6 mmol/mol) Mg/Ca ratio values obtained herein were notably higher than elemental values previously measured on several calcite carbonate mollusk species using ICP-based techniques. Ferguson et al. [39] and Freitas et al. [56] obtained molar concentration values lower than 30 mmol/mol from modern limpets and scallops collected in the southern Iberian Peninsula. Likewise, Graniero et al. [65] and García-Escárzaga et al. [69] analyzed Patella vulgata limpets from Atlantic Europe, reporting elemental ratios (from ca. 10 to ca. 30 mmol/mol) lower than obtained herein from P. depressa species by LIBS. We attribute this higher concentration of Mg measured by LIBS techniques to a possible species-dependent variability, but also to the CF-LIBS algorithm itself. The results could be affected by the self-absorption of Ca lines, as self-absorption tends to lead to an underestimation of the concentration of major elements in the chemical composition. The current implementation of the algorithm tries to reduce this effect with a careful selection of the emission lines, but the effect is difficult to remove entirely [79].
Another noticeable aspect of the sequences was their high variability, with apparent high-frequency fluctuations on top of the expected yearly sinusoidal variations. Although much research using the LIBS technique has tended to explain general intrinsic variability in LIBS spectra as a consequence of fluctuations in the varying plasma properties of each laser shot [80], we attribute this behavior to the actual composition heterogeneity of the biogenic carbonates on a micrometric scale, as has been suggested elsewhere [58,61,69]. As other studies report, the mineralization process in limpets may not only be dependent on seawater temperature, but also on other biological and environmental factors, such as mollusk ontogeny, intra-annual variation in shell growth rates, individual physiological processes, seawater salinity, and/or food supply conditions [60,64,81,82]. To evaluate the reliability of measured Mg/Ca ratios, one limpet (LAN545) cross-section was measured twice along the same path. Although the ablated material was not strictly the same because the second series was measured on material located 20 micrometers deeper than the first series, a high correlation between the two series was expected. In this case, the results obtained were not processed using the SSA algorithm in order to compare actual raw Mg/Ca ratio series. As can be seen in Figure 7, the two LAN545 series measured under the same experimental conditions showed the same trend through time, reporting a high correlation (R 2 = 0.78) between them. Overall, this indicates that, while patterns can be interpreted with confidence, some intrinsic variability is inherent in using the LIBS technique, as well as a local heterogeneity in Mg distribution that is not only related to temperature dependence [58,61,69,70].  [65] and García-Escárzaga et al. [69] analyzed Patella vulgata limpets from Atlantic Europe, reporting elemental ratios (from ca. 10 to ca. 30 mmol/mol) lower than obtained herein from P. depressa species by LIBS. We attribute this higher concentration of Mg measured by LIBS techniques to a possible species-dependent variability, but also to the CF-LIBS algorithm itself. The results could be affected by the self-absorption of Ca lines, as self-absorption tends to lead to an underestimation of the concentration of major elements in the chemical composition. The current implementation of the algorithm tries to reduce this effect with a careful selection of the emission lines, but the effect is difficult to remove entirely [79]. Another noticeable aspect of the sequences was their high variability, with apparent high-frequency fluctuations on top of the expected yearly sinusoidal variations. Although much research using the LIBS technique has tended to explain general intrinsic variability in LIBS spectra as a consequence of fluctuations in the varying plasma properties of each laser shot [80], we attribute this behavior to the actual composition heterogeneity of the biogenic carbonates on a micrometric scale, as has been suggested elsewhere [58,61,69]. As other studies report, the mineralization process in limpets may not only be dependent on seawater temperature, but also on other biological and environmental factors, such as mollusk ontogeny, intra-annual variation in shell growth rates, individual physiological processes, seawater salinity, and/or food supply conditions [60,64,81,82]. To evaluate the reliability of measured Mg/Ca ratios, one limpet (LAN545) cross-section was measured twice along the same path. Although the ablated material was not strictly the same because the second series was measured on material located 20 micrometers deeper than the first series, a high correlation between the two series was expected. In this case, the results obtained were not processed using the SSA algorithm in order to compare actual raw Mg/Ca ratio series. As can be seen in Figure 7, the two LAN545 series measured under the same experimental conditions showed the same trend through time, reporting a high correlation (R 2 = 0.78) between them. Overall, this indicates that, while patterns can be interpreted with confidence, some intrinsic variability is inherent in using the LIBS technique, as well as a local heterogeneity in Mg distribution that is not only related to temperature dependence [58,61,69,70].

Mg/Ca Ratio Series vs δ 18 O shell Profiles and Instrumental SST (T meas )
The Mg/Ca ratio series obtained showed very well-defined temporal patterns, covering a variable number of seasonal cycles depending on the age of the shell (Figure 6a-d). A cycle refers to that part of the series between one minimum or one maximum value and the next one. Each cycle is therefore composed of a maximum, a minimum, and a variable number of intermediate values. A variable number of cycles were documented for each specimen. LAN541, LAN554, and LAN559 exhibited elemental fluctuations corresponding to approximately one cycle, while LAN545 covered a time span longer than two cycles. The number of cycles reported by the elemental series for each limpet was very consistent with the actual number of years of the life-span of these specimens, which was previously decoded through a detailed sclerochronological investigation carried out on these four same specimens combining a stable oxygen isotope analysis and a study of the growth patterns on every shell [72].
According to the results obtained from the previous trace element investigations carried out on biogenic carbonates [56,68], and especially those conducted on limpet species [39,57,58,69], higher and lower Mg/Ca ratio values represent warmer and colder SST, respectively, and each cycle corresponds to a complete year. This approach was also confirmed by the results derived from the sclerochronological study conducted on these four shells [72], since the cross-section areas where the minima and maxima Mg/Ca ratio values were located coincided very well with those limpet portions deposited when coldest and warmest seawater temperatures were reached, respectively. A comparison between the obtained Mg/Ca ratio series and instrumentally measured SST (T meas ) at the place where the mollusks grew ( Figure 6e) demonstrated a strong agreement between the two datasets through time, covering the expected number of summer and winter seasons through each mollusk's lifespan. Finally, the Mg/Ca ratio series obtained herein were also compared with δ 18 O shell values previously published for these same specimens [72]. The results reported a very high correlation between both profiles (R 2 = 0.78-0.87) (Figure 8). Therefore, LIBS-derived Mg/Ca ratio series from P. depressa clearly reported SST variations during the lifespan of the mollusk, covering a variable time span between one year (LAN559, Figure 6d) and almost three annual cycles (LAN545, Figure 6b).
In order to quantify the seawater temperature dependence of Mg incorporation into the carbonate matrix, maxima, minima, and average Mg/Ca ratios obtained for each annual cycle, and maxima, minima and average T meas for each year were compared. Linear regression analysis showed very high correlations in each of the four shells (R 2 = 0.91-0.95) (Figure 9a), even when all four shells were considered together in the same analysis (R 2 = 0.86) (Figure 9b). Results obtained herein, therefore, for the first time, show that the P. depressa mollusk species incorporates magnesium into its calcium carbonate matrix in dependence with SST in northern Spain. Furthermore, the LIBS-derived Mg/Ca ratio series have enough resolution to correctly reflect the general evolution of SST variations on a year-round basis, with clearly distinguishable warmer and colder periods. This has implications for both paleoclimate reconstruction as well as for archaeological studies looking to determine seasonal collection strategies of this mollusk.

Sea Surface Temperature Reconstruction and Paleoclimate Implications
Reconstructing environmental conditions and climate change during the past, particularly at a seasonal resolution, is key for better understanding the impact of external changes on different aspects of human behavior [4,5,14]. Previous sclerochronological studies carried out on mollusks collected from northern Spain have revealed that δ 18 O shell values obtained from Phorcus lineatus [51,73], P. vulgata [74], P. depressa [72], and Mytilus galloprovincialis [83] can all be considered as effective high-resolution SST recorders. However, to reconstruct SST from δ 18 O shell values, information about δ 18 O sw is also required, the value of which is effectively unknown for the past. Different methodological approaches have been applied in order to overcome this limitation using, for example, the current δ 18 O sw for SST reconstruction during the Roman period [42], or even the Early Holocene period [84]; deciphering δ 18 O sw from available salinity datasets [41,85]; or correcting the current δ 18 O sw according to known prehistoric sea level changes [86].  Figure 3b) and the sclerochronological investigation previously conducted on these same specimens [72].
In order to quantify the seawater temperature dependence of Mg incorporation into the carbonate matrix, maxima, minima, and average Mg/Ca ratios obtained for each annual cycle, and maxima, minima and average Tmeas for each year were compared. Linear regression analysis showed very high correlations in each of the four shells (R 2 = 0.91-0.95) (Figure 9a), even when all four shells were considered together in the same analysis (R 2 = 0.86) (Figure 9b). Results obtained herein, therefore, for the first time, show that the P. depressa mollusk species incorporates magnesium into its calcium carbonate matrix in dependence with SST in northern Spain. Furthermore, the LIBS-derived Mg/Ca ratio series have enough resolution to correctly reflect the general evolution of SST variations on a year-round basis, with clearly distinguishable warmer and colder periods. This has implications for both paleoclimate reconstruction as well as for archaeological studies looking to determine seasonal collection strategies of this mollusk. To try and overcome this issue, Mg/Ca ratio analysis has been proposed as a reliable proxy independent of δ 18 O sw variations [39,56,68]. However, very few studies have reported a robust dependence between Mg/Ca values and seawater temperatures when carbonate was deposited by the mollusks, and not every species appears to precipitate Mg in equilibrium with environmental surroundings, while the nonquantified effects of other biological and/or environment factors add further dimensions of uncertainty [58][59][60][61]64,65,81]. Therefore, more analyses involving new shell taxa and species collected from coastal areas that were previously not considered are required to better understand trace element precipitation mechanisms. During the last few years, LIBS has emerged as an optimal technique for rapidly investigating if trace element ratios measured on different mollusk species can be used as a climate proxy [57,68,69]. Nevertheless, previous investigations also noticed some limitations of LIBS. These studies showed differences in the Mg/Ca ratio values (from emission line ratios) obtained between specimens, which were caused by inconsistent LIBS experimental conditions and precluded the development of an equation for estimating SST from Mg/Ca values [69]. Here, the LIBS-derived variability was reduced using the CF-LIBS technique, resulting in very similar regression lines between each measured specimen of P. depressa (Figure 9a). This allowed us to tentatively develop the first equation to estimate SST from LIBS-derived Mg/Ca ratios (T Mg/Ca ) considering all maxima, minima, and average values derived from the four shells (Figure 9b): T Mg/Ca = 0.5492 × Mg/Ca ratio value − 8.5532 (1) Nevertheless, further investigation, increasing the number of shells analyzed and improving the CF-LIBS technique algorithm, is still required in order to obtain a more accurate equation.

Sea Surface Temperature Reconstruction and Paleoclimate Implications
Reconstructing environmental conditions and climate change during the past, particularly at a seasonal resolution, is key for better understanding the impact of external changes on different aspects of human behavior [4,5,14]. Previous sclerochronological studies carried out on mollusks collected from northern Spain have revealed that δ 18 Oshell values obtained from Phorcus lineatus [51,73], P. vulgata [74], P. depressa [72], and Mytilus galloprovincialis [83] can all be considered as effective high-resolution SST recorders. However, to reconstruct SST from δ 18 Oshell values, information about δ 18 Osw is also required, the value of which is effectively unknown for the past. Different methodological approaches have been applied in order to overcome this limitation using, for example, the current δ 18 Osw for SST reconstruction during the Roman period [42], or even the Early Holocene period [84]; deciphering δ 18 Osw from available salinity datasets [41,85]; or correcting the current δ 18 Osw according to known prehistoric sea level changes [86].
To try and overcome this issue, Mg/Ca ratio analysis has been proposed as a reliable proxy independent of δ 18 Osw variations [39,56,68]. However, very few studies have re-

P. depressa Season of Collection and Archaeological Implications
Beyond paleoclimate conditions, determining the season in which archaeological mollusks were collected by prehistoric foragers also has significant implications for better understanding past subsistence strategies, as well as providing valuable information relating to the seasonality of site occupation and mobility patterns of past societies [46,[87][88][89]. The LIBS technique can offer key valuable advantages for future studies. Not only does LIBS avoid sample pretreatments, but it also avoids highly time-consuming sampling procedures in order to extract carbonate powders along the axes of shell growth [50,74,90,91]. García-Escárzaga et al. [69] showed that it offers a 20× reduction in the time required to measure a specimen, with respect to traditional ICP-based methodologies. In addition, the LIBS technique also offers a notable cost reduction relative to isotope ratio mass spectrometry techniques. The time and cost reductions enable a notable increase in the number of specimens analyzed for each stratigraphic unit, offering more statistically consistent results through time and between sites. Moreover, as can be seen from the data present here, and in other studies, the ability of Mg/Ca ratio series to decipher past climate changes on an intra-annual basis means that the approach is able to reveal seasonality in mollusk collection patterns [57,58,69].
To decipher the season of collection of each specimen, the most recent (youngest) calcium carbonate deposited by the mollusk, as well as the sequence trend prior to the mollusk's death, were employed [44,47,87]. Mg/Ca ratio values representative of the last mollusk growth in the four shells always reported a value higher than 50 mmol/mol. However, slight differences were observed in the series trend during the weeks or months prior to the collection event. While LAN541 showed a harvest event during a period of cooling SST, just after the annual SST maximum was reached, LAN545, LAN554, and LAN559 reported a collection event during a warming trend in SST, just before reaching the maximum value documented during the previous annual cycle. Only LAN559 reported a final value (54 mmol/mol) similar to the maximum observed during the previous summer (56.5 mmol/mol), which suggests a collection near the warmest annual SST. The sclerochronological study previously carried out on these same four shells revealed that LAN545 and LAN554 specimens stopped their growth for one and two months, respectively, in summer 2012. Both specimens resumed carbonate deposition at the end of September 2012, just before the collection event [72]. This explains why LAN545 and LAN554 did not record the annual maximum SST. In spite of these growth stoppages in LAN545 and LAN554, the last Mg/Ca ratio value of both specimens was always located in the upper quartile of the total annual range reported from each individual [44,87,92], which corresponded to the summer season.
Overall, if the real collection date of these specimens (i.e., end of summer) were unknown, the series obtained would have allowed us to decipher, with a high temporal resolution, the period of the collection in the case of LAN541 and LAN559, which clearly reported a collection during the end of the warmer season. Although LAN545 and LAN554 did not reach the annual maximum SST, the Mg/Ca ratio series also correctly suggested the general season of collection (i.e., summer). As also occurred with the isotopic series obtained from these same two specimens (LAN545 and LAN554) [72], a collection at the end of summer could have also been correctly deduced if Mg/Ca ratio series were combined with growth line analyses in order to take into account possible growth stoppages close to the edge. In summary, our results highlight that the season of collection for P. depressa can be accurately, rapidly, and inexpensively determined using LIBS-derived Mg/Ca ratio series and multiproxy sclerochronological investigations, combining trace element analyses and incremental shell growth patterns. The LIBS technique, therefore, offers a significant opportunity to effectively increase the number of shells measured for each archaeological context, obtaining thus more accurate archaeological conclusions.

Conclusions
For the first time, we applied the CF-LIBS technique to biogenic calcium carbonate to estimate molar concentrations of chemical elements. Our results reported slight differences with molar concentration data derived from ICP-based methodologies, suggesting that algorithms developed to conduct CF-LIBS must be improved in further studies, particularly, to deal with the self-absorption effect that tends to overestimate the concentration of minority elements. Nevertheless, Mg/Ca ratio sequences derived from P. depressa using this novel LIBS technique correctly reflected annual variations in SST, demonstrating that P. depressa incorporates magnesium into its calcium carbonate matrix in dependence with SST in northern Spain. The high correlation between Mg/Ca series and SST at the place where the mollusks grew allowed us to obtain a tentative equation to estimate SST from measured Mg/Ca ratios, although further studies are still required utilizing more samples and improving the CF-LIBS algorithm. Together then, this novel methodology represents a valuable approach for paleoclimate applications in the region, as well as a means of estimating of the season when mollusks from archaeological sites were collected by past humans with implications for the study of human-environment interactions.