Calcification response of planktic foraminifera to environmental change in the western Mediterranean Sea during the industrial era

The Mediterranean Sea sustains a rich and fragile ecosystem currently threatened by multiple anthropogenic impacts that include, among others, warming, pollution, and changes in seawater carbonate speciation associated to increasing uptake of atmospheric CO₂. This environmental change represents a major risk for marine calcifiers such as planktonic foraminifera, key components of pelagic Mediterranean ecosystems and major exporters of calcium carbonate to the sea floor, thereby playing a major role in the marine carbon cycle. In this study, we investigate the response of planktic foraminifera calcification in the northwestern Mediterranean Sea on different timescales across the industrial era. This study is based on data from a 12-year-long sediment trap record retrieved in the in the Gulf of Lions and seabed sediment samples from the Gulf of Lions and the promontory of Menorca. Three different planktic foraminifera species were selected based on their different ecology and abundance: Globigerina bulloides, Neogloboquadrina incompta, and Globorotalia truncatulinoides. A total of 273 samples were weighted in both sediment trap and seabed samples.

The results of our study suggest substantial different seasonal calcification patterns across species: G. bulloides shows a slight calcification increase during the high productivity period, while both N. incompta and G. truncatulinoides display a higher calcification during the low productivity period. The comparison of these patterns with environmental parameters indicate that controls on seasonal calcification are species-specific. Interannual analysis suggests that both G. bulloides and N. incompta did not significantly reduce their calcification between 1994 and 2005, while G. truncatulinoides exhibited a constant and pronounced increase in its calcification that translated in an increase of 20 % of its shell weight. The comparison of these patterns with environmental data reveals that optimum growth conditions affect positively and negatively G. bulloides and G. truncatulinoides calcification, respectively. Sea surface temperatures (SSTs) have a positive influence on N. incompta and G. truncatulinoides calcification, while carbonate system parameters appear to affect positively the calcification of three species in the Gulf of Lions throughout the 12-year time series.

Finally, comparison between sediment trap data and seabed sediments allowed us to assess the changes of planktic foraminifera calcification during the late Holocene, including the pre-industrial era. Several lines of evidence indicate that selective dissolution did not bias the results in any of our data sets. Our results showed a weight reduction between pre-industrial and post-industrial Holocene and recent data, with G. truncatulinoides experiencing the largest weight loss (32 %–40 %) followed by G. bulloides (18 %–24 %) and N. incompta (9 %–18 %). Overall, our results provide evidence of a decrease in planktic foraminifera calcification in the western Mediterranean, most likely associated with ongoing ocean acidification and regional SST trends, a feature consistent with previous observations in other settings of the world's oceans.

1 Introduction

Growing population and its linked human activity since the industrial period (defined according to Sabine et al., 2004 from 1800 and therein) has caused an increase in carbon dioxide, the ecological and economic consequences of which are considered a major threat (IPCC, 2022). Atmospheric CO₂ concentrations during the Pleistocene and Holocene ranged from 200 to 280 parts per million (ppm) (Loulergue et al., 2007; Lüthi et al., 2008; Parrenin et al., 2007), but these values have increased dramatically since the onset of the industrial period, exceeding the threshold of 400 ppm in 2015 for the first time for at least the last 800 000 years (Lüthi et al., 2008). This increase is significantly more important since the 1950s, when rapid atmospheric changes due to human activity took place, a process referred to as the “Great Acceleration” (Head et al., 2022a). Since then, between, 25 % and 30 % of anthropogenic CO₂ has been absorbed by the world's ocean (Sabine et al., 2004). The ocean uptake of atmospheric CO₂ causes a drop in both pH and carbonate ion concentration (Barker and Elderfield, 2012), lowering seawater alkalinity; this process is commonly known as ocean acidification (OA), and it is expected to affect all areas of the ocean and to have a wide impact on marine life (Davis et al., 2017; Figuerola et al., 2021; Orr et al., 2005). One of the main questions about recent environmental change is how different ecosystems and regions in global oceans are going to react to the ongoing increase of anthropogenic atmospheric carbon dioxide.

A large body of evidence indicates that ocean acidification has substantial and diverse effects on the distribution and fitness of a wide range of marine organisms (Kroeker et al., 2013; Meier et al., 2014; Moy et al., 2009). For example, some fleshy algae and diatom species have been shown to increase their growth and photosynthetic activity at enhanced CO₂ concentrations (Kroeker et al., 2013). In turn, most calcifying organisms such as calcifying algae, corals, pteropods, coccolithophores, and foraminifera are negatively affected by this process, often showing a reduction in their abundance, calcification, and growth rates (Kroeker et al., 2013; Orr et al., 2005).

Planktic foraminifera are a group of marine single-celled protozoans that produce calcareous shells. Their distribution across the water column is conditioned by factors that include, but are not limited to, food availability, temperature, salinity, and sunlight (Schiebel and Hemleben, 2005). These organisms are considered to play a key role in marine carbon cycle and carbonate production, accounting for between 32 % and 80 % of the deep ocean calcite fluxes (Schiebel, 2002). Depending on their ecology and feeding strategies, these organisms can be algal (dinoflagellates) symbiont bearing or not symbiont bearing and be spinose or non-spinose. Planktic foraminifera represent a useful tool for palaeoecological and palaeoceanographic studies, as the abundances of different species and their geochemical signature allow reconstructing sea surface temperatures (SSTs) and water column physical and chemical properties (Lirer et al., 2014; Margaritelli, 2020; Schiebel and Hemleben, 2017).

Previous studies suggest that planktic foraminifera are sensitive to ocean acidification (OA). Laboratory experiments indicate that when carbonate ion concentration decreases, shell weight and calcification decrease too in a variety of species (Bijma et al., 2002; Lombard et al., 2011). Species that host symbionts have been described showing a higher tolerance to dissolution due to the capacity of algal symbionts to alter immediate seawater chemistry (Lombard et al., 2009). Moy et al. (2009) documented a decrease of 30 %–35 % shell weight in the planktic foraminifera Globigerina bulloides during the industrial era in the subantarctic Southern Ocean, most likely induced by anthropogenic-driven ocean acidification. A recent study by Fox et al. (2020) showed that non-spinose (Neogloboquadrina dutertrei) foraminifera species exhibit a more pronounced calcification reduction than the spinose (Globigerinoides ruber) species in response to increasing CO₂. The main difficulty for studying the impact of OA on foraminifera (and any calcifying organisms) resides in finding long-term continuous records in order to be able to evaluate possible changes in shell calcification (Fox et al., 2020).

In order to assess the impact of recent environmental change on planktic foraminifera, in this work we present data from Planier sediment trap (data from 1993 to 2006) (Rigual-Hernández et al., 2012) and from seabed sediments from three different sites located in both the Gulf of Lions and the promontory of Menorca. The Mediterranean Sea is a semi-enclosed sea with a high saturation state for calcite (Álvarez et al., 2014). It is often considered as a “miniature ocean” and a “laboratory basin” (Bergamasco and Malanotte-Rizzoli, 2010), which makes it a valuable zone to study the response of marine calcifying organisms to environmental change.

The advantage of sediment traps is that they can provide data coming from annual fluxes, avoiding the effects of seasonal abundance and ontogeny and making interannual comparisons more reliable (Jonkers et al., 2019). Three different planktic foraminifera species, each of which characterized by contrastingly different depth habitats and ecologies, were selected for our analysis: Globigerina bulloides, a spinose opportunist surface dweller that lies above the thermocline; Neogloboquadrina incompta, a non-spinose temperate surface dweller; and Globorotalia truncatulinoides, a non-spinose deep-dwelling species which migrates through the water column with a complex life cycle (Schiebel and Hemleben, 2017). Our aims for this study are: (i) to compare two widely used foraminifera weighing and size-normalization techniques and provide a baseline of modern foraminifera weight data and calcification in the western Mediterranean against which future changes in foraminifera calcification can be assessed; (ii) to document seasonal and interannual trends in the planktic foraminifera calcification of the three planktic foraminifera species; and (iii) to evaluate possible changes in shell calcification through the Holocene to the present day by comparing shell weights of the foraminifera collected by the traps with those of the seafloor sediments.

2 Study area

The Mediterranean is a semi-enclosed sea and is considered a concentration basin (Bethoux et al., 1999) with a negative hydrological budget: fresh water inputs do not compensate the overall basin evaporation. The surface oceanic waters that enter the Mediterranean through the Strait of Gibraltar and spread towards the eastern basin compensate this negative balance. The waters of Atlantic origin (AW) circulate according to a cyclonic circuit (anticlockwise) along the Mediterranean rim (Fig. 1a). In the northwestern basin, this along-slope current, called the northern current (NC), is formed by the convergence at the level of the Ligurian Sea (LS) of the waters flowing on both sides of Corsica. The general circulation in this sub-basin forms a cyclonic pattern, flowing westward along the Gulf of Lions, bypassing the Balearic Sea and finally closing its circuit eastward along the North Balearic Front (NBF) (Fig. 1b).

Figure 1 (a) Study area location in the Mediterranean Sea and general surface circulation (b) geographic setting of the Gulf of Lions and location of studied sites. Red diamond shows the position of the Planier site sediment trap and core top (PLA ST/CT). Black diamonds represent the location of the seabed sediments samples analysed from Lacaze–Duthiers canyon (LCD SC) and Menorca promontory (MR 3.1A). Yellow diamond represents the location of the Dynamics of Atmospheric Fluxes in the MEDiterranean Sea (DYFAMED) site, located 200 km upstream of Planier station position. Arrows represent the surface circulation (Millot, 1999). The topographic model was downloaded from the GEBCO database.

Moreover, the Mediterranean is recognized as a sensitive region to increasing atmospheric CO₂ due to the fast turnover time of its waters (Béthoux et al., 1999) and the fast penetration of anthropogenic CO₂ (Schneider et al., 2007). Sea surface temperatures are predicted to increase by 1.5–2 ^∘C by the end of the century, a faster rate than the global average (Lazzari et al., 2014). pH is expected to decrease according to the global average (0.3–0.4 units by 2100) or even exceed the global trend (Hassoun et al., 2015). The Mediterranean Sea is also affected by other stressors, which impact marine organisms in many ways (Lejeusne et al., 2010). Finally, it is also a region shaped by human development, and its associated activities interact with environmental changes (Mediterranean Experts on Climate Change; MedECC, 2020).

The Gulf of Lions is located in the northwestern part of the Mediterranean Sea, and its morphology presents a continental slope with an array of complex submarine canyons (Rigual-Hernández et al., 2012) (Fig. 1b). Vertical mixing, generated by intense surface cooling and evaporation, occurs in winter in the Gulf of Lions driven by cold, dry northern winds, resulting in dense water on the shelf and offshore (Durrieu de Madron et al., 2013; Houpert et al., 2016; Millot, 1990). This winter mixing recharges surface waters with nutrients. This enrichment with increased solar radiation stimulates primary production in spring. Increasing heat fluxes during spring and summer cause water mass stratification and nutrient depletion, which lasts until late summer, until fall cooling breaks the stratification of the water column and causes a fall bloom (Heussner et al., 2006; Monaco et al., 1999; Rigual-Hernández et al., 2012). River inputs are the main source of suspended particles in the Gulf of Lions, and the Rhone River represents the most important river in the northwestern Mediterranean; however, other sources include Saharan dust deposits and biological production (Heussner et al., 2006; Monaco et al., 1999). Overall, the oceanographic setting of the Gulf of Lions is an exception to the general oligotrophy of the Mediterranean Sea.

3 Material and methods

3.1 Sediment traps, core tops, and sediment cores

A series of deployments of sediment traps mooring lines in the Gulf of Lions continental margin was initiated in 1993 within the framework of several French and European projects (PNEC, Euromarge-NB, MTP II-MATER, EUROSTRATAFORM) and the monitoring of two sites, Planier and Lacaze–Duthiers stations (Fig. 1), continues as a component of the MOOSE (Mediterranean Ocean Observing System for the Environment) programme (Coppola et al., 2019). Planier station (43^∘02^′ N, 5^∘18^′ E) is located at the northeastern end of the Gulf of Lions, in the axis of the Planier Canyon. The sediment trap used for this work was located at around 530 m water depth in a water column of ∼ 1000 m. Further details of the mooring design can be found in. Heussner et al. (2006). Planktic foraminifera fluxes for the 1993 to 2006 period were documented by Rigual-Hernández et al. (2012). Here, we used the samples from the latter study for our weight and calcification analysis. This sediment trap is used here as a baseline of the planktic foraminifera dwelling in the modern Mediterranean Sea. Moreover, we analysed a set of core top and sediment cores collected from several locations of the northwestern Mediterranean that are considered to represent foraminifera assemblages sedimented during the Holocene era (Table 1).

3.2 Sediment core samples processing

A total of two sediment samples from Planier core top, seven from Lacaze–Duthiers sediment core and 40 from Minorca sediment core, were weighed (Table 1). Dry bulk sediment samples from all sites were weighed using a Sartorius CP124S balance (precision of 0.1 mg).

The samples were then wet-sieved in order to separate the < 63 µm fraction and dry sieved to separate the bigger fractions (> 150 and > 300 µm). The sediment washing was carried out with potassium phosphate-buffered solution (pH of 7.5) in order to optimize foraminifera preservation. Each fraction was oven-dried at a constant temperature (50 ^∘C) and then weighed. The > 150 µm fraction was used for identification, counting, and shell morphometric and weight analyses.

Table 1 Description of the core tops used in this study. Data for Planier core top (PLA CT) and Lacaze–Duthiers sediment core (LCD SC) are available in Heussner et al. (2006), and data concerning Menorca sediment core (MR 3.1.A) can be found in Cisneros et al. (2016). Conventional ¹⁴C ages, 1σ uncertainties, local reservoir, and the calibrated age have been rounded according to convention. See Sect. 3.7 for details concerning “Bomb ¹⁴C”.

3.3 Ecology and life cycle of G. bulloides, N. incompta, and G. truncatulinoides

G. bulloides is a spinose surface to subsurface dweller (Schiebel and Hemleben, 2017), found in the upper 60 m of the water column. This species has affinity for temperate to subpolar waters and upwelling systems in lower to midlatitudes (Azibeiro et al., 2023; Bé et al., 1977). In terms of its seasonal distribution, it has been documented to be associated to enhanced productivity periods in mid to high latitudes (Chapman, 2010; Schiebel and Hemleben, 2005). No symbiont algae are hosted by this species, and, contrary to most spinose species, its diet is mainly algae based (Schiebel et al., 2001). G. bulloides shows an opportunistic feeding and strategy, leading to a high abundance in the foraminifera assemblages preserved in the sedimentary record. This is despite tests that have documented it to be more susceptible to dissolution than the average planktic foraminifera species (Dittert et al., 1999).

N. incompta is a surface dweller abundant in subpolar to temperate water masses across all the ocean basins (Kuroyanagi and Kawahata, 2004). It is a non-spinose species and does not carry symbiont algae. In North Atlantic waters, N. incompta is a major component of foraminifera assemblages from late spring to late fall, and generally it is the dominant foraminifera species during late summer when maximum shoaling of mixed layer depths occur (Schiebel and Hemleben, 2000). It shows a minor presence in low latitudes and during periods of enhanced nutrient supply, N. incompta is outnumbered by other more opportunistic species (Schiebel et al., 2002).

G. truncatulinoides is considered the deepest dweller among the extant planktic foraminifera, with living specimens documented below 2000 m (Schiebel and Hemleben, 2005). Considered a widespread species, it can be found from subpolar to subtropical water masses (Schiebel and Hemleben, 2017). It is a non-spinose species (Margaritelli, 2020), and it does not carry any symbiont algae (Takagi et al., 2019). An important aspect to highlight about this species is its complex life cycle (Margaritelli et al., 2022). It reproduces once a year in the upper water column during late winter, when mixing of the water column allowed the migration of juveniles to the surface waters (Lohmann and Schweitzer, 1990; Schiebel et al., 2002). The former authors speculated that nutrient availability and the avoidance strategies to predation could explain its life cycle. Then, the adults migrate downward the water column (Rebotim et al., 2017) and spend the rest of their life cycle developing an additional calcite layer in cooler waters below the thermocline (Lohmann and Schweitzer, 1990; Wilke et al., 2009). Around 70 % of G. truncatulinoides calcification has been estimated to take place around the thermocline, while the remaining 30 % takes place in surface waters (LeGrande et al., 2004).

The abundance of these three species has been previously studied in the Gulf of Lions by Rigual-Hernández et al. (2012). The latter study showed that both G. bulloides and N. incompta displayed their maximum abundances during the spring bloom, while the G. truncatulinoides abundance was maximum during early winter. On the other hand, minimum abundances were reached during late spring and summer for G. bulloides and N. incompta, respectively, and G. truncatulinoides displayed a minimum abundance during fall.

3.4 Foraminifera picking and mass and size estimations

Different sizes were selected depending on the maximum availability of each species: 250–300, 200–250, and 400–500 µm for G. bulloides, N. incompta, and G. truncatulinoides, respectively. For the latter species, both coiling morphotypes were selected, although the right coiling was substantially less abundant, representing less than 3 % in our counts, a feature consistent with the literature that indicates a low presence of right-coiled specimens (Margaritelli et al., 2020, 2022).

A total of 273 foraminifera samples were picked for this study, 126 coming from the sediment trap and 147 from the three sediment cores and core tops (Table 2). However, these numbers represent the total of samples analysed, but the unique sample number is lower, as not all the sediment trap samples presented the three species in high enough numbers to perform the picking. The species were analysed in size fractions in order to estimate the efficiency of sieve fractions and the impact of size and morphometric parameters on the foraminifera weight and calcification.

The mean weight of each available sediment trap sample was obtained by weighting between 15 to 45 specimens of G. bulloides (mean: N=27), 5 to 25 of N. incompta (mean: N=15), and 5 to 25 of G. truncatulinoides (mean: N=13). Concerning the analyses of the core top and sediment core samples, between 15 and 25 G. bulloides and N. incompta (mean: N=20 for both) and between 9 and 25 G. truncatulinoides (mean: N=18) were picked.

Each foraminifera sample was then exposed to gentle ultrasonication (50 Hz) for 5 to 75 s (depending on the species and the degree of visual uncleanliness) in methanol in order to clean the shells. The samples were then left to dry in a temperature-controlled oven at 50 ^∘C. One out of three analysed samples were weighted before and after cleaning in order to assess potential shell mass losses and shell preservation due to ultrasonication. Our results indicate that this method has little impact on shell preservation, with around 95 % of the total foraminifera conserved in good conditions. Weight loss between non-cleaned samples and cleaned samples is a mean of 0.5 to 3 µg (between 6 % and 32 % of the sample's total measured weight) depending on the species, mainly due to the presence of clay and non-calcite material in the shells, which justifies this cleaning process (see Fig. S6 in the Supplement).

The weightings were carried in the micropaleontology laboratory of the Geology Department at the University of Salamanca using a Sartorius ME5 balance (precision of 0.001 mg). This method allowed us to obtain foraminifera sieve-based weight (SBW) by dividing the average shell weight per sample (5–45 tests) by the total number of foraminifera within each sample. The lowest number of individuals selected per sample was five, in order to maximize the number of samples available for our study. According to Beer et al. (2010a), the higher the number of individuals, the more reliable the SBWs are. Here we aim to compare SBW results with a measured weight technique. Measured techniques are acknowledged to be reliable with a lower number of individuals, therefore a minimum of five individuals were selected in order to compare the two techniques.

However, it has been described that the traditionally-used sieve fractions method is considered unreliable because of the effect of morphometric parameters on the foraminifera weight (Beer et al., 2010b). In order to remove the size effect on the weight, the mean SBW was normalized to the mean diameter and area of the planktic foraminifera to obtain measurement-based weights (MBW). Morphometric parameters were measured using a Nikon SMZ18 stereomicroscope equipped with a Nikon DS-Fi3 camera and NISElements software. These measurements were carried out on the same shells that were weighted. Foraminifera shells were positioned in order to obtain the maximum area of each individual; in this case, the umbilical side (aperture facing upwards) was measured for the three species.

MBW_area and MBW_diameter were calculated according to the following formula (Aldridge et al., 2012; Beer et al., 2010a), where parameter accounts for area or diameter:

$${\mathrm{MBW}}_{\mathrm{sample}}={\displaystyle \frac{{\text{mean SBW}}_{\text{sample}}\times {\text{mean parameter}}_{\text{size fraction}}}{{\text{mean parameter}}_{\text{sample}}}}.$$

Size fraction accounts for the mean of the parameter (area or diameter) measured in all the sites studied, while sample accounts for the mean of the parameter in the particular sample being measured. The advantage of these measurements is that the resulting MBW is being given with a weight unity (µg), thereby allowing direct comparison with other studies (Beer et al., 2010a, b), and it is useful for estimating their contributions to marine biogeochemical cycles.

Correlations between SBW and MBW_area against area are displayed in Fig. 2. The reason for this comparison is to show the relation between size and weight and to avoid the effect of having the bigger specimens displaying the heaviest weight and impacting the mean weight (and therefore the calcification indicator) of the sample.

Finally, in order to compare weight patterns from the sediment trap with weights from core tops and sediment cores and overcome the seasonality effect, MBW were flux-weighted. Mean monthly MBWs values from each species were multiplied by the corresponding mean monthly flux and then divided by the total annual flux of the corresponding species. For these calculations, the flux data from each species estimated for the > 150 µm fraction from Rigual-Hernández et al. (2012) were employed.

3.5 Environmental data

Foraminifera fluxes and abundances together with chlorophyll a were taken from Rigual-Hernández et al. (2012) for the entire time span of the analysed samples. Both fluxes and abundance come from direct sediment observation from the Planier site, while chlorophyll a data were obtained from SeaWiFS monthly measurements through NASA's Giovanni programme on a $\mathrm{0.2}\times \mathrm{0.2}{}^{\circ }$ area around the mooring location. SeaWiFS measurements started in 1997 and were used due to the lack of direct chlorophyll measurements in our samples. Sea surface temperature (SST) was recovered from the NOAA database with the same gridding as the data from the NASA's Giovanni programme.

Salinity, nutrient concentrations (nitrates and phosphates), and carbonate system parameters data were collected from the DYFAMED database (http://www.obs-vlfr.fr/dyfBase/index.php, last access: 26 December 2022) (Coppola et al., 2008, 2021). The DYFAMED site is located around 200–220 km (Fig. 1b) east of the sediment trap location (43^∘25^′ N, 7^∘52^′ E), in the Ligurian Sea. From an oceanographic view, its situation is upstream of the NC circulation and can be considered representative of seasonal and interannual variability of biological and water column properties of the open-ocean waters in the northwestern Mediterranean (Heussner et al., 2006; Meier et al., 2014). Alkalinity and total carbon measurements were available for years 1998 to 2000 and mid-2003 to 2005. Missing values comprised in these years were replaced with values obtained from linear regression of the measurements from above and below. The CO2SYS macro has been used to reconstruct the [CO₂], [CO${}_{\mathrm{3}}^{\mathrm{2}-}$], [HCO${}_{\mathrm{3}}^{-}$] and pH values from the measured total alkalinity and dissolved inorganic carbon. The constants used were the CO₂ dissociation constant by Mehrbach et al. (1973) refit by Dickson and Millero (1987), the KHSO₄ by Dickson (1990), and the seawater scale for pH.

3.6 Statistical analysis

In order to have uninterrupted monthly environmental values from the DYFAMED site during available measurements, a resampling every 10 d has been carried out with the QAnalySeries programme.

Independence and correlation between the area of the different species SBWs and MBW_area was tested using a Pearson linear correlation test with an R script (see Supplement).

Seasonal correlation analyses were carried out with the Statistica programme. A p< 0.05 was used in order to consider a correlation as significant. The number (N) of correlations depended on data availability and was 10 for G. bulloides, 9 for N. incompta, and 12 for G. truncatulinoides.

It should be noted that the analysis of interannual trends was hindered by gaps in the sediment trap record and by the low number of specimens during some sampling intervals. Therefore, interannual trends in planktonic foraminifera calcification should be interpreted with caution.

The influence of environmental variables upon MBW_area was assessed using general additive models (GAMs) (fitted using the gam function from the mgcv R package). Due to data limitation, the GAMs could not be fitted to multiple independent variables, so potential effects of interacting environmental variables could not be assessed. Each model tested the dependence of the different MBW upon a single independent variable: month or year, to evaluate seasonal and interannual trends; the flux of each species, to test effects of ecological variability; and a suite of environmental variables, to determine impacts of various aspects of ocean chemistry on the calcification. Smooth functions of these measured quantities were used as the single independent variable within the GAMs, which were fitted using the default settings of the gam function, a Gaussian family and identity link function, and the GCV.cp smoothing method. GAM results quantified the significance of the effect of each independent variable upon MBW.

In order to investigate the difference between the MBW data from the sediment trap and the core top/sediment cores, a non-parametric two-way Mann–Whitney test has been performed. This test determines if there are significant differences in the medians of data sets without making assumptions about the data distributions. A p value < 0.05 has been used to consider the median of two data sets differently.

3.7 Radiocarbon dating

Between 50–100 individuals of well-preserved G. bulloides shells (> 150 µm) were picked for radiocarbon analyses. Radiocarbon (¹⁴C $/$ ¹²C) was measured as CO₂ with a gas ion source in a mini carbon dating system (MICADAS) at the Laboratory of Ion Beam Physics from ETH Zürich. The employed automated method consists of initial leaching of the outer shell to remove surface material with 100 µL of ultrapure HCl (0.02 M) and the subsequent acid digestion of the remaining carbonates with 100 µL of ultrapure H₃PO₄ (85 %) (Wacker et al., 2013). Therefore, no cleaning was applied after the picking, contrary to the samples used for mass and size measurements. Marble (IAEA-C1) was used for blank correction, and results were corrected for isotopic fractionation via ¹³C $/$ ¹²C isotopic ratios.

Conventional radiocarbon age for a sample of 14–14.5 cm from the MR 3.1.A site was calibrated with the online calibration program CALIB (Stuiver and Reimer, 1993) using the Marine20 curve, which applied a marine reservoir correction of 550 ¹⁴C years (Heaton et al., 2020) to the corresponding ¹⁴C age and error. Additionally, a local reservoir effect (Stuiver and Braziunas, 1993) of −165 × 95 ¹⁴C years was considered. This local reservoir was calculated as the average of the eight nearest points to the sample location from the marine reservoir correction database (Reimer and Reimer, 2001), whose values have already been corrected for the Marine20 curve. ¹⁴C ages from samples 0.5–1 cm from both PLA CT and LCD SC lied out of the range for calendar calibration, implying these samples contain some bomb ¹⁴C and cannot be considered pre-industrial (Table 1, see “bomb ¹⁴C”). In order to have an estimation of the time span that could be covered by these dates, the same marine and local reservoir corrections were applied to the most recent ¹⁴C date that could be corrected (i.e. 603 ¹⁴C yr BP). As the F¹⁴C for this sample was < 1 (see Table S1 in the Supplement), this means that the ¹⁴C found in these samples is not dominated by the bomb carbon. Here we propose a 110–50 cal. yr BP age for these samples. Then, these samples could be considered post-industrial. The detailed results of the calibration and the ¹⁴C dating can be found in the Supplement in Figs. S1 and S2.

Finally, it is important to consider that these ¹⁴C ages represent mean average values. Therefore, time integration within each sample and the effects of bioturbation could cause a variation on the foraminifera real ages (Dolman et al., 2021).

Both the samples and dates obtained are detailed in Table 1. Planktic foraminifera present in the dated samples that were not selected for radiocarbon dating were also analysed following the methodology described previously.

Table 2 Minimum, mean, maximum, and standard deviation values of shell area, diameter, and SBW for G. bulloides, N. incompta, and G. truncatulinoides at all studied sites. The last three columns show the number of samples, the mean number (N) of individuals analysed per sample, and the total number of individuals measured for each site.

4 Results

4.1 Shell morphometric parameters and shell-weight normalization

Overall, the mean values for both diameter and area correspond to the mean narrowed size fraction used during the picking, but morphometric parameters show some variability between the studied sites. Standard deviation of both area and diameter values for the three species are higher in the sediment trap record than in seafloor sediments, with mean values (of all three species) of 82 % higher for area and 69 % higher for diameter. SBW exhibits the same pattern, as both area and diameter standard deviation is a mean 130 % higher in the Planier sediment trap. Regarding the variability across the seafloor samples, Planier core top exhibits a greater area and diameter values (about 40 % to 50 % increase for the three species) compared to those of the other two sediment cores, probably due to the fewer samples analysed (Table 2).

The Planier sediment trap results (Table 2) show a higher standard deviation for both area and diameter for the three species, i.e. 76  % and 68 % higher for G. bulloides compared to the data from core tops, 78 % and 54 % for N. incompta, and 81 % and 73 % for G. truncatulinoides.

Because of the lack of precision of the initial individual picking, carried out with a micrometre installed in the microscope, the selection is not totally accurate. Due to this issue, one-third of the of the total measured foraminifera (i.e. 1645 of 4694) were out of the desired size fraction, of which 12 % were bigger ($\mathrm{580}/\mathrm{4694}$) and 23 % were smaller ($\mathrm{1065}/\mathrm{4694}$). Nonetheless, only 0.02 % were more than 20 % out of the selected size range ($\mathrm{64}/\mathrm{4694}$ more than 20 % bigger and $\mathrm{29}/\mathrm{4694}$ more than 20 % smaller). Mean size difference for the foraminifera out of the size fraction is around 7 %. Results vary according to the site and the species. Individuals of 50 % from the Planier sediment trap ($\mathrm{1046}/\mathrm{2077}$) and 26 % of the individuals coming from the core tops ($\mathrm{692}/\mathrm{2617}$) were out of range. G. bulloides showed a 45.5 % (53.2 % in the sediment trap and 39.3 % in the core tops samples) of individuals out of selected size fraction, while this value was 21.5 % (22.2 % in sediment trap, 21.1 % in sediment cores) for N. incompta and 35 % for G. truncatulinoides (53.4 % in sediment trap, 16.7 % in sediment cores).

Even though a narrow size class was selected for each species (see Sect. 3.4), a clear influence of the area on the SBW was found in our data set (Fig. 2).

Figure 2 SBW in µg and MBW_area in µg against the mean test area in µm² for foraminifera samples in the Planier sediment trap. Dark blue dots correspond to G. bulloides, black dots to N. incompta, and red dots to G. truncatulinoides.

In particular, SBW shows a positive correlation with area: 0.33 < r² < 0.53 (Fig. 2). This indicates that the SBW is dependent on the size of the specimens within the selected size range. Thus, to isolate the component of variation in foraminifera shell thickness that represents a change in calcification and does not occur as a direct result of changes in shell size, normalization of the shell weight was performed following the formula detailed in Sect. 3.3 (Beer et al., 2010a). After normalization, MBW_area shows no correlations with area: $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ < r² < 0.12 (Fig. 2). Note that the weight variations in our data set are quite considerable, especially for G. truncatulinoides, probably due to the wider size fraction. Diameter does show correlation with SBW: 0.33 < r² < 0.5, and it shows a non-negligible correlation with MBW_diam: 0.2 < r² < 033. Our data demonstrate that SBW correlates more strongly with MBW_diam than with MBW_area for the three species: 0.9>0.48 for G. bulloides, 0.89>0.52 for N. incompta, and 0.97>0.81 for G. truncatulinoides. These values are consistent with previous studies (Beer et al., 2010).

Table 3 Pearson correlation test results for the three species' correlation between area (µm²) and both SBW (µg) and MBW_area (µg). Here ci stands for “confidence interval”. Significant r values (0 < ci < 1) are set in bold.

Furthermore, a Pearson correlation test (see Sect. 3.6) has been carried out in order to assess the influence of area on SBW and MBW_area (Table 3). Results showed that the SBWs from the three species correlated positively and significantly (0 < ci < 1) with their corresponding areas (0.57 < r < 0.72). Concerning the MBWs, no significant (0 > ci > 1) correlations with the area are observed (−0.15 < r < 0.2). Therefore, these correlations further highlight the fact that SBW values are significantly influenced by shell area, while MBW_area values appeared to be independent of the area.

Differences between SBW and both MBW_area vary depending on the species: SBW is slightly heavier for G. bulloides, heavier for N. incompta, and lighter for G. truncatulinoides. The mean standard deviation for all three species is around 8 %: 7.8 % for G. bulloides, 6.4 % for N. incompta, and 13 % for G. truncatulinoides. We take these values as the error adjustment for SBW in the different size fractions (250–300, 200–250, and 400–500 µm, respectively). It is difficult to compare these results with other studies, as size fractions and species are often different, but these error estimates are in the same order of magnitude as some other MBW published in core top records and sediment traps (de Moel et al., 2009; Moy et al., 2009).

These findings highlight the fact that the use of sieve fractions does not provide enough control on the influence of morphometric parameters in test weight. Morphometric variations described in Table 1 indicate that the typically used sieve fractions may be unreliable due to the number of individuals out of the desired fractions and the variability within the size range. The correlations between SBW and shell area are consistent with previous studies (Aldridge et al., 2012; Beer et al., 2010a) and underscore the importance of isolating the component of variation in foraminifera shell thickness that represents a change in calcification and does not occur as a direct result of change in shell size. Thus, the shell weight was size-normalized after Beer et al. (2010a) by isolating the influence of isometric scaling on wall thickness and calcification density.

Moreover, both MBW_area (Fig. 2) and MBW_diam, in either the sediment trap data and core top data, do not correlate with area and diameter ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ < r² < 0.33 and 0.001 < r² < 0.2, respectively), indicating that size does not have an influence on these values. This suggests that our size-normalization procedure adequately removes the size influence (Fig. 2), and therefore, our MBW data represent a robust parameter reflecting test wall thickness and calcification intensity not influenced by test size (Table 3). Therefore, MBWs can be considered as a reliable calcification intensity proxy.

Based on all the above, from this point we will focus our discussion on the MBW_area to discuss the foraminifera shell weight variability on seasonal, interannual, and pre- to post-industrial Holocene timescales.

Figure 3 Mean MBW_area (µg) values across the years and for a composite year for G. bulloides, N. incompta, and G. truncatulinoides in the Planier sediment trap. Light-green shaded area represents the high productivity period in the study zone (Rigual-Hernández et al., 2012).

4.2 Seasonal variations of foraminifera calcification in the NW Mediterranean

MBW_area values were calculated for the three species to illustrate the seasonal variability of these parameters (Figs. 3 and 4). Samples have been assigned to their corresponding month according to the mean cup sampling date.

The mean MBW_area for the three species in the Planier sediment trap are 4.29 µg (± 0.45 µg for G. bulloides, 4.04 µg (± 0.4 µg) for N. incompta, and 23.25 µg (± 6.2 µg) for G. truncatulinoides. The seasonal variations in shell calcification differ according to the species.

In the case of G. bulloides, maximum annual calcification values are reached during winter and early spring: 5.05 in December and 5.02 µg in March. January displays the highest mean value: 4.78 µg. Minimum values are reached during summer: 3.72 µg in June, which is also the month that exhibits the lowest mean MBW_area. Overall, there is a 1 µg seasonal difference in calcification between maximum and minimum values, which corresponds to a 24.5 % change in the mean MBW_area value. Mean seasonal standard deviation is 0.47 µg.

N. incompta shows a maximum in calcification in late spring to mid-summer: a maximum value of 5 µg is reached in May, while July is the month that displays the highest mean value (5 µg). Lowest values are reached in early spring: 2.96 µg in March, while January displays the lowest mean value (3.85 µg). Thus, the annual mean seasonal amplitude is 1.15 µg, which translates into a 28 % seasonal MBW_area variability. Standard deviation is 0.28 µg.

Finally, G. truncatulinoides displays a seasonal maximum MBW_area value in late summer–autumn, with a maximum reached in October: 35.07 µg, while November is the month that shows the highest mean MBW_area value (32.85 µg). The lowest value is reached in March: 13.57 µg, and April is the month that shows the lowest mean value: 18.45 µg. Seasonal MBW_area difference is 14.3 µg, a 60 % variability. Mean typical seasonal deviation is 3.7 µg.

4.3 Interannual MBW_area trends

Trends throughout the 12-year record are represented in Figs. 3 and 5. In order to obtain representative data for each year, maximize data availability of each species, and avoid the impact of months with insufficient specimens on the interannual trends, only MBWs from the productive period (January to May) of each year analysed were included.

G. bulloides MBW_area showed a slight decrease from 1994 to 2000 and a slight increase from 2000 to 2006. Over the studied interval, the lowest value is reached in the year 2000 and the highest in 2004. Lowest mean annual values were reached during years 2000 and 2005 (3.9 and 3.85 µg, respectively).

On the other hand, N. incompta MBW_area showed a slight calcification reduction with the highest variability in recent years. Both maximum and minimum values are displayed in recent years: 2004 and 2005, respectively. Mean yearly MBW_area values reach a maximum in 2003 (4.4 µg) and a minimum in 2005 (3.2 µg).

Finally, G. truncatulinoides MBW_area displayed a different pattern, with an overall steep calcification increase throughout the record. Minimum calcification values are observed in 1996, which is also the year with the lowest mean MBW_area (16.5 µg) observed. Maximum value is displayed in 2003, and its mean value is also the highest of the record (26.1 µg).

All environmental parameters showed variations across the years. Sea surface temperatures (SSTs) displayed a slight but constant decrease over the years, while salinity showed a slight increase, mainly since 2002. From late 2000 until late 2002, phosphate and nitrate concentrations were exceptionally low (Fig. 5). This feature has already been described in the Gulf of Lions (Meier et al., 2014). Between the two periods for which direct in situ carbonate system parameters measurements were available, 1998 to 2000 and 2003 to 2005 (Fig. 5), CO${}_{\mathrm{3}}^{\mathrm{2}-}$ dropped by 10–15 µmol kg⁻¹, dissolved organic carbon (DIC) increased by 40 to 60 µmol kg⁻¹, leading to a pH decrease of 0.02 to 0.025.

4.4 Sediment trap, core top, and sediment core MBW patterns

Foraminifera weights analysed in core tops and sediment cores from the NW part of the Mediterranean (Fig. 6) and radiocarbon dating allowed a further insight on foraminifera calcification during the Holocene.

Flux-weighted MBWs (see Sect. 3.4) from Planier sediment trap for the three planktic species were 4.1 µg for G. bulloides, 3.9 µg for N. incompta, and 22.3 µg for G. truncatulinoides (Fig. 6).

Data from Planier core top showed higher mean MBW_area values: 5.3, 4.65, and 35.4 µg. ¹⁴C dating carried out in this core top was out of the calibration range (see Sect. 3.7 for more details), implying that this sample could be considered post-industrial. Compared to the flux-weighted MBWs from the sediment trap, G. bulloides weight has been reduced by 1.2 µg, N. incompta by 0.75 µg, and G. truncatulinoides by 12–13 µg.

Located west of Planier site, Lacaze–Duthiers sediment core mean MBWs were: 4.99 µg for G. bulloides, 4.14 µg for N. incompta, and 32.9 µg for G. truncatulinoides. ¹⁴C analysis displayed a post-industrial age (see Sect. 3.7) for this sample, and corresponding MBWs from this sample for G. bulloides, N. incompta, and G. truncatulinoides, respectively, were: 4.7, 4.3, and 34 µg. Overall, compared to the data from the sediment trap, this corresponds to a 0.6 µg weight loss for G. bulloides, 0.4 µg for N. incompta, and 12.2 µg for G. truncatulinoides.

Finally, in the Gulf of Minorca, northwest of the Planier site, Minorca sediment core mean MBWs were: 5.4 µg for G. bulloides, 4.5 µg for N. incompta, and 36.3 µg for G. truncatulinoides (Fig. 6). ¹⁴C dating on this core top was carried out on an intermediate depth (see Sect. 3.7) due to the lack of availability of enough specimens in the core top, and it displayed a date of 1560 calendar years BP (Table 1). Corresponding MBWs for this sample were 5.4, 4.9, and 38.2 µg for the three species. Therefore, the weight reduction compared to the sediment trap flux-weighted MBWs are 1.3 µg for G. bulloides, 1 µg for N. incompta, and finally 16 µg for G. truncatulinoides.

Table 4 Mann–Whitney variance test results between the MBW_area of the different sites for the three species. Significant values (p < 0.05) are set in bold.

A Mann–Whitney variance test (see Sect. 3.6) was carried out in order to analyse the variance between the different MBW_area data sets (Table 4) from the different sites. MBW_area data from the sediment trap appeared to have a significantly different variance compared to the MBW_area from Menorca sediment core for the three species (3.13 $\times {\mathrm{10}}^{-\mathrm{15}}$ < p < 2.59 $\times {\mathrm{10}}^{-\mathrm{5}}$), however, differences between the sediment trap data and with the Lacaze–Duthiers sediment core were only significant for G. bulloides and G. truncatulinoides (p=0.003). Concerning differences between the Planier sediment trap and the underlying core top, no significant differences were observed for any of the species (0.11 < p < 0.438), most likely due to the small number of samples from the latter site: only two samples were available. Note that the differences between the sediment cores MBW_area data sets differed according to the site and species. In the case of G. bulloides, no significant differences were observed between Planier core top, Lacaze–Duthiers sediment core, and Menorca sediment core. In the case of N. incompta and G. truncatulinoides, differences between Lacaze–Duthiers and Menorca sediment core are significant (0.013 < p < 0.03), although on lower orders of magnitude compared to the differences between the sediment trap and sediment cores data sets (Table 4). This demonstrates that the difference between the sediment trap MBW_area data set and the seabed sediments MBW_area data sets was greater than the difference between the different seabed MBW_area data sets.

5 Discussion

5.1 Seasonal controls on planktic foraminifera shell calcification in the NW Mediterranean

As described in Sect. 4.2, the seasonal variability of MBW_area displays important differences across the three species analysed. The different seasonal pattern in MBW_area is reflected by the lack of correlation between the seasonal patterns of MBW_area of the different species, i.e. $r=-\mathrm{0.23}$ (p>0.05) between G. bulloides and N. incompta and r=0.16 (p>0.05) between G. bulloides and G. truncatulinoides. Only the seasonality of N. incompta MBW_area and G. truncatulinoides MBW_area share some similarities, as reflected in the significant and positive correlation (r=0.66; p < 0.05). In order to examine the main controls on foraminifera seasonal calcification in the Gulf of Lions, here we compare the seasonal variability of planktic foraminifera calcification with foraminifera fluxes previously estimated for the Planier sediment trap (Rigual-Hernández et al., 2012) satellite data for the studied site and a suite of environmental parameters measured at the DYFAMED site (see Sect. 3.4). Furthermore, GAM have been generated for all three species (see Figs. S3, S4, and S5 in the Supplement) and the environmental parameters considered here in order to give a further insight on the potential factors controlling the MBW_area. These models showed that G. bulloides and G. truncatulinoides seasonal calcification trends are significant (p=0.05 and p=2.4 $\times {\mathrm{10}}^{-\mathrm{5}}$, respectively). On the other hand, N. incompta seasonal trend does not appear to be significant (p=0.14).

Table 5 Correlation matrix of seasonal (monthly) test weights and the environmental parameters from Planier (sediment trap and satellite data) and DYFAMED site (see Sect. 3.4). Significant correlations (p < 0.05) are set in bold.

Here, we first approach seasonal shell calcification by considering the optimum growth conditions (OGCs). Previous studied have defined these conditions on a wide variety of ways: abundance of foraminifera, the chlorophyll-a concentration, and even nutrients concentration (de Villiers, 2004; Schiebel et al., 2001; Schiebel and Hemleben, 2017). Therefore, we aim to explore the impact of these parameters as OGC on the shell calcification.

Among all the environmental parameters, de Villiers (2004) suggested that shell calcification, and therefore MBWs, is primarily controlled by the OGC that can be defined as the most suitable environmental conditions for the development of a given planktic foraminifera species. Based on the latter study, it could be expected that favourable environmental conditions for foraminifera growth would lead to both greater shell fluxes and enhanced shell calcification. Our correlation analysis shows that only G. truncatulinoides displays a significant (and negative) correlation with its flux ($r=-\mathrm{0.66}$; p>0.05). GAM results (see Figs. S3, S4, and S5 in the Supplement) support these observations, with shell flux showing a stronger effect on the calcification for G. truncatulinoides compared to the other two species fluxes.

According to the OGC theory, species calcification patterns vary according to the species and their fluxes. Interestingly, G. truncatulinoides calcification correlates negatively and significantly with all three species fluxes, a pattern opposite to what the OGC theory predicts (de Villiers, 2004), i.e. optimum ecological niche is associated with enhanced calcification. Thus, a possible explanation reconciling our observations with the OGC theory may be that G. truncatulinoides tends to prioritize energy allocation toward growth and reproduction at the price of reduced calcification. It is thought to reproduce once a year in winter in subtropical waters, and it has been speculated that nutrient availability and the lack of predation could explain this strategy. During this interval, the other major species display low abundances in the water column, which could allow G. truncatulinoides to reproduce due to the lack of competition. N. incompta calcification displays a similar pattern, a negative correlation with all three species but with a lower level of significance. It's MBW_area correlates negatively and significantly (p < 0.05) with G. bulloides flux, but its fluxes correlate positively and significantly with the latter species fluxes (see Table S2 in the Supplement). This is interesting, as it may highlight interspecific relations. First, this could lead to the assumption that when G. bulloides dominates the assemblages, N. incompta also displays a high abundance (Rigual-Hernández et al., 2012). Then, it could show that when conditions are favorable, N. incompta reproduces at a higher rate at the price of thinner shells (Table 1). This agrees with the N. incompta life cycle, which is known to be outnumbered by opportunistic species when nutrient supply is high (Schiebel et al., 2002) but dominate the assemblages when stratified waters are set. Therefore, when conditions are favourable or when in cohabitation with opportunistic species, it could focus on its reproduction. Note that G. truncatulinoides and N. incompta MBW_area correlate positively and significantly (p < 0.05), showing a similar calcification pattern on a seasonal scale.

An alternative proxy for OGC that may be considered is chlorophyll a concentration. Chlorophyll is considered an indicator of the algal biomass concentration, which is known to represent a large part of some foraminifera species diet (Schiebel and Hemleben, 2017). However, our data only showed a significant correlation of chlorophyll a with G. truncatulinoides calcification. A stronger trend would be expected under the OGC theory for G. bulloides, as algae are a vital part of its diet (Hemleben et al., 1989; Schiebel and Hemleben, 2017). This lack of correlation between G. bulloides and chlorophyll a has already been described (Weinkauf et al., 2016). We speculate that G. bulloides may preferentially feed on certain groups of phytoplankton whose changes in seasonal abundance in the photic zone do not necessarily follow the seasonal pattern of total chlorophyll concentration (Marty et al., 2002). Also, note that the chlorophyll a data presented here only represent the conditions in the surface layer. GAM results further support these observations (see Fig. S5 in the Supplement), with chlorophyll a showing a significant impact on G. truncatulinoides calcification. This observation indicates that optimum calcification conditions for G. truncatulinoides are reached at times of minimum annual algal biomass concentration in the photic zone. It is possible that, due to its deeper habitat (Schiebel and Hemleben, 2017), G. truncatulinoides feeds on phytoplankton dwelling in subsurface levels of the water column. In fact, a deep chlorophyll maximum is known to develop during large part of the year in the northwestern Mediterranean (Estrada et al., 1993), but its presence is not detected by satellites. This interpretation is in agreement with earlier work by Pujol and Vergnaud Grazzini (1995) who found peak abundances of this species during the summer below the thermocline.

Figure 4 G. bulloides, N. incompta, and G. truncatulinoides seasonal mean monthly MBW_area variations compared with Planier environmental data and the resampled seasonal signal of environmental parameters from the DYFAMED site across a composite year.

Previous studies have described that, in those settings where foraminifera abundance covaries with nutrient concentrations, then nutrients are probably a better OGC proxy than chlorophyll concentrations (Schiebel et al., 2001). In turn, the correlation of nutrients (nitrates and phosphates) with fluxes were positive for all three species, although only significative (p < 0.05) for G. truncatulinoides abundance (r=0.58 and 0.59 for nitrates and phosphates, respectively). Although here we have first described the OGC as species fluxes and then as the chlorophyll a concentration, it is important to remember that the niche and favourable conditions meant to be described by the OGC for each species are multidimensional.

Note that nitrate and phosphate concentration variations were closely linked to each other (r=0.876; p < 0.05), making it difficult to determine if the resulting effect on foraminifera calcification is due to the effect of a single driver or to the combination of both. Our work shows that nutrient concentrations (both nitrates and phosphates) do not correlate significantly with any of the three species MBW studied, and this observation is supported by the GAM results which do not show any significant effect of nutrients concentrations on the calcification.

Previous studies have suggested that salinity may have an influence on foraminifera calcification (Zarkogiannis et al., 2022). However, our data suggest that the role of salinity on calcification in our study region is unlikely, since its seasonal amplitude is tiny (0.1 PSU; Fig. 4). This idea is supported by the lack of correlation between salinity and MBW_area for the three species studied (Table 5) and the GAM results.

Temperature (sea surface temperature) has been described as a major factor that controls the size (Schmidt et al., 2004) and porosity (Burke et al., 2018) of planktic foraminifera, therefore it could represent a major control factor on shell calcification in the NW Mediterranean. In particular, calcification could be positively linked to temperature through different mechanisms: (i) warmer temperatures have been shown to increase enzymatic activity and therefore enhanced growth and calcification rates (Spero et al., 1991); and (ii) Lombard et al. (2011) stated that higher temperatures could also increase feeding and ingestion rates, but it remains unclear if this could result in a calcification rate increase. Our data revealed that SST correlates positively and significantly with N. incompta and G. truncatulinoides calcification (r=0.69 and 0.67, respectively; p < 0.05). GAM results also displayed a positive and the most significative effect of the SST on these two species. These findings highlight that SSTs are one of the main factors affecting N. incompta and G. truncatulinoides calcification among the parameters considered here. Finally, in addition to having an impact on the size and calcification of the planktic foraminifera, temperature is well known as a major control of the carbonate system, due to the increased solubility of atmospheric CO₂ at lower temperatures, and therefore it could have an indirect effect on foraminifera calcification by affecting the carbonate system.

Data for the carbonate system were only available for years 1998 to 2000 and 2003 to 2005, and therefore gaps comprised in these years were filled with estimates using the CO2SYS macro (see Sect. 3.6 for more details). However, note that the data available for these parameters were relatively smaller compared to the other parameters and may have prevented detection of other significant relationships. The relationship between CO${}_{\mathrm{3}}^{\mathrm{2}-}$ and MBW has been described in previous studies (Barker and Elderfield, 2002; Marshall et al., 2013), and the bulk of evidence indicates that foraminifera MBWs to be positively linked with CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations (Aldridge et al., 2012; Osborne et al., 2016). However, it appeared that planktic foraminifera response to CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentration was not uniform and varied across species (Beer et al., 2010b; Lombard et al., 2010). The trends between carbonate system parameters and MBWs were similar to those observed when comparing MBWs with temperature, highlighting the covariations between these two parameters (Fig. 4). Our data showed that CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations were only significantly correlated with G. truncatulinoides MBW_area (r=0.73; p < 0.05), implying that carbonate availability may represent a key control on this species in the northwestern Mediterranean. On the other hand, GAM result (see Figs. S3, S4, and S5 in the Supplement) did not show a significant impact of any carbonate system parameters for any of the three species calcification. As stated previously, the lack of data could have prevented the detection of further trends, but considering the seasonal patterns of carbonate system parameters, a potential role of the CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentration could be expected.

In summary, seasonal correlations, trends, and GAM showed that the environmental parameters that displayed the highest correlation with MBW_area vary according to the species. G. bulloides calcification appeared to be affected mainly by the OGC and interspecific relations. N. incompta calcification showed to be mainly positively linked to the SST. Finally, G. truncatulinoides calcification was positively linked with the SST and potentially CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentration, while OGC displayed a negative effect on its MBW_area. The combined effect of these parameters seems to control foraminifera calcification in the Gulf of Lions; however, it should be considered that covariation between these parameters is strong, and therefore it is difficult to isolate the effect of a single parameter. Moreover, it is likely that the ecology and life cycle of the species could also be a major factor affecting the response of the species calcification to the environmental parameters variations. Our results are in agreement with earlier studies that stated that OGC (de Villiers, 2004), SST, and CO${}_{\mathrm{3}}^{\mathrm{2}-}$ (de Villiers, 2004; Marshall et al., 2013; Osborne et al., 2016) concentrations are the main factors that impact calcification in planktic foraminifera, while the calcification response to those parameters is species-specific, which is in agreement with the work of Weinkauf et al. (2016).

5.2 Interannual trends in planktic foraminifera calcification

As stated previously, the Mediterranean Sea is a sensitive zone to atmospheric CO₂ accumulation (Ziveri, 2012) and is experiencing ongoing ocean acidification. On an interannual timescale, different studies (Beer et al., 2010b; Osborne et al., 2016) have shown that sea surface warming and carbonate system parameters are the most likely parameters to control calcification on key calcifying phytoplankton species such as the coccolithophore Emiliania huxleyi organisms (Meier et al., 2014). However, data sets from sediment traps that cover a wide span of years and in which foraminifera weights have been analysed are rare (Kiss et al., 2021), therefore it is difficult to place our results in a more global context. Our GAM results (see Figs. S3 and S4 in the Supplement) showed that both G. bulloides and N. incompta interannual patterns were non-significant. This is not surprising as the calcification trends for these two species did not display a clear and marked variation over the years, excepting a small mean calcification reduction (Fig. 3) and minimum calcification values in 2004 and 2005 (Fig. 3 and Figs. S3 and S4 in the Supplement).

Figure 5 Interannual mean MBW_area (µg) values for G. bulloides, N. incompta, and G. truncatulinoides from the high productivity period (see Sect. 2) and Planier and DYFAMED environmental data variations across the record. Black lines represent the trends from the MBW_area and resampled data. DIC represents dissolved inorganic carbon.

Notably, the trend in G. truncatulinoides is opposed to the previous two species and shows a steady and steep increase throughout our record. Over the analyzed time span, its MBW increased around 20 % (equivalent to an increase of ∼ 5 µg). According to the GAM results, the interannual calcification trend for this species is significant (see Fig. S5 in the Supplement). If this calcification increase continues on current trends, then the average MBW of G. truncatulinoides will double by 2024. Analysis of present G. truncatulinoides populations is urgently needed to assess if the observed trend held true during the last 2 decades. It is important to note that while G. truncatulinoides seems to exhibit a positive correlation with CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentration on a seasonal scale, no clear correlation was found with the interannual changes of CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentration. This feature is also supported by the GAM results. A similar enhancement in shell calcification has been described in the Balearic Sea for G. truncatulinoides in high-resolution sediment cores (Pallacks et al., 2020) but also in Globorotalia inflata. Taken together, our observations and the study mentioned above suggest that deep dwellers are unaffected by the recent ocean acidification and changes in the carbonate system and that the recent change in one or several environmental drivers may be stimulating the calcification of these species.

Here, we theorize that the interannual patterns presented in Figs. 3 and 5 mainly reflect the seasonal changes in the regional oceanographic setting. As described previously (see Sect. 2 for more details), the Gulf of Lions is influenced by a strong seasonality. The recent SST decrease could be linked to an enhancement in water mixing, as cold and deep salty water reach up to the surface. This mechanism would be less intense during years 2000 to 2002, corresponding to a SST increase along with a salinity decrease and absolute minimums in nutrients concentrations (Fig. 5); as water stratificates, these are consumed by primary production. Finally, in recent years, water mixing seems to be reactivated, as SST keeps decreasing and nutrients concentrations increase again. This mechanism also affects the carbonate system parameters, as water mixing brings to the surface deeper DIC-enriched waters, coupled with a [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] reduction. Our data show that alkalinity patterns display similar tendencies to DIC; however, until the second time span covered by carbonate system data, alkalinity variations are proportionally higher than DIC variations (see Supplement), suggesting a water mixing phenomenon. On the other hand, DIC variations turn to be higher than alkalinity variations from 2003 to 2005, suggesting an additional effect of carbon inputs on the carbonate system not reflected in the alkalinity data.

Note that SSTs, despite showing a positive and significant correlation with N. incompta and G. truncatulinoides on a seasonal scale and the GAM showing a positive and significant effect on the calcification, did not follow the same pattern as the latter species. This observation implies that other mechanisms or parameters than the ones considered here may be affecting the MBW_area in recent years.

Figure 6 MBW_area in µg and area in µm² comparison in the sediment trap (PLA ST), Planier core top (PLA CT), and both Lacaze–Duthiers (LCD SC) and Minorca sediment core (MR 3.1.A). Black dots represent data from the sediment trap, while lighter colours represent data from the different seabed sediments. Red lines represent the flux-weighted values from the sediment trap. Note that the age provided represents the dating (see Sect. 3.7) results carried out in specific samples of each site (see Sect. 4.3).

5.3 Holocene core top data comparison

The comparison of the well-preserved assemblages of planktic foraminifera in the pre-industrial and industrial Holocene-aged surface sediments with those collected by a long-sediment trap record offers a unique opportunity to assess the impact of recent environmental change on the calcification of calcareous zooplankton in the Mediterranean Sea (Fig. 6). However, when comparing data from sediment traps and seabed sediments, the possible role of calcite dissolution must be taken into account.

Calcite dissolution in the water column and/or on the sea floor could be invoked as a source of variability between the sediment trap and surface sediment data sets (e.g. Dittert et al., 1999). Therefore, in order to obtain meaningful interpretations from our data sets, it is important to assess the possible role of dissolution in the preservation of planktic foraminifera shells. Several lines of evidence suggest, however, that calcite preservation does not represent an important source of bias in our study area. Firstly, the Mediterranean Sea is supersaturated with respect to calcite (Millero et al., 1979) and the location of all the analyzed samples is much shallower than the location of the calcite saturation horizon (Álvarez et al., 2014), therefore, calcite dissolution seems unlikely (Schneider et al., 2007). Secondly, several sediment trap studies have documented that calcareous plankton experience negligible dissolution in their transit from the surface ocean to the sea floor (Beaufort et al., 2007; Moy et al., 2009; Rigual-Hernández et al., 2020). Thirdly, SEM and microscopic observations of all three species in samples from both the sediment traps and sediment cores showed no sign of dissolution and foraminifera were well preserved (see Fig. S6 in the Supplement). These arguments suggest that calcite dissolution does not represent an important control in the weight of the planktic foraminifera shells in the analysed samples. However, it has been documented that when dissolution takes place, the thinnest shells are affected first (Berger, 1970), while the heaviest and more calcified specimens remain. In our study, the specimens from the sediment trap were lighter that the ones from the sediment cores, which corresponds to the previous statement. This is important to acknowledge, as the individuals from the seabed sediment could only represent more calcified and solid individuals, while the most fragile and less calcified individuals may not have been preserved. Therefore dissolution cannot be completely ruled out here as a possible source of variability between the surface sediment and sediment trap data sets.

Overall, the lower shell weights of the foraminifera collected by the traps suggests that the three planktic foraminifera species may have experienced a reduction in their calcification since pre-industrial times to post-industrial and recent Holocene. While the shell weight of each species measured in the sediments shows some variability across seabed sediments (Fig. 6), our data suggest an overall reduction of 18 %–24 % for G. bulloides, 9 %–18 % for N. incompta, and 32 %–40 % for G. truncatulinoides. It is important to note that the range of shell weight variability across core tops and sediment cores (4.5–6.7 and 0.37 µg typical deviation for G. bulloides, 3.8–4.9 and 0.23 µg typical deviation for N.incompta, and 29.5–40.9 and 2.6 µg typical deviation for G. truncatulinoides) is substantially lower than the difference with the sediment trap data (3–5 and 0.5 µg typical deviation for G. bulloides, 2.9–5.2 and 0.5 µg typical deviation for N. incompta, and 12–35 and 6 µg typical deviation for G. truncatulinoides), implying that the shell weight of recent foraminifera populations for the three species is lower than anywhere in the NW Mediterranean in the pre-industrial and post-industrial times. The source of the variability across core tops and sediment cores is most likely caused by the different age of the samples, ranging from 1560 cal. yr BP at Minorca mid-depth (Table 1) sample to post-industrial at Planier and Lacaze–Duthiers core tops, and the different environments associated to the location of each core top.

A non-parametric two-way Mann–Whitney test (see Sects. 3.6 and 4.3) showed that the sediment trap MBW_area data set was significantly different (p < 0.05) from MR 3.1.A and non-different from PLA CT for all three species studied here (Table 5).

Something to consider when comparing recent sediment trap data with pre-industrial Holocene data is the life cycle of the species. As all the species analysed presented a lighter weight in the sediment traps, the degree to which the different specimens responded vary. The greatest weight reductions were observed for G. truncatulinoides populations, while N. incompta populations exhibited the lowest weight loss.

Previous work stated that those species hosting photosynthetic algal symbionts exhibit a higher tolerance to environmental changes that may affect their calcification (Lombard et al., 2009). This is due to the fact that these symbionts can modify the sea water chemistry that is in a close range to the shell, allowing a calcification enhancement. Of the species studied here, none are known, with the possible exception of G. bulloides, to be symbiont bearing species, therefore, they are among the most vulnerable foraminifera species to any sea water chemistry change.

It has been described that some morphotypes of G. bulloides host bacterial endobionts in their cytoplasm (Bird et al., 2017). The later work showed that high amounts of Synechococcus, a cyanobacteria, were found in morphotype Id specimens of G. bulloides from the California coast. Although no such observations have been reported on morphotype Ib, the dominant G. bulloides morphotype in the Mediterranean Sea (Schiebel and Hemleben, 2017), this could be relevant, as bacterial photosynthetic activity would interact on the close range seawater chemistry by removing ¹²CO₂ and therefore impacting the ¹³C $/$ ¹²C ratios in the surrounding dissolved CO₂. The work by Moy et al. (2009) on the Southern Ocean showed a 30 %–35 % calcification reduction for G. bulloides during the industrial era. Our study shows that such a similar reduction in G. bulloides MBW_area (i.e. a mean 20 % taking into account the three sites studied) has also taken placed in the Mediterranean Sea.

Even though the species studied were different in Fox et al. (2020), and that shell thickness was analysed, the latter work showed a massive shell reduction for N. dutertrei (around 75 %) and a smaller reduction for G. ruber (around 20 %). Mean N. incompta weight reduction in this study is around 15 %; despite that the life cycles are different between these species, our results come in the same line.

Data for G. truncatulinoides calcification comparison between pre-industrial and post-industrial Holocene are scarce. One of the few available studies is one by Pallacks et al. (2020) on the western Mediterranean Sea using pre-industrial data and recent foraminifera weight data obtained from high-resolution core tops. G. truncatulinoides showed a 24 % weight reduction, which is a lower reduction than what is shown in our study (around 35 % MBW_area decrease) but showing a similar trend. Taken together, all these observations suggest that a decrease in major planktic foraminifera calcification is not only a regional feature but a global-scale process.

On a more regional scale, Hassoun et al. (2015) documented the ongoing changes in seawater carbonate speciation in the Mediterranean waters. In the latter work, the distributions of anthropogenic CO₂ showed that all Mediterranean water masses have already experienced ocean acidification. This effect was more pronounced in the intermediate to deep masses (300–500 and > 500 m, respectively) in the western basin, which translated into a minimum pH reduction of 0.1 in this part of the Mediterranean. As stated previously, over the years in which carbonate parameters were retrieved from the DYFAMED database, pH was reduced, DIC showed a marked increase, and [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] displayed a decrease. Taken together these observations and our data, it is possible that the observed changes in foraminifera calcification could have been partially driven by the ongoing ocean acidification in the Mediterranean.

Moreover, the largest calcification reduction is observed between the seabed sediments and the sediment trap; this means that the highest calcification reduction has taken place between post-industrial Holocene and recent Holocene (i.e. the reduction between LCD SC/PLA CT and PLA ST) (Fig. 6). This could be explained with the “Great Acceleration theory”. The Great Acceleration is a term used to describe the trends in CO₂ emissions and the associated temperature changes as consequences of the human impacts on the atmosphere since the 1950s (Head et al., 2022a, b).

However, other important changes in physical and chemical parameters co-occur with ocean acidification, and therefore should be also considered. Based on the seasonal and interannual patterns of the SST in the Gulf of Lions (Figs. 4 and 5), temperature trends could also be invoked as a likely parameter to affect calcification here. As shown by the correlations (Table 5) and the GAM results, SSTs are one of the most likely parameters to affect calcification on different timescales. However, on a pre-industrial to post-industrial timescale, the effect of the SST on the foraminifera calcification on longer timescales may be hard to evaluate due to the effect of the latter on the carbonate system parameters such as CO₂ and CO${}_{\mathrm{3}}^{\mathrm{2}}$ concentration in the water. But note that the Mediterranean is considered to be warming at a faster rate than the global average (Hassoun et al., 2015; Lazzari et al., 2014).

Calcification data from the sediment trap have been flux-weighted (see Sect. 3.4) in order to be compared with the sedimentary calcification data; therefore, these data could be affected by a change in the incoming foraminifera flux (de Moel et al., 2009). In this line, the Gulf of Lions presents a marked seasonality (see Sect. 2), and both the mass fluxes (Heussner et al., 2006) and foraminifera fluxes present strong seasonal variations. Parameters such as the North Atlantic Oscillation index, the river runoff, and the intensity of the seasonal water cascading process have been suggested to play a role on planktic foraminifera production and export (Rigual-Hernández et al., 2012). The later study shows that most of the species flux showed a yearly uni-modal distribution, but the flux values and distribution remained fairly constant over the years. This highlights that, in our study zone, a major change in the foraminifera flux affecting the flux-weighted calcification value is unlikely.

In summary, our results suggest that the interactive effect of rising ocean acidity, and enhanced SST (regionally amplified in the NW Mediterranean; Hassoun et al., 2022, 2015), represent the most likely responsible factors for the MBW differences between the pre-industrial and post-industrial to recent Holocene. However, the analysis of seasonal and interannual trends indicates that the influence of these parameters is species-specific and varies across the studied timescales. This implies that the controls of planktonic foraminifera are complex and that factors other than ocean acidification and warming are likely to also account for part of the variability observed between sediment trap and seabed sediments.

5.4 Influence of environmental variability on MBW_area across different timescales

Our results show that the influence of environmental parameters over the different timescales studied is not constant and depends on the species, the environmental driver, and the timescale.

In the case of G. bulloides, our data suggest that OGC and interspecific relationships seem to affect its MBW on a seasonal scale, and the carbonate system seems to play a major role while on an interannual and on pre/post-industrial timescales. N. incompta calcification seems affected by OGC, interspecific relationships, and SST on a seasonal scale, while on longer timescales, the carbonate system appears to play a preponderant role. Finally, G. truncatulinoides calcification seems positively correlated to carbonate system and SSTs and negatively with the OGC on a seasonal scale. However, these patterns seem to have an opposite effect on an interannual scale, as G. truncatulinoides calcification shows a clear increase, while carbonate system parameters become less and less favourable for calcification. In turn, on a pre/post-industrial Holocene timescale, its MBW_area seems to be affected by regional processes such as OA and warming.

Factors such as changes in the regional oceanographic processes (Cisneros et al., 2019; Durrieu de Madron et al., 2017) affect the physical and chemical properties of the water column and hence could impact the life cycle of the species studied here. Also, while G. bulloides can either present regular or encrusted forms, N. incompta and G. truncatulinoides are crust-forming species. In our study, G. bulloides individuals are mainly regular forms, but encrusted individuals were identified in both the sediment trap and seabed sediments. It is out of the scope of this work to focus on the effect of the crust on the species MBW, however, the study by Osborne et al. (2016) showed that encrusted G. bulloides individuals are around 20 %–30 % heavier than the regular ones.

6 Conclusions

The variability in shell calcification of three planktic foraminifera species (G. bulloides, N. incompta, and G. truncatulinoides) was studied in the northwestern Mediterranean Sea at different timescales using sediment trap and seabed samples. The analysis of 273 samples and more than 4000 individuals revealed the following:

• i.

  The sieve-based weight (SBW) method is not a reliable tool as calcification indicator due to the influence of morphometric parameters on foraminifera weight. The measured-based weight (MBW) technique, on the other hand, shows little to negligible influence of the morphometric parameters and therefore can be considered a reliable calcification proxy.

• ii.

  Analysis of the seasonal variability of planktic foraminifera calcification revealed important differences between species. G. bulloides exhibited peak calcification during winter, N. incompta during mid-summer, and G. truncatulinoides during late summer to autumn.

• iii.

  Interannual analysis suggests that G. bulloides and N. incompta did not display any significant pattern between 1994 and 2005; on the other hand, G. truncatulinoides displays a constant and steady calcification increase over recent years.

• iv.

  Sediment trap and seabed sediment data comparisons between pre-industrial, post-industrial, and recent Holocene assemblages showed that all three species experienced a calcification reduction. Modern G. bulloides, N. incompta, and G. truncatulinoides individuals were 18 %–24 %, 9 %–18 %, and 32 %–40 % less calcified, respectively.

• v.

  Finally, correlations with environmental parameters and GAM indicate that optimum growth conditions positively affect G. bulloides and negatively affect G. truncatulinoides calcification, respectively. Sea surface temperatures positively affect both N. incompta and G. truncatulinoides calcification. Finally, CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentration is also a likely parameter to positively influence planktic foraminifera calcification in the northwestern Mediterranean. However, calcification appeared to be species-specific and varies depending on the timescale studied. This may suggest that other parameters than the ones studied here may play a role in foraminifera calcification.

As planktic foraminifera represent roughly about 50 % of pelagic calcite production (Schiebel, 2002) in the world's oceans, and therefore are an important component of the marine carbon cycle, a reduction in the calcification of their shell could induce important changes in the future carbon cycle with feedbacks on climate. Our results call for increasing efforts in monitoring planktic foraminifera calcification in the Mediterranean in order to determine if the trends suggested by our data will be sustained over time.