Abstract
A mid- to late-Holocene paleohydrological reconstruction from the northeast Atlantic is proposed through the study of a high-resolution sedimentary record from the northern continental shelf of the Bay of Biscay (BoB). Three foraminiferal species dominate the assemblages with Rosalina globularis showing an overall decrease in absolute and relative abundances from ~7 to 0.4 cal. ka BP, whereas the opposite trend is observed for Cibicides refulgens and Lobatula. These long-term patterns are interpreted as a response to the overall cooling trend and/or the progressive deepening of the water column because of the relative sea-level (RSL) rise. Foraminiferal δ18O and grain-size analyses show a significant shift around 3.5–2.5 cal. ka BP toward a heavier isotopic signature and finer sediments. We mainly link this change to enhanced contribution of continental freshwaters and fine sediments after the near-stabilization of the RSL rise. By reducing coastal accommodation spaces, this led to a better channelization of river outflows and probably to the formation of the modern winter thermohaline front. Superimposed on these long-term patterns, our data highlight strong millennial-scale variability (1250-year peak). Such cyclicity is consistent with several records tracing changes in rainfall and storminess regimes in northern Europe, and the dynamics of the subpolar gyre (SPG). We suggest a millennial time-scale control of a NAO-like (North Atlantic Oscillation) climatic process modulating continental humidity (and the associated river discharges) and SPG dynamics through wind stress. Spectral analyses reveal an additional 500-year frequency peak implying a possible solar forcing.
Introduction
Holocene continental shelf deposits in the Bay of Biscay (BoB) have the potential to capture both continental and open marine signals (e.g. Durand et al., 2018; Garcia et al., 2013; Mojtahid et al., 2013; Naughton et al., 2007). The continental signal originates mainly from the Dordogne and Garonne Rivers (referred to in the following as the Gironde Rivers) and the Loire River (Figure 1a). Their runoff is conditioned by rainfall regimes over their catchment areas (Figure 1a). The marine signal is captured because the oceanic circulation in the BoB is under the influence of a branch of the North Atlantic Current (NAC). The modern variability of the Atlantic Meridional Overturning Circulation (AMOC) and the associated inflow from the NAC to higher latitudes are of great importance in modulating climate over Europe (e.g. Clark et al., 2002; Jackson et al., 2015; Staines-Urías et al., 2013; Thornalley et al., 2009). Particularly, the expansion/contraction of the North Atlantic subpolar gyre (SPG) modulates the exchange of heat, salinity, and freshwater with the Nordic seas (Born and Stocker, 2014; Häkkinen and Rhines, 2004), thus pre-conditioning the strength of the AMOC through the production of deep and intermediate water masses (e.g. Hátún et al., 2005; Sherwin et al., 2012; Yong-Qi and Lei, 2008). Along the French Atlantic coast, the slowing down of the relative sea-level (RSL) rise from middle to late Holocene (García-Artola et al., 2018; Goslin et al., 2015) led to the reduction of coastal accommodation spaces, which in turn influenced the export of freshwaters and sediments to the ocean (Clavé et al., 2001; Durand et al., 2016; Spencer et al., 1998). In addition, internal and external climate forcing influenced the runoff of the French rivers as well as the strength of the AMOC. It is fairly well known that rainfall over Europe is governed, besides local hydrometeorological control, by oscillating internal climatic forcing such as the North Atlantic Oscillation (NAO) (Hurrell, 1995; Hurrell et al., 2001; Labat, 2008). Several studies indicate that the hydrological variability of the main French watersheds (Seine, Loire, Garonne, and Rhône) has a high degree of coherence with the NAO (Chevalier et al., 2014; Milliman, 2001; Shorthouse and Arnell, 1997). Other studies are more doubtful about the dominant effect of the NAO on the hydrology of the French rivers and/or evoke probable connections with other climatic mechanisms such as the AMO (Atlantic Multidecadal Oscillation) and EAP (Eastern Atlantic Pattern) (e.g. Boé and Habets, 2014; Costoya et al., 2017). The study of Tréguer et al. (2014) reveals that the modern variability (AD 1998–2013 period) in sea surface temperature (SST) in coastal waters of Brittany is connected to the NAO while the EAP exerts a major influence on discharges of the adjacent rivers. Meanwhile, at annual, decadal to probably millennial time-scales, the NAO is also known as the dominant mode of climate variability controlling the intensity of winter westerlies across the North Atlantic (Hurrell, 1995; Hurrell et al., 2008; Fletcher et al., 2012; Thompson et al., 2003; Wanner et al., 2008). Typically, in the northern North Atlantic, stronger westerly winds (typical of a positive NAO state) increase the depth of the mixed-layer and strengthen the SPG (Sarafanov, 2009; Sarafanov et al., 2010). Superimposed on this internal climatic forcing, solar activity has also a major imprint on Holocene climate variability (e.g. Debret et al., 2009; Soon et al., 2014).

(a) A map illustrating the general surface circulation of the BoB (Koutsikopoulos and Le Cann, 1996) and the extent of the Loire and Gironde catchment areas. The locations of core CBT-CS11 (this study) and the other cited cores in the text are indicated: VK03-58/59bis (Naughton et al., 2007); KV14-bis (Durand et al., 2018); ProsecanKS10b (Mojtahid et al., 2013); PP10-07 (Brocheray et al., 2014). Map performed using the online EMODnet Portal for Bathymetry (http://www.emodnet-bathymetry.eu/); (b) zoom on the study site, and the location of the nearby cores St G (Duchemin et al., 2008) and vk16 (Estournès et al., 2012). ENAW: Eastern North Atlantic Water; EPC: European Poleward Current.
Here, we investigate Holocene variability as recorded in sediments from the mid-shelf mud belt ‘La Grande Vasière’ (south of the tip of the French Brittany; Figure 1), in response to global (e.g. AMOC dynamics, NAO-like processes, solar forcing) and regional forcing (e.g. sea-level rise, continental inputs). We use benthic foraminiferal assemblages and stable isotopes measured on their shells supplemented by grain-size x-ray fluorescence (XRF) analyses. Several studies based on a similar multi-proxy approach have been conducted successfully in various marginal environments to track Holocene variability (e.g. Bartels-Jónsdóttir et al., 2006; Mojtahid et al., 2015; Polyak et al., 2002; Xiang et al., 2008), but only very few have been performed in the BoB because of the generally low preservation of Holocene sediments (Durand et al., 2018; Garcia et al., 2013; Mojtahid et al., 2013). One of the main interests of our study is that the exceptional sedimentation rate (~70 cm ka‒1) characterizing our coring site (a depocenter of fine sediments; Dubrulle et al., 2007) permits studying the past ~7 kyr at a high time-resolution (~19–130 years), unprecedented in the BoB. This allows robust spectral analyses, which highlight for the first time in marine sediments from the BoB a clear millennial variability, linking our local record to larger scale mechanisms such as North Atlantic hydrological patterns.
Study area
Sedimentary characteristics
‘La Grande Vasière’ is a mid-shelf mud belt (McCave, 1972) extending over more than 225 km long and 40 km wide on the Armorican shelf between 80 and 110 m water depth (Lesueur et al., 2001) (Figure 1). Characterizing it as a mud belt is, however, not entirely satisfactory because it has been described as a thin and broken blanket of muddy sands in which the silt and clay content does not exceed 30% (Lesueur et al., 2001). The quaternary sequence of the Armorican continental shelf starts with a regional erosion surface inherited from river incisions of Eocene limestones during an early Pleistocene regression (~700–500 ka BP; Proust et al., 2001). This sequence is subdivided into (a) the ‘lower unit’ which consists of several meters of bioclastic coarse to medium sands deposited in a ~0–50 m water depth range, (b) the ‘upper unit’ which consists of several meters of coarse silts to fine gray sands deposited in a ~50–100 m water depth range (Andreieff et al., 1971), and (c) the ‘surficial layer’ which consists of few decimeters of muddy and bioturbated silts (Bourillet et al., 2006b). The two upper units, that concern our sediment core (Figure 2), form the mud belt ‘La Grande Vasière’ and were deposited over the past ~9 kyr BP (Bourillet et al., 2006a; Naughton et al., 2007) in a context of three major phases of RSL rise (Goslin et al., 2015): (a) a relatively rapid increase between ~8 and ~5.5 cal. ka BP (‒10 and ‒3 m RSL), (b) a near-stabilization between ~5.5 and ~2.5 cal. ka BP (around ‒3/‒2 m RSL), and (c) a slight rise in the most recent period. In a context of RSL high stand, such as since ~5.5 cal. ka BP, major sources of fine sediments (silts and clays) to ‘La Grande Vasière’ are fluvial supplies carried in suspension, essentially from the Loire and Gironde Rivers and, to a lesser extent, from the Adour, Vilaine, and Charente Rivers (Guillaud et al., 2008; Jouanneau et al., 1999). Today, the sandy fraction is mainly sourced from coastal erosion occurring mostly during winter storm events, whereas direct sandy fluvial contributions are considered negligible (Castaing and Jouanneau, 1987). In our core, the continuous sedimentation (14C chronological coherency) and the absence of erosion surfaces, especially in the first half of the record dominated by fine sands (Figure 2), may point out to either a non-negligible direct fluvial input at the time or to distal tempestite deposits characterized by fine sands and weakly developed erosional surfaces (Aigner, 1985).

Synthetic sedimentological log and age model of core CBT-CS11.
Surface circulation in the BoB
The central part of the BoB is influenced by the Eastern North Atlantic Water (ENAW) transported by the NAC following an anticyclonic circulation (Charria et al., 2013; Pingree, 1984; Pollard and Pu, 1985; Figure 1a). The circulation becomes cyclonic along the continental slope during autumn and winter with the European Poleward slope Current (EPC) transporting northward ENAW of subtropical origin (Garcia-Soto and Pingree, 2012) (Figure 1a). The EPC, which flows from the surface to 1000 m depth along the continental slope, turns equatorward in summer (Charria et al., 2013). The slope boundary flow appears unstable, particularly where there are abrupt changes in topography (Pingree and Le Cann, 1992a). This results in the intrusion of slope waters into the central BoB and the consequent formation of surface-intensified meso-scale structures (~100 km in diameter) called SWODDIES for Slope Water Oceanic eDDIES (Dickson and Hughes, 1981; Koutsikopoulos and Le Cann, 1996; Madelain and Kerut, 1978; Pingree and Le Cann, 1992b) (Figure 1a). Over the shelf, weak currents occur from April to September (<2.5 cm s−1) and a poleward circulation takes place from October to March with intense currents reaching 10–15 cm s−1 on the Aquitaine shelf (Charria et al., 2013). This so-called residual circulation is controlled by wind, tides, and water density (Figure 1a) (Pingree and Le Cann, 1989).
Hydrographical pattern of the continental shelf and surface productivity
The hydrographical structure of the continental shelf of the BoB is highly seasonal (e.g. Castaing et al., 1999; Vincent and Kurc, 1969). One of the main features is the formation of a winter thermohaline front along the ~100-m isobath that lasts for several months (Castaing et al., 1999). This front results from the outflow of freshwaters mainly from the Loire and Gironde estuaries (Lazure and Jégou, 1998). The Loire and Gironde catchment areas are mostly characterized by a simple hydrological regime of oceanic-rainfall type (Joly et al., 2010), leading to maximum runoff during winter and early spring (>1500 m3 s‒1) and minimum values during summer (<300 m3 s‒1) (hydro.eaufrance.fr; Castaing and Allen, 1981). In response to high river runoff, temperature and salinity gradients in the BoB continental shelf are present from the inner shelf to the 100-m isobath (from ~10 to 12°C and ~30 to 35.6 psu; Lazure et al., 2006; Lazure and Jégou, 1998). During summer, thermal stratification takes over with increasing temperatures (Puillat et al., 2004). Continental suspended matter derived from estuaries spreads out on the Armorican shelf by means of surface (i.e. turbide plume) and intermediate nepheloid layers (Castaing et al., 1979; Weber et al., 1991), which are transported by the northward residual circulation. Suspended particles may be trapped by the winter thermohaline front, and are deposited at ~100-m water depth as a mud belt (Lazure and Jégou, 1998). It is postulated that the position of this mud belt is linked to the long-term stable location of the thermohaline front that separates oceanic waters from coastal waters (Castaing et al., 1999). During summer, another front, the Ushant tidal front, takes place north of the study area in the Iroise Sea (Cadier et al., 2017). The dynamics of these fronts (i.e. strength of T-S gradient) concentrate nutrient inputs, which in turn stimulate the seasonal production of phytoplankton in the northern BoB (e.g. Harlay et al., 2011; Pingree and Garcia-Soto, 2014; Sharples et al., 2009; Wollast and Chou, 2001).
Material and methods
Coring method
Piston core CBT-CS11 (47°46.429N, 4°25.308W; 73 m water depth; 3.87 m long) was recovered about 4.5 km off the tip of south Brittany during the CABTEX cruise (June 2010, CALYPSO corer, R/V Pourquoi pas?, IFREMER; Dussud, 2010; Figure 1). A 6.26-m long core barrel was used with 1.68 tons of head weights to achieve a reasonable balance between the corer diameter (10 cm) and the corer weight (2.205 tons). The corer was lowered to 1-m distance from the sediment water interface before being released by a trigger. Before penetration, some of the supernatant water was sucked out and the corer over-penetrated the sediment (recovery ratio = 107%).
Chronology
A total of 23 14C dates, of which 21 were measured using French national ARTEMIS 14C AMS facilities and two dates measured at Poznán Radiocarbon Laboratory, were obtained (Table 1). For dating, we used mostly the gastropod Turritella sp. When not found, other mollusks or mixed benthic foraminifera were used (Table 1). Absolute 14C dates were corrected for the geographically nearest reservoir age of 324 ± 24 years (Mangerud et al., 2006). By using the Bchron package in R program (R Core Team, 2014), corrected AMS 14C dates were converted to calendar ages (cal. ka BP) using the IntCal13 calibration curve (Reimer et al., 2013) and an age-depth model was computed based on 18 14C dates, five dates being considered as outliers from the main trend (Figure 2). Most of the 14C outliers are present in the upper 50 cm of the core. This might be the result of (a) physical disturbance during extraction, (b) bioturbation, and/or (c) the different nature of the dating material (Table 1). Therefore, the upper 50 cm (~ the last 0.38 kyr) are not discussed in the following.
AMS radiocarbon ages of core CBT-CS11.
Corrected 14C ages were calibrated using the atmospheric calibration curve IntCal13 (Reimer et al., 2013).
Reservoir correction inferred from Mangerud et al. (2006).
14C age showing age inversion or incoherent results in comparison with the general tendency. SacA: code for ARTEMIS facilities; Poz: code for Poznan laboratory.
XRF and grain-size analyses
The bulk intensity of elements (Ca, Sr, Fe, Ti, Al, Si, Br, Cl) was measured using an Avaatech XRF core scanner at 1-cm resolution at IFREMER (France). XRF data were measured with a 10 s count time, by setting the voltage to 10 kV (no filter) and intensity to 600 mA, except for Sr (voltage 30 kV Pd thick filter; 1000 mA). Because XRF scanning measures element intensities on a volumetric basis (i.e. sensitive to variations in grain size and density), XRF raw counts were normalized by applying the total number of counts (TNC) excluding Rh (x-ray tube) following the method of Bahr et al. (2014). To allow more robust statistical analyses, XRF ratios were converted to log-ratios (Tjallingii et al., 2007; Weltje and Tjallingii, 2008). It is generally accepted that Ti, Fe, Al, and Si are related to the siliciclastic components of the sediment and vary with the terrigenous fraction (Jansen et al., 1998). Because sediments of ‘La Grande Vasière’ originate mainly from the Loire and Gironde Rivers and upper shelf environments, increased intensities of these elements are interpreted as increased supplies of siliciclastic material of continental origin (e.g. Chiessi et al., 2009; Haug et al., 2001; Nizou et al., 2011; Peterson et al., 2000). In marine sediments, the coarse sediment fraction is generally enriched in Ti, while Al is mostly associated with fine-particle clay minerals, so that the ratio Ti/Al changes can represent grain-size variations (e.g. Govin et al., 2012). In a siliciclastic margin, Ca and Sr often represent the carbonate content of the sediment (e.g. Govin et al., 2012; Rothwell et al., 2006), with Sr indicating the aragonite content (e.g. Bayon et al., 2007). The biophilic halogen bromine (Br) is used as a tracer of organic matter in sediments, and the ratio Br/Cl permits to isolate organic Br from salt-bounded Br (Ziegler et al., 2008).
CaCO3-free grain-size analyses were performed for 69 samples with a sample resolution of about 5 cm (~100-year resolution) using a laser diffraction particle size analyzer (Malvern™ Mastersizer 3000; LPG). Sediment decarbonatation was performed with 1 M hydrochloric acid.
Benthic foraminiferal assemblages
Fifty-five sediment samples (on average every 7 cm; ~130 years resolution) were analyzed for their benthic foraminiferal content. Samples were dried and weighed to obtain dry weight. Samples were washed through 150 μm and 63 µm sieves. Although small sized and/or elongated adult species may compose a large part of the foraminiferal assemblage, the 63–150 µm fraction was not considered in this study. However, the inspection of some samples from the small fraction revealed the presence of mainly juveniles of the species found in the larger fraction. Because of the high foraminiferal abundances, samples were split into subsamples (aliquots), using an Otto Microsplitter. More than ~200 specimens were picked out from a single aliquot under a binocular microscope, stored in separate Chapman slides (for each sample), and identified at a species level. In order to correct for the taphonomical loss of some agglutinated taxa because of their rapid disintegration, the few specimens of the only non-fossilizing species Eggerella scabra (Fontanier et al., 2003) were removed before calculating absolute and relative abundances. Diversity (Shannon index H) was calculated using the software PAST (PAleontological STatistics; Version 2.14; Hammer et al., 2001) following the formula, with = the proportion of individuals belonging to the ith species.
Benthic foraminiferal oxygen and carbon stable isotopes
Stable oxygen and carbon isotopes were measured on the epibenthic species Cibicides refulgens in 113 sediment samples (~3 cm spacing; ~64 years resolution). For each sediment level, about 15 specimens of C. refulgens (150–250 μm) were hand-picked, cleaned in a methanol ultrasonic bath for a few seconds, then roasted under vacuum at 380°C for 45 min to remove organic matter, prior to isotopic analyses. The δ18O and δ13C (‰VPDB) were measured at the PSO (IUEM, Brest) using the IRMS platform: a Delta V mass spectrometer coupled with a GasBench II preparation line for benthic species (~1 test). The external reproducibility (1σ) of an internal standard calibrated with NBS19 is ±0.05‰ and 0.15‰ for δ13C and δ18O, respectively.
Statistical and spectral analyses
We performed a principal components analysis (PCA) and a Linear r (Pearson) multiple correlation analysis on all our proxies using PAST software package (Hammer et al., 2001). Because our variables have different units, we used the correlation matrix (i.e. normalized variance–covariance) for PCA calculation (Hammer et al., 2001). As all proxies have different sampling resolution, data were resampled (i.e. interpolated between data points) at continuous and even spaced time-resolution according to the minimum time-resolution obtained by the proxy with the lowest resolution (i.e. H = 30 years).
In order to explore the frequency in the temporal variability of selected parameters (δ18O, δ13C, Shannon H, and XRF-Ti, Br/Cl, Sr/Ti), a Continuous Wavelet Transform (CWT) (Morlet wavelet) and REDFIT analyses were performed using the PAST software package (Hammer et al., 2001). Before spectral analyses, data were resampled (i.e. interpolated between data points) at continuous and even spaced time-resolution (according to the minimum time-resolution; XRF data = 5 years, stable isotopes = 10 years, and H = 30 years).
Results
Grain-size and XRF analyses
Grain-size distribution shows a sharp drop in the median particle size (D50) changing from coarser (~75 µm) to finer sediments (~55 µm) between ~3.5 and 2.5 cal. ka BP (Figure 3a). XRF-Ti/Al shows generally the same pattern with higher values before ~3.5 cal. ka BP (Figure 3b). XRF-Ti, which is positively correlated to Fe, Si, and Al (Figure 3c), shows a great variability along the record oscillating between 0.006 and 0.011 cps/TNC (Figure 3d). XRF-Br/Cl values follow the same trend as Ti and oscillate between 0.10 and 0.40 cps/TNC (Figure 3e). XRF-Sr and Ca increase until ~5 cal. ka BP and remain generally stable after (Figure 3f). The ratio Sr/Ti varies between 80 and 140 cps/TNC (Figure 3g) and follows a similar pattern as XRF-Sr and Ca but with amplified oscillation.

(a) CaCO3-free grain size D50 (µm); (b–g) XRF data. The tick lines in XRF data represent the 3-point moving average. The dashed horizontal lines are placed on the average value. LIA: ‘Little Ice Age’; MCA: Medieval Climate Anomaly; RWP: Roman Warm Period; HTM: Holocene Thermal Maximum.
Foraminiferal assemblages and stable isotopes
Absolute foraminiferal abundances range between ~134 and 631 ind g‒1 of dry sediment, and generally increase along the record (Figure 4b). Shannon diversity (H) values show also an overall increase and range between 2.2 and 3 (Figure 4c). Percentages of the three major species (>10% in at least one sample) are presented in Figure 4d and e. Rosalina globularis shows a general decrease from 35% in the early part of the record to less than 10% at ~0.38 cal. ka BP. Lobatula lobatula and C. refulgens show the opposite trend with increasing values over time from ~10% to 20%. The percentages of the minor species (5–10% in at least one sample) are shown in Figure S1 (supplemental material, available online).

(a) Summer insolation at 47°N and a simplified Holocene reconstruction of sea level in Brittany (Goslin et al., 2015); (b) total benthic absolute densities (ind g‒1); (c) diversity index (Shannon H); (d–e) relative abundances (%) of the three dominant benthic foraminiferal species; (f–g) oxygen and carbon stable isotopes measured on the benthic species Cibicides refulgens. The tick curves represent the 3-point moving average. The dashed horizontal lines are placed on the average values. LIA: ‘Little Ice Age’; MCA: Medieval Climate Anomaly; RWP: Roman Warm Period; HTM: Holocene Thermal Maximum.
The stable oxygen isotope values of C. refulgens (δ18O C. refulgens ) show a great variability with a major change around 3.5–2.5 cal. ka BP (Figure 4f). Before ~3.5–2.5 cal. ka BP, δ18O C. refulgens oscillates between ~‒0.25‰ and 0.5‰, whereas after ~3.5–2.5 cal. ka BP, values generally range between 0.4‰ and 1‰. The stable carbon isotope values of C. refulgens (δ13C C. refulgens ) show also a great variability, oscillating between ~1.1‰ and 1.5‰ (Figure 4g).
Discussion
Evolution of benthic foraminiferal communities from middle to late Holocene
The downcore foraminiferal composition is largely dominated by epiphytic/epibenthic calcareous species, mostly R. globularis, L. lobatula, and C. refulgens (Langer, 1988, 1993) (Figure 5). Coastal areas in Brittany are known for their abundant and diverse marine vegetation extending from supralittoral to circalittoral depths (e.g. Golléty et al., 2010; Hart, 2016; Stiger-Pouvreau et al., 2014). It is well known that these latter provide abundant niches worldwide for motile or (temporarily) fixed epiphytic foraminifera (Barras et al., 2014; Hart, 2016; Johnson and Scheibling, 1987; Langer, 1993). The steep inner to mid-shelf (0–100 m) bathymetric profile (Figure 1b) makes our study area close to the upper coastal habitats, explaining the large abundance of epiphytes and the difference with the nearby site G located more distantly from the continent (Figure 1b; Duchemin et al., 2008). On the larger southern Brittany inner shelf (Bay of Etel; Figure 1b, vk16), a very similar Holocene fossil assemblage dominated by epiphytes was found in a core retrieved at ~30 m water depth (Estournès et al., 2012).

Statistical analyses: (a) PCA analysis based on correlation between proxies; (b) graph representing the results of the linear multiple correlations between proxies. Colors refer to Pearson correlation r, with blue/red referring to positive/negative correlation. The size and the color intensity of the round symbols are proportional to the strength of the correlation. Non-significant correlations (p > 0.05) are not represented.
Along the record, the evolution of the three main species follows two main long-term trends. R. globularis shows a general decreasing trend from middle to late Holocene (Figure 4d), whereas the opposite is observed for L. lobatula and C. refulgens (Figure 4d and e). In the BoB, R. globularis is commonly found in phytal environments in the vicinity of estuaries (Martínez-García et al., 2013; Pascual et al., 2008). R. globularis is a temporary attached species (e.g. to algae, marine grasses) secreting an organic glue to fix its test to the substrate, but it can free itself when searching for food or during sexual reproduction (Delaca and Lipps, 1972). Rosalina is one of the few genera that calcifies a special float chamber in a temporary planktonic phase, immediately before and during gametogenesis (Banner et al., 1985; Sliter, 1965). As such, this species can be easily transported by currents. This might explain the strong positive correlation between R. globularis and the D50 (i.e. strong hydrodynamics) especially before ~4.5–5 ka BP characterized by lower RSL stand (Figure 4a). However, with relative abundances reaching 30–40%, the transport alone is hardly conceivable. According to Sliter (1965) and Saraswat et al. (2011), sexual reproduction and thus also the floating chambers only occur at temperatures above 18°C. Therefore, the general cooling trend in summer SSTs recorded from middle to late Holocene in the BoB (from ~22 to <18°C; Zumaque et al., 2017), which is coherent with the decreasing trend in summer insolation (Figure 4a), might also explain the long-term decreasing trend of R. globularis (Figure 4d). In addition to temperature, the progressive deepening of the study site in response to RSL rise (~7 ± 1 m since 7 ka; Goslin et al., 2015; Figure 4a) might have moved shoreward the habitat of R. globularis. L. lobatula and C. refulgens habitats are not restricted to upper shelf environments but can extend to bathyal depths (García-Gallardo et al., 2017; Hayward et al., 2010). When attached on plants, L. lobatula and C. refulgens are permanently sessile (Langer, 1993; Langer and Gehring, 1993). By being able to inhabit greater water depths (Hayward et al., 2010), L. lobatula and C. refulgens might have adapted more easily to the deeper habitat. In the end, it is very likely that this long-term ecological benthic foraminiferal pattern results from the overlay of these environmental parameters.
On a shorter millennial time-scale, high density and diversity of benthic foraminiferal faunas are coherent with periods of (a) fine sedimentation (i.e. low D50 and log-Ti/Al) relatively rich in organic matter (i.e. high log-Br/Cl) and carbonate content (i.e. high Ca and Sr) and (b) high δ18O and δ13C values that we interpret in the following as more saline and potentially more productive waters (Figures 4 and 5). These features are coherent with our current knowledge of the ecological requirements of benthic foraminifera in the BoB. Indeed, fine-grained sediments, high labile organic content, low hydrodynamics, and normal open ocean salinity values (~35) have all been shown to influence positively the density and diversity of benthic foraminifera (e.g. Fontanier et al., 2002; Jorissen et al., 1995; Mojtahid et al., 2009; Saraswat et al., 2015).
Hydrological features implied from the foraminiferal isotopic signature
The core top values of δ13C and δ18O of the epifaunal C. refulgens (~1.6‰ and ~0.9‰, respectively; Figure 6) are in agreement with the modern isotopic signature of surface waters (0–100 m) in the BoB (Fontanier et al., 2006), implying that this species calcifies in equilibrium with bottom waters. The δ13C of Cibicides species reflects the δ13C of dissolved inorganic carbon (DIC) (e.g. Duplessy et al., 1984; Schmittner et al., 2017). The δ18O of foraminiferal calcite reflects the temperature of calcification and the δ18O of the seawater in which it calcifies, which in turn depends on salinity (e.g. Ravelo and Hillaire-Marcel, 1999). It is therefore impossible to deconvolve temperature and salinity effects without an additional independent temperature/salinity-proxy. On one hand, the significant positive correlation between depleted foraminiferal δ18O values and terrestrial source proxies (XRF-Ti, Fe) (Figure 5) may indicate a link to continental/freshwater input, and thus to a possible dominant salinity control over the δ18O signature during periods of high river runoff. That said, the δ18O–Fe/Ti link is partially biased by the difficulty to differentiate between the continental sources of Ti and Fe in our sediments (i.e. direct river discharges or upper shelf erosion during storm events). On the other hand, continental waters are depleted in both δ18O (water) and δ13C (DIC) in opposite to marine waters. Duplessy (1972) found that δ18O (water) and δ13C (DIC) in the Loire estuary are highly quasi-linearly correlated over the entire mixing range of marine and freshwater end members. Interestingly, the Loire freshwater has a very negative δ13C (~ –10‰) because most of the DIC stems from the remineralization of organic carbon in ground waters. Therefore, any variation in freshwater input to the core site, even though relatively distant from the Loire river mouth, should therefore leave an imprint in the δ13C of epibenthic foraminifera. In our record, δ18O and δ13C show a significant, but weak (r = 0.25; Figure 5) positive correlation, which might corroborate a partial control of continental freshwater input over the foraminiferal isotopic signature. That said, the weak δ18O and δ13C correlation implies that other factors are involved, especially for δ13C that is affected by local primary production (enriched in 13C). Therefore, periods of relatively high δ13C can indicate a more productive water column, perhaps triggered by fluvial nutrient inputs.

Comparison to intervals/episodes defined from other studies at the regional to global scale: (a) concentration of lithic grains from the work of Bond et al. (1997); (b) 3-point average curve of benthic foraminiferal δ18O at our study site. The pre and post ~2.5 cal. ka BP periods have separated scales in order to ignore the potential thermohaline front; (c) winter and summer sea surface salinities estimated from dinocysts in core MD95-2002 (Zumaque et al., 2017); (d) 3-point average density difference proxy from core Rapid-12-1K (Thornalley et al., 2009) calculated using derived Mg/Ca and δ18O temperatures and salinities of Globigerina bulloides and Globorotalia inflata; (e) sea surface temperature anomalies calculated using transfer function on planktic foraminiferal assemblages in core ODP Hole 658C (deMenocal et al., 2000). The interpretative horizontal line separating red and blue colored areas are placed around the average between maximum and minimum values of the total variability.
The major deviation recorded in the δ18O from ~3.5 to 2.5 cal. ka BP is not recorded in the δ13C (Figure 4f and g). This may imply an additional temperature control resulting probably from a major change in the hydrological configuration of the study area and probably of the BoB. At a regional scale, the 3.5–2.5 cal. ka BP time period coincides with the maximum of the Holocene transgression (e.g. Goslin et al., 2015; Ters et al., 1968) (Figure 4a). During periods of sea level high stand, freshwater coastal wetlands were inundated and coastal bays and rias were flooded, acting as sinks for fluviatile sediments (e.g. Morris and Mitchell, 2013; van Maren et al., 2016). Once the downstream valleys filled and the rivers channelized, suspended sediments and freshwaters were not anymore trapped by these coastal systems and therefore exported offshore leading to a significant reduction in coastal accommodation spaces (Clavé et al., 2001; Durand et al., 2016; Spencer et al., 1998; Traini et al., 2015). The resulting overall increase in freshwater load into the BoB certainly affected the hydrographical structure on the continental shelf leading eventually to the establishment of the modern winter thermohaline front located around 100 m water depth as previously suggested by Castaing et al. (1999) and Dubrulle et al. (2007). This front is the eastern boundary of warm and salty oceanic waters and the western boundary of colder and less salty shallow coastal waters (Castaing et al., 1999). According to Castaing et al. (1999), continental waters arriving at our site are materialized by major river plumes carrying organic-rich fine-grained sediments. In our core, this can be reflected by the finer grain-size distribution recorded after ~3.5–2.5 cal. ka BP (Figure 3a). Superimposed to this purely hydrological scenario, we cannot discard other possible factors such as climate forcing and early human impact triggering and/or amplifying the observed change. Indeed, the time period from ~3.5 to 2.5 cal. ka BP records a major storm event reported in northern Europe (Sorrel et al., 2012) and western Brittany (Van Vliet-Lanoë et al., 2014), and corresponds with increasing rainfall over the Loire and Gironde catchment areas (Naughton et al., 2007). Furthermore, several palynological studies have shown that human activity (e.g. forest clearing, livestock, crop farming) was already present at the time in the Loire catchment area (e.g. Visset et al., 2002), but it is only since ~2.0 cal. kyr BP that this activity impacted partly the terrigenous supply through soil erosion (Naughton et al., 2007).
After detrending the isotopic signal from the possible thermohaline front effect (Figure 6b), the δ18O variability of our oxygen isotopic signal seems also to be consistent with open oceanic North Atlantic records, at least for some specific periods (6.6, 5.9, 4.6, 4, 3.3, 2.5 cal. ka BP; Figure 6). This may further imply a probable link to larger scale hydrological patterns, such as those identified in the North Atlantic for the AMOC (Wanner et al., 2008). For example, at the Meriadzek Terrace (site MD95-2002 at ~300 km west from CBT-CS11; Figure 1), dinocyst assemblages and stable isotope data reveal alternating low and high salinity episodes since 11 ka BP (Figure 6c) (Zumaque et al., 2017). The low salinity episodes match dense sea-ice cover in the Labrador Sea (Solignac et al., 2004). Although the time-resolution is lower than in our study, we observe a general coherency between δ18O C. refulgens and the reconstructed sea surface salinities produced by Zumaque et al. (2017), except around 5.1 and after ~2 cal. ka BP (Figure 6). Interestingly, our low salinity events seem to match the cooling pulses defined by Bond et al. (1997) and characterized by episodes of increasing ice-rafted debris in the northern latitudes (Figure 6a). Furthermore, the overall good consistency between the heavier/lighter values of benthic δ18O signal (i.e. high/low salinity) at our site with periods of weak/strong SPG circulation (Thornalley et al., 2009) (Figure 6d) seems to corroborate the previous hypothesis. Indeed, strong SPG circulation results in a more east-west-oriented gyre, leading to more freshwater contribution (mostly from Greenland ice meltwaters) to the Atlantic inflow and a consequent weakening of the AMOC (Hátún et al., 2005). The opposite is observed when weak SPG results in a more north-south-oriented gyre (Staines-Urías et al., 2013; Thornalley et al., 2009). Periods of intensified AMOC lead to an active northward transport of warm and salty subtropical waters, which is also consistent with SST records from western subtropical coasts of Africa (deMenocal et al., 2000) (Figure 6f). Such connection between the BoB and the Labrador Sea has been suggested by Mary et al. (2017) at sites ProsecanKS10b and PP10-07 (Figure 1). Thus, we can also hypothesize that our isotopic signal might be also modulated by open water salinity changes (driven by AMOC variability).
Millennial-scale variability and controlling factors
Our data show a considerable high frequency variability superimposed on longer term trends. From middle to late Holocene, spectral analyses carried out on XRF-Ti, XRF-Br/Cl, XRF-Sr/Ti, δ18O, and δ13C signals revealed a statistically significant 1250-year frequency peak (Figure 7), but not a pervasive 1500-year cycle as reported in several studies from the North Atlantic (e.g. Bond et al., 1997; Debret et al., 2007; Thornalley et al., 2009; Wanner et al., 2008). The only exception is for Shannon index presenting a 1360–1760-year peak, which is probably because of the lower sampling resolution (Figure 7). Although Bond et al. (1997) concluded that pacings of the Holocene and glacial events are the same statistically, and together constitute a cyclic signal centered on 1470 ± 532 years, they specified that the mean pacing for Holocene events is 1374 ± 502 years. A similar periodicity of 1300 years was recorded in sediments of the Loire River mouth (core KV14-bis; Figure 1a), although the significance of the frequency peak is not very robust since this core covers only the past 2600 years (Durand et al., in press). Also, a periodicity of 1300 years has been reported by Pisias et al. (1973) and (Moros et al., 2006) in sediment cores from the North Atlantic. Debret et al. (2009) argued that frequencies with periods of 1250/1350 years are more common in the Southern ocean (e.g. Delmonte et al., 2002) whereas the 1500/1600-year cycles are more frequent in the North Atlantic. Debret et al. (2009) and Soon et al. (2014) conclude that the periodicities from 1250 to 1600 years most likely reflect an internal forcing which is distinct from the fundamental solar modes at 2300-year (Hallstattzeit), 1000-year (Eddy), and 500-year (unnamed) periodicities. Spectral analysis of the Holocene density stratification record of Thornalley et al. (2009) from the South Iceland rise (Figure 6d) shows one broad peak centered at 1500 years. These authors argued that the heat transport toward the northern latitudes has undergone large-amplitude millennial variability modulated by SPG dynamics. In addition to oceanic processes, a similar millennial cyclicity is also documented in continental records (Denton and Karlén, 1973; Solomina et al., 2015), and periods of intensified storminess in northern Europe (Sorrel et al., 2012), which strongly suggests an ocean–atmosphere connection. In our study site, NAO-like processes can represent a probable climatic forcing controlling both precipitations over the main French rivers’ catchment areas and therefore river runoff into the BoB, and wind stress that modulates SPG dynamics over the North Atlantic. Although the Holocene millennial climate variability is difficult to reconcile with the meteorological NAO dynamics that favors annual to decadal time-scale variations (Hurrell et al., 2003), North Atlantic Holocene records and climate model simulations show that millennial-scale hydrographic changes exhibit spatial patterns that reproduce fairly well the seesaw associated with the NAO (Jessen et al., 2011; Rimbu et al., 2003; Schmidt et al., 2004; Schulz et al., 2007; Thornalley et al., 2009). In addition to the NAO, the AMO (Mann et al., 2009) has been identified as a coherent pattern of oceanic oscillatory changes in North Atlantic SSTs at multidecadal time-scales, displaying a spatial pattern resembling that of the atmospheric NAO mechanism (Mann et al., 2009; Sutton and Hodson, 2005). This link has been hypothesized in previous studies from the BoB (e.g. Durand et al., in press; Garcia et al., 2013; Mary et al., 2017; Mojtahid et al., 2013; Zumaque et al., 2017). However, every region and every proxy record has a different sensitivity to these mechanisms that are further modulated by local changes in the sedimentation (e.g. sedimentation rate, diagenesis, etc.).

Continuous Wavelet Transform (CWT) analyses (periodicities) and REDFIT analyses (power spectra) performed on (a) log-Ti, (b) log-Br/Cl, (c) log-Sr/Ti, (d) Shannon H, (e, e1, e2) δ18O, and (f) δ13C. For CWT, the cone of influence is plotted to show the region where boundary effects are present and the signal power (squared correlation strength with the scaled mother wavelet) is shown in color scale.
There appear to be other significant higher frequency oscillations corresponding to multi-century periods of which the 450/500 years is a recurrent peak (Figure 7). This periodicity was also recorded in core KV14-bis (Figure 1a) in front of the Loire estuary and a strong link to solar forcing was hypothesized (Durand et al., in press). In fact, a ~500-year Holocene cycle in solar activity was inferred from variations of global atmospheric 14CO2 production (Stuiver et al., 1997). It has been strongly suggested that the ~500-year periodic solar irradiance plays a key role in cyclic oceanic and atmospheric change including NAO-like processes.
Conclusion
We present a paleoenvironmental history for the northern BoB for the past 7 kyr. The benthic foraminiferal fauna is dominated by three epiphytic/epibenthic species (R. globularis, L. lobatula, and C. refulgens). Through the record, absolute and relative densities of the three species display a long-term mid- to late-Holocene pattern with R. globularis showing a general decrease whereas the opposite is observed for L. lobatula and C. refulgens. The pattern of R. globularis can be explained by the overall SST cooling and/or the progressive deepening of the water column due to RSL rise. By being able to inhabit greater water depths, L. lobatula and C. refulgens might have adapted to deeper habitats with time. In addition to these long-term benthic foraminiferal patterns, a millennial-scale variability is recorded by all our proxies. This variability seems to be closely linked with changes in precipitations over the catchment areas of the main French rivers and/or with AMOC variability, itself modulated by SPG dynamics. We further suggest a NAO-like process controlling both humidity over the main rivers’ catchments areas and SPG dynamics through wind stress. Solar forcing seems to be superimposed to this millennial variability.
Foraminiferal δ18O and grain-size analyses show a major shift toward heavier isotopic values and finer sedimentation around 3.5–2.5 cal. ka BP. We hypothesize a link with the near-stabilization of sea-level rise since ~3.5–2.5 cal. ka BP, resulting in the reduction of estuarine accommodation spaces and the channelization of the main rivers outflow. This likely caused extended freshwater and suspended sediment input into coastal waters of the BoB, triggering therefore the formation of the present-day winter 100-m isobath thermohaline front. This time period also coincides with increased storminess in northern Europe and Brittany, and increased rainfall on the Loire River catchment area. These factors may have triggered and/or amplified the observed change by affecting the hydrology of the BoB.
Supplemental Material
Supplemental material for Millennial-scale Holocene hydrological changes in the northeast Atlantic: New insights from ‘La Grande Vasière’ mid-shelf mud belt
Supplemental material, Supplemental_Material for Millennial-scale Holocene hydrological changes in the northeast Atlantic: New insights from ‘La Grande Vasière’ mid-shelf mud belt by Meryem Mojtahid, Matthieu Durand, Pierre-Olivier Coste, Samuel Toucanne, Hélène Howa, Jean Nizou, Frédérique Eynaud and Aurélie Penaud in The Holocene
Footnotes
Funding
This study was supported by the French CNRS INSU and is a contribution to the 2013–2014 LEFE-IMAGO project ‘HCOG2’ (‘Forçages climatiques Holocène et répercussions Côtières et Océaniques dans le Golfe de Gascogne’; coord. A. Penaud) and by the Nantes Saint-Nazaire seaport (project ‘DETAiLH’: DETermination de l’état pré-Anthropique de l’estuaire de la Loire et évolution Historique de son fonctionnement, coord. M. Mojtahid). 14C dates were obtained, thanks to French national ARTEMIS 14C AMS facilities. We would like to thank Oanez Lebeau (IUEM, Brest) from the ‘Pôle Spectrométrie Océan’ (IUEM, CNRS, IFREMER) for her help with isotope analyses. We also thank IFREMER for core acquisition and laboratory facilities. This study represents a contribution to the French ANR project HAMOC. This work was also supported by the ‘Laboratoire d’Excellence’ LabexMER (ANR-10-LABX-19) and co-funded by a grant from the French government under the program ‘Investissements d’Avenir’. Raw data are available as
or by contacting the first author (e-mail:
Supplemental Material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
