Multi-proxy reconstruction of climate changes in the Lower St. Lawrence Estuary, Canada, during the Middle and Late-Holocene

The micropaleontological and palynological content, and geochemical and isotopic composition of a marine sediment core collected off Pointe-des-Monts in eastern Québec, Canada, reveal regional palaeoclimatic and paleoceanographic conditions in the Lower St. Lawrence Estuary over the last ~8200 years. The pollen and spore content allows comparison with the terrestrial palynostratigraphy, whereas dinoflagellate cysts and benthic foraminifera are used to reconstruct sea-surface conditions and bottom water properties, respectively. The dinocyst-based reconstructions indicate shifts between estuarine and oceanic conditions with important changes in sea-surface temperature, salinity, and primary productivity. Both the dinocyst assemblages and the quantitative sea-surface estimates highlight a distinct transition at ca. 4200 cal years BP. It is notably marked by a change towards higher salinity, which suggests reduced freshwater discharge, hence lower precipitation in the watershed, during the Late-Holocene. The isotopic composition (δ18O and δ13C) and assemblages of the benthic foraminifera indicate centennial to millennial frequency variability of bottom water properties, over a general trend towards decreasing temperatures and increasing ventilation from the beginning of the Middle Holocene until the last century. Since then, reverse trends with abrupt warming and decreasing dissolved oxygen content in bottom water have been observed.


Introduction
Palaeoecological studies from the Gulf of St. Lawrence area, Quebec, Canada, as well as paleoceanographical investigations in the Estuary and the Gulf of St. Lawrence have shown that important fluctuations in terrestrial climate (e.g.Fréchette et al., 2021;Magnan and Garneau, 2014;Perrier et al., 2022;Sauvé, 2016), vegetation (e.g.Asnong and Richard, 2005;Fréchette et al., 2021;Lavoie and Filion, 2001;Marcoux and Richard, 1995;Richard, 1994;Sauvé, 2016), sea-surface conditions (e.g. de Vernal et al., 1993;Dhahri, 2010;Levac, 2003;Wu et al., 2022) and bottom water properties (e.g.Dhahri, 2010;Genovesi et al., 2011;Thibodeau et al., 2013) have characterized the Holocene.In the nearshore context of the Gulf of St. Lawrence, close relationships between surface ocean temperatures and climate conditions can be hypothesized.On one hand, the heat budget is modulated by insolation and exchanges between the surface ocean and the atmosphere.In particular, the winter sea-ice cover regulates the energy budget through its albedo and influences the exchanges of energy and water vapour at the ocean-atmosphere interface (Faucher et al., 2004;Wu et al., 1997).On the other hand, precipitation within the watershed determines the freshwater budget, salinity, and water mass stratification (Bourgault and Koutitonsky, 1999), while evaporation at the sea surface is a source of moisture to the atmosphere.
The St. Lawrence Estuary (SLE) and Gulf of St. Lawrence (GSL) furthermore constitute a transitional marine environment between the St. Lawrence River and the North Atlantic Ocean.Strong stratification of the water column characterizes these environments, where surface waters are characterized by low salinity due to freshwater inputs from the St. Lawrence River and its tributaries, while the intermediate and deep layers result from the mixing of variable proportions of waters from the Gulf Stream and the Labrador Current.Hence, both the SLE and GSL are marked by a strong stratification fostering low thermal inertia and a very pronounced seasonal cycle in the surface layer, which records freezing conditions in winter and warming up to 18°C in summer (Saucier and Chassé, 2000).
In this article, we aim to document the relationship between hydrographic conditions in the SLE and GSL versus the regional climate during the present interglacial.In order to reconstruct simultaneously climate, vegetation, and oceanic changes, we have analysed the palynological content of a core collected off Pointe-des-Monts at the mouth of the Lower St. Lawrence Estuary (LSLE) in Quebec (Figure 1).Special attention was paid to the cysts of dinoflagellates (hereafter referred to as dinocysts), which are planktonic organisms with a distribution closely related to sea-surface conditions such as temperature, salinity, and sea-ice cover (de Vernal and Marret, 2007;de Vernal et al., 2020).We also analysed the pollen and spores originating from the onshore vegetation to allow comparison with the terrestrial palynostratigraphy (de Vernal et al., 1993).We furthermore examined the benthic foraminiferal assemblages and the isotopic composition (δ 13 C and δ 18 O) of a species (Bulimina marginata) present throughout a large part of the record.These tracers allow us to reconstruct changes in the properties of bottom waters which are governed by the mixing of the North Atlantic Water with the Labrador Current and thus record changes in the ocean circulation along the continental margins of eastern North America (e.g.Bilodeau et al., 1994;Boyle and Keigwin, 1987;Genovesi et al., 2011).

Physiography and hydrography of the study area
The coring site is located in the Laurentian Channel at the seaward limit of the LSLE, 22.26 km southeast of Pointe-des-Monts (Figure 1).
The SLE constitutes a transitional environment between the St. Lawrence River and the northwest North Atlantic Ocean (Dickie and Trites, 1983;Koutitonsky and Bugden, 1991).The dominant bathymetric feature of the SLE and GSL is the Laurentian Channel, a deep submarine valley >250 m deep, that extends over 1280 km from the continental shelf off Nova Scotia and Newfoundland in the northwest North Atlantic to Tadoussac, upstream in the SLE.Seaward, the GSL forms a semi-enclosed sea of ~226,000 km 2 .Southeastward, it is connected to the North Atlantic Ocean through the ~104 km-wide Cabot Strait, where the Laurentian Channel bathymetry reaches ~408 m.Northeastward, it is connected to the Labrador Sea, through the ~16 km-wide Belle Isle Strait, with a ~60 m-deep sill.
Three water layers characterize the LSLE and GSL (Dickie and Trites, 1983;Koutitonsky and Bugden, 1991).The surface water layer features low salinity (27.0-31.5 psu) due to freshwater from the St. Lawrence River and tributaries surrounding the Estuary and Gulf.This layer is less than 50 m deep.It has low thermal inertia and thus records large amplitude changes of temperature throughout the year, from partial freezing in winter, up to 18°C in August.Underneath, the intermediate layer extends down to about 150 m.This layer develops from winter cooling and isopycnal penetration of cold waters from the Labrador Sea, through the Cabot and Belle-Ile straits.It has a mean temperature and salinity of ~0°C and ~32 psu respectively (Gilbert and Pettigrew, 1997).By the end of December to early January, the surface and intermediate layers partly merge, due to the cooling and vertical mixing of surface waters (Gilbert and Pettigrew, 1997).The warmer (4°C -6°C) and saltier (35 psu) bottom water layer flows from Cabot Strait landward in the Laurentian Channel, up to Tadoussac, at the eastern limit of the upper SLE (Figure 1).This bottom water results from the mixing of the Labrador Current Water (LCW) and North Atlantic Current Waters (NACW) (Gilbert et al., 2005).At the study site (Figure 1), summer sea-surface temperature and salinity average 10.30 ± 2.89°C and 27.69 ± 1.40 psu, respectively, according to the World Ocean Atlas 2013 (Locarnini et al., 2013;Zweng et al., 2013).The seasonal sea-ice cover lasts for about 3 ± 0.5 month per year (Galbraith et al., 2019), and the annual productivity averages ~925 g C m −2 , based on MODIS data spanning the 2002-2017-time window (http://daac.gsfc.nasa.gov).

Sampling and laboratory treatment
The study core COR0602-42PC (49°7.18′N -67°16.80′W; water depth 332 m) hereafter referred to as P42 has been collected in the Laurentian Channel off Pointe-des-Monts during the 2006 cruise of the Coriolis II.This core is 7.28 m long and was retrieved using a modified version of a Benthos™ piston corer of 10 cm diameter and 9 m long.The core was stored in 1-meter-long plastic liner sections and placed in a cold room at 4°C at the Institut des sciences de la mer de Rimouski (ISMER).
After splitting, the sediment core was subsampled at 1 cm intervals and analysed at 8 cm intervals for its palynological content.5 cm 3 of wet sediment subsamples were weighed, then airdried and weighed again.The dry weight was used for calculating concentrations of palynomorphs and foraminifera, as well as the dry bulk density of sediment subsamples.The subsamples were then sieved on 106 μm and 10 μm mesh sieves.Palynological preparation followed the procedures of de Vernal et al. (1996).The 10-106 µm fraction was treated with hydrochloric acid (HCl 10%) and hydrofluoric acid (HF 49%) at room temperature to dissolve carbonates and silica particles.One tablet with a known number of Lycopodium spores (batches #483216 and #177745, see Supplemental Table T1) was added to each sample before chemical treatments for concentration calculations (Stockmarr, 1971).The spiked spores are distinguishable from in-situ Lycopodium clavatum with their darker colour.A final sieving was done through a 10 μm mesh sieve and a drop of glycerine jelly was added to the residue, which was mounted between the slide and cover slide for further analyses under microscope.The fraction >106 µm was dried out and kept for foraminifer studies.

Chronology of the core
The age model for the core was initially established by Barletta et al. (2010) from five AMS 14 C measurements of shell fragments.For the present study, one shell fragment and three samples of mixed benthic foraminifera were also analysed at the Keck-CCAMS facility (Irvine, USA) and the A.L. Lalonde AMS Laboratory (University of Ottawa, Canada).The rarity of biogenic carbonate prevented more 14 C measurements from being made for a more precise chronology of the core.Thus, in total, nine AMS 14 C dates (Table 2) were available.They were calibrated with the Marine20 calibration curve (Heaton et al., 2020) and a regional offset (ΔR) of 25 ± 72 years, calculated from the weighted mean of the four reference points (McNeely et al., 2006) from the LSLE reported in the database of Reimer and Reimer (2001).We used the Bacon package (Blaauw and Christen, 2011) available in R software to build an age model.

Palynological and micropaleontological analyses
A total of 91 samples were analysed.Palynomorph identification and counts were performed using a transmitted light microscope at 400× magnification.All palynomorphs were systematically counted.They mostly include pollen grains, spores, dinocysts, and organic linings of foraminifera.They also include Halodinium, which is likely related to ciliates (Gurdebeke et al., 2018).Here, we present the results from pollen and spore analyses in addition to dinocyst counts.At least 300 pollen grains and spores per sample were identified at the genus or family level using the nomenclature of Richard (1970), McAndrews et al. (1973), and Moore et al. (1991).About 300 dinocysts per sample were identified at the species or genus level using the nomenclature of Rochon et al. (1999)  The identification of benthic foraminifera used the taxonomy of Loeblich and Tappan (1988) and Rodrigues et al. (1993).We illustrate in a diagram the percentages of taxa that form ⩾3% of the assemblages in at least one sample per core.Some samples contain a low number of benthic foraminifer tests, especially at the bottom of the core.The assemblages containing less than 20 specimens were not included.Detailed results of systematic counts for all taxa can be found in Supplemental Table T1.
Multivariate analysis was performed on dinocyst, pollen, and benthic foraminifer assemblages using the CANOCO 5 software (Ter Braak and Šmilauer, 2012) for defining the zonation.A square root transformation was applied to pollen and benthic foraminifer percentages and a logarithmic (base 10) transformation was applied to dinocyst percentages.The transformation of our raw data allowed us to down-weigh the influence of overwhelmingly dominant data in the analysis.Data analysis in CANOCO 5 suggested that the linear method and principal component analysis (PCA) were appropriate for assessing affinities within the dinocyst, pollen, and benthic foraminifera data sets with gradient lengths of 0.9, 0.8, and 2.7 SD units respectively.The PCA analysis was performed on benthic foraminiferal taxa recording ⩾3% of the assemblage in at least one sample.See Supplemental Table T2 for assemblage data in percentages and results of PCA.

Reconstruction of sea-surface conditions with the modern analogue technique.
In order to reconstruct sea-surface conditions from dinocyst assemblages, we used the Modern Analogue Technique (MAT) as described by Guiot and de Vernal (2007).MAT is based on the similarity between fossil and modern spectra with the assumption that fossil populations evolved in comparable environmental conditions to their modern analogues (Guiot and de Vernal, 2007).The reference modern dinocyst database used here is a set from 1968 sites mostly located in the Northern Hemisphere (de Vernal et al., 2020).It includes 71 taxa and covers a wide range of environmental conditions with respect to i) summer and winter sea-surface temperatures (SST), ii) sea-surface salinity (SSS), iii) sea-ice cover duration as expressed in months per year, and iv) primary productivity (de Vernal et al., 2005;Guiot and de Vernal, 2007).We used MAT implemented in the R software (http://www.r-project.org/)following the procedure described in Guiot et al. (2013).
Prior to the search for analogues, the dinocyst percentage data were transformed using a natural logarithm to increase the weight of accompanying taxa and reduce the weight of cosmopolitan taxa.The five best analogues were selected from the reference database to calculate past sea-surface conditions from their average, weighted inversely to their distance.The distance of the five best analogues was used as an index of the reliability of the reconstruction.The threshold distance for significant non-random analogue with the n = 1968 database is 1.2.Validation exercises indicate root-mean-square errors of prediction of ±1.5 month year −1 for the sea-ice cover, ±1.9°C for summer SST, ±1.9 psu for SSS, and 125 g C m −2 year −1 for annual primary productivity.Modern SST and SSS data were compiled from the World Ocean Atlas 2013.The sea-ice cover was compiled from the 1955 to 2012 data provided by the National Snow and Ice Data Center.It is expressed in months per year with >50% of sea-ice concentration considered as sea-ice-present.Primary productivity data result from the vertical generalized production model (VGPM) algorithm applied to the 2002-2017 Moderate Resolution Imaging Spectroradiometer (MODIS) chlorophyll data (http://daac.gsfc.nasa.gov)and are expressed in g C m −2 year −1 .The results of reconstructions can be found in Supplemental Table T5.

Geochemical analyses
Oxygen (δ 18 O) and carbon (δ 13 C) isotopes analyses on benthic foraminiferal tests.The identification of calcareous benthic fora- minifera was done using a stereoscopic binocular microscope at 40x magnification.In each sample, a minimum of eight tests of Bulimina marginata were hand-picked for stable isotope ( 18 O and 13 C) analyses at the Geotop-UQAM Research Center.The tests were first combusted under vacuum at 250°C for approximately 1 h.The analyses were then performed on a Micromass Isoprime™ isotope ratio mass spectrometer in dual inlet mode, coupled to a MultiCarb™ preparation system.CO 2 was extracted at 90°C by acidification with 100% phosphoric acid (H 3 PO 4 ).Measurements, made against an internal reference carbonate material calibrated on the VPDB scale and calculation of δ‰ values followed Hillaire-Marcel et al. ( 2021).The analytical reproducibility was routinely better than 0.05‰ (one sigma) for both carbon and oxygen isotopes.Results can be found in Supplemental Table T4.
Sediment carbon and nitrogen analyses (%C org , C org /N, δ l3 C).Between 7 and 10 mg of dried and crushed sediments from each sample were analysed following the HCl-fumigation method, at Geotop-UQAM (Hélie, 2009).This method consists of placing the sediment in tin cups and analysing it with a Carlo Erba™ NC 2500 elemental analyser to determine the weight percent of total carbon (C tot ) and total nitrogen (N tot ) contents.The organic carbon content (C org ) was measured from sediments acidified after fumigation.Subsamples were weighed into silver cups and placed in a clean Teflon tray.The tray was put in a closed glass container with 4 N HCl beaker for 24 h.During this interval, the HCI that saturated the air in the container reacted with the carbonate minerals.After fumigation, the silver cups were closed, wrapped in tin cups, and then analysed with a Carlo Erba™ NC Table 2. Radiocarbon dates and calibrated ages with Marine20 calibration curve (Heaton et al., 2020) 2500 elemental analyser.The results are given in percent (%) of the total weight.Replicate measurements of organic analytical reference materials determined the analytical precision, which is estimated at ±0.1% for C org and ±0.3% for N. The carbon isotope composition of C org (δ 13 C org ) was measured with a Carlo Erba™ elemental analyser in-line with a Micromass Isoprime™ ratio mass spectrometer.Data are reported in the δ ‰ notation versus VPDB, as above, but with an analytical uncertainty of ±0.1‰.

Chronology and lithology
The sediment recovered in core P42 contains mainly fine hemipelagic mud.As described by Barletta et al. (2010), the mean grain size ranges between 2.25 and 6 µm.Both the grain size and magnetic susceptibility parameters show relatively uniform distribution through time, which suggests little changes in the sedimentation regime and thus a relatively constant sedimentation rate for the entire record.The age model based on AMS 14 C dates indicates an apparent steady average sedimentation rate of ~0.875 mm year −1 for the last ~8200 years (Figure 2).However, the dry bulk density of the sediment increases with depth (see Figure 2, Supplemental Table T1), which suggests compaction at depth and hence a slight upward decline in sedimentation rate.All concentrations are reported versus dry weight to avoid any bias due to compaction.A critical issue in most marine records is to estimate the mean 14 C age of the core top mixed layer and the depth of this layer, which are needed to assess the smoothing function linked to mixing by bioturbation.This function depends on the time range considered (e.g.Boudreau, 1986).As documented in Supplemental Figure F1, 228 Th/ 232 Th and 210 Pb data from a box core (49°05.24N, 67 °26.22W; water depth 311 m) collected at closeproximity (12 km) to core P42, show 228 Th-excess versus 228 Ra down to ~4.5 cm, 228 Ra-deficit versus 232 Th down to ~8 cm and, 210 Pb in excess versus its supported fraction, down to ~17 cm (Jennane, 1992).The time ranges defined by these disequilibria (~5 half-lives of the shortest-lived isotope) are ~10, ~30, and ~100 years, respectively.When confronted with the mean sedimentation rate estimated from 14 C data (i.e.0.875 mm year −1 ), potential mixing down to ~17 cm, possibly more at 14 C time scale due to non-local mixing (Boudreau, 1986), would produce a smoothing window of about 200 years, possibly more.Assuming relatively uniform mixing processes and good recovery of the sediment including the surface, all age estimates indicated below should be considered with a broad uncertainty of about ±100 years.

Biogenic content of sediment
Pollen concentrations range from 7200 to 140,000 grains g −1 with accumulation rates averaging about 2700 grains cm −2 year −1 (Figure 3).Throughout the core, the dinocyst concentrations vary from 685 cysts g −1 at the base (8180 cal years BP) to 12,672 cysts g −1 (core top) with an average of ~4000 cysts g −1 .From 8180 cal years BP to 3100 cal years BP, dinocyst concentrations are low (3134 cysts g −1 ) but recorded an increase between 5200 and 4600 cal years BP (5000 cysts g −1 ).From 3100 cal years BP to the  Blaauw and Christen (2011) and the 'Marine20' calibration curve (Heaton et al., 2020), with a regional offset (ΔR) of 25 ± 72 years (see text).Semi-transparent blue zones represent calibrated age ranges as in Table 1.The red dashed line represents the median values of the age model, and the grey dashed lines represent 95% confidence ranges.Right: dry bulk density (g cm −3 ) versus depth, calculated as the ratio of the mass of dried sediment to its total wet volume.top, a steady increase in cyst concentrations is observed (1600 cysts g −1 to 12,600 cysts g −1 , averaging 6000 cysts g −1 ).The concentration of organic linings of benthic foraminifera ranges from 0 to 9000 linings g −1 with values >5000 linings g −1 around 4940 and 950 cal years BP.Concentrations of benthic foraminiferal tests are however low (<75 tests g −1 ), especially from 8180 to 6700 cal years BP (<5 tests g −1 ).The dominance of pollen and spores over dinocysts in the palynological assemblage of core P42 is characteristic of estuarine environments (de Vernal and Giroux, 1991; see dinocysts/pollen graph in Figure 3).
The C org content of the sediment shows a slow and continuous increase from 0.27% at the base of the core to 1.17% near the surface.Such an increase in C org may be associated either with a recent increase in organic matter fluxes or by the diagenetic effect increasing with depth in the sedimentary column (Muzuka and Hillaire-Marcel, 1999).The absence of major variations in organic carbon content suggests regular sedimentary inputs in agreement with the age model suggesting uniform sedimentation rates (Figure 2).The δ 13 C org averages -22.90‰ and corresponds to values typical of marine organic matter (OM) with a low contribution of terrestrial OM (Muzuka and Hillaire-Marcel, 1999).The δ 13 C org curve shows two notable changes.A decrease from -21.9 to -23.2‰ occurred between ~8180 and ~7580 cal years BP.It would represent an increase in terrestrial organic matter input as suggested by a gradual increase in pollen concentrations.From 1370 cal years BP to core top, δ 13 C org values increase from −23.1 up to −22.2‰.The organic carbon to total nitrogen ratio (C:N) varies from 4.16 to 9.08, averaging 8.37.There is a slightly increasing trend from core bottom to core top with two depletion episodes, one between 8180 and 7580 cal years BP, and the second from ~560 cal years BP to core top (~200 cal years BP).Both the δ 13 C org and the C:N data suggest a higher proportion of marine over terrestrial organic matter, except at the bottom and top of the core.

Pollen and spore assemblages
The pollen assemblages are largely dominated by coniferous tree taxa (Picea, Pinus, Abies; Figure 4a).Miscellaneous deciduous tree taxa include mainly Acer, with low occurrences of Carpinus, Juglans, Tilia, and Fagus.Alnus includes Alnus crispa and Alnus rugosa.Among spores, Lycopodium spp.includes Lycopodium annotinum and Lycopodium clavatum.The trilete spores include those other than Lycopodium.
Pollen assemblage zones (PI, PII and PIII) were defined based on downcore variations of sample scores on the first and second principal components (PC1 and PC2), which explain respectively 48% and 25% of the variance (Figure 4b).PC1 is characterized by a distinct opposition between Pinus, which has a very low negative loading, and Picea with a very high positive loading.PC2 shows an opposition between conifer taxa on the negative side and all deciduous taxa on the positive side.Among the deciduous taxa, Betula has the highest loading on PC2.
Zone PIII (<3500 cal years BP) records an increase of PC1 sample scores with an increase of Picea (from 31% to 56%) compensated by a decrease of Pinus (from 49% to 24%).In this zone, Picea gradually dominates over Pinus.PC2 sample scores are generally negative but drastically increase over the last 1150 years, which corresponds to a steep decline in Abies (7 to 0%), followed by a maximum of Alnus (7%) and Salix (7%).Finally, Ambrosia-type occurred in the uppermost part of the core, after ca.270 cal years BP.
The overall pollen assemblages of core P42 are characterized by a gradual increase of Picea and Abies percentages and a decrease of Pinus over the last 8180 years.Such a trend is also documented from pollen records of adjacent terrestrial areas (e.g.Asnong and Richard, 2005;Lavoie and Filion, 2001;Marcoux and Richard, 1995;Sauvé, 2016).
The composition of assemblages is consistent with that reported from other marine records of the GSL and SLE (e.g. de Vernal et al., 1993;Dhahri, 2010;Thibodeau et al., 2013).The record is dominated by Nematosphaeropsis labyrinthus, Operculodinium centrocarpum and the cyst of Pentapharsodinium dalei in addition to Islandinium minutum and Brigantedinium spp.Such an association bears a postglacial signature, younger than ca.8500 cal years BP at the regional scale, as described by de Vernal et al. (1993) based on the occurrence of N. labyrinthus.
In core P42, two main assemblage zones with a well-defined transition can be distinguished based on downcore variations of PC1 and PC2, which explain respectively 49% and 13% of the variance (Figure 5b).PC1 is marked by a strong opposition between the cyst of Protoperidinium americanum, Spiniferites elongatus, Brigantedinium spp., Selenopemphix quanta on the positive side, and I. minutum, P. dalei, Spiniferites ramosus on the negative side.On the PC2 axis, the cysts of P. americanum and P. dalei and I. minutum have positive loadings, while all other taxa have negative loadings.
Zone DI spans from the base of the core to ca. 4200 cal years BP.It is characterized by generally positive scores of PC1 that gradually decrease.It corresponds to decreasing percentages of S. quanta (from 25% to 4%) and Brigantedinium spp.(from 30% to 15%) while I. minutum percentages increase (2%-31%).
A transition at ca. 4200 cal years BP is well marked by a steep decrease in both PC1 and PC2 sample scores, which is related particularly to the disappearance of the cyst of P. americanum, together with the decrease in S. quanta.In Zone DII, the negative PC1 scores correspond to high percentages of I. minutum (14%-48%), which dominates the heterotrophic taxa, and low percentages of S. elongatus (1%-6%), Brigantedinium spp.(4%-19%) and S. quanta (0%-6%).A small recurrence peak of the cyst of Protoperidinium americanum (up to 11%) is centred around 950 cal years BP, thus contributing to high PC2 scores.

Reconstruction of sea-surface conditions
Sea-surface conditions reconstructed from MAT (Figure 6) appear reliable as they are based on high dinocyst counts (>300) and because all samples have close modern analogues (distance <0.8).The analogues are mainly from the GSL and ESL.The reconstructed summer sea-surface temperature (SST) does not depart much from the modern values (10.30 ± 2.89°C).They vary mostly from 9 to 12°C and average 10.9°C throughout the record.However, a distinctive cold interval (8.7-9.9°C) is recorded from ~8180 to ~7290 cal years BP.It is about 1°C cooler than the modern average value.From ~7290 to about 4600 cal years BP, the summer SST progressively increased from 9.6 to 11.8°C.This warming trend was interrupted by a cool event of 1.5°C centred at ~4200 cal years BP.In the upper section, there is a quasi-cyclical pattern of summer SST with maxima (>11.5°C)recorded around 3750-3000 cal years BP and 2000-1050 cal years BP, and then a cooler interval from ~950 to ~490 cal years BP.Summer SSS shows large variations following a pattern similar to that of summer SST, except from ca. 4090 to 1480 cal years BP.
The analogues suggest a mean SSS of 28.6 psu, with values as low as 25.7 psu at ~8100 cal years BP to values around 30 psu during the ~4090-1480 cal years BP interval.Reconstructed seaice cover averages 3.3 months year −1 .It varies generally in parallel with summer SST, indicating freezing conditions of up to 5.3 months year −1 at summer SST minimum around 2130 cal years BP.Reconstructed productivity ranges from 331 to 1037 g C m −2 year −1 , with an average of 674 g C m −2 year −1 , and appears inversely related to SSS.The transition at ca. 4200 cal years BP in the assemblages leads to the reconstruction of rapidly increasing SSS and decreasing productivity, relatively low SST, and long sea-ice season.
Zone FI (8180-6440 cal years BP) is characterized by highly variable but generally positive PC1 scores, while PC2 scores fluctuate around 0. The assemblages contain many rare minor taxa, and the concentrations are too low to reach statistically significant counts as most samples contained no more than 20 specimens in total and are not illustrated in Figure 7. Very low concentrations of foraminifer tests (1-22 tests g −1 ) may indicate low benthic productivity, but the high recovery of organic linings reaching up to 2710 g −1 (Figure 3) may indicate calcium carbonate dissolution in the sediment.Zone FI is dominated by B. marginata.
The transition from zone FII to zone FIII (3630-200 cal years BP) is marked by an abrupt increase of PC2 sample scores, which remain high thenceforth.This corresponds to a rapid increase of G. auriculata, I. nocrossi and N. labradorica in compensation to a sharp decline of E. exilis and quasi-disappearance of O. umbonatus.Also of note, is the decrease around 430-200 cal years BP of B. marginata, N. labradorica, and E. clavatum relative to the increase of low-oxygen indicators B. subaenariensis (0 to 40%) and E. exilis (4 to 15%).

Isotopic composition of Bulimina marginata tests
Isotopic analyses were performed on tests of B. marginata, which is the most common taxon of the entire record, but many samples contained too few tests (<8) to allow reliable stable isotope analyses.Therefore, the record is discontinuous especially below 536 cm (>6350 cal years BP).Nevertheless, some trends can be observed.
The δ 13 C values of B. marginata vary slightly, ranging from −0.05 to 0.39‰ with an average of 0.17‰ (Figure 8).They show sub-millennial variability within an overall range well-below isotopic equilibrium with the total dissolved organic carbon of the North Atlantic (about +1‰ today, cf.Becker et al., 2016).This offset points to mixing with some 13 C-depleted CO 2 released to bottom water and in the pore waters of sediment by the oxidation of the organic matter (see de Vernal et al., 1991;Lucotte et al., 1991).A positive shift of δ 13 C org during the last millennium indicates enhanced fluxes of marine organic carbon to the sea floor, a signal which seems to be captured by the δ 13 C of B. marginata.The δ 18 O of B. marginata ranges from less than 3‰ in the lower part of the core to well above 3‰ in the upper part of the core, except for the last two centuries.A linear regression analysis confirms the statistically significant increasing trend with a p-value of 0.0012 (see Figure 8).The driver of this trend could be a progressive change in the isotopic composition and/or of the temperature of bottom waters filling the deep Laurentian Channel.Whereas a general cooling of the northwest North Atlantic waters has been observed in the Middle and Late-Holocene (Sachs, 2007), a stronger contribution of the cooler Labrador Current (see Figure 8; Lochte et al., 2019Lochte et al., , 2020) ) to the deep-water mass filling the Laurentian Channel could be invoked to explain the 18 O-encrichment throughout the recorded interval.This is also supported by the relative increase of N. labradorica in the benthic foraminiferal assemblages during the same time period.
The comparison of δ 18 O in benthic foraminifera from our study site with that of core MD99-2220 in the SLE (Figure 8) shows converging values during the last ~1000 years and overlap during short intervals at ca. 4000 and 2000 cal years BP.Otherwise, a ~0.2‰ offset persisted throughout most of the Middle to Late-Holocene between the two sites with slightly 18 O-enriched compositions at the site of core MD99-2220, situated about 112 km upstream from core P42.Slightly cooler conditions could thus have prevailed at the MD99-2220 site.Nonetheless, a strong reverse trend is observed near the core top at the two sites.Its precise duration is unclear using the Holocene age models of the two cores.However, based on Thibodeau et al. (2010), this reverse trend should date from the early 19th century, that is, close to the ~150 cal years BP duration observed in Figure 8.

Comparison of the proxies
The PCA analyses highlight with their first(s) component(s) the most important changes in assemblages, which we assume to be related to environmental factors.In Figure 9, we compare the downcore variation of the first component (PC1) for the three micropaleontological tracers of climate and hydrographical conditions used in our study, which include pollen, dinocyst, and benthic foraminifer that respectively relate to land vegetation, surface waters, and bottom waters.The three tracers have recorded similar trends of their PC1, especially at the beginning of the Middle Holocene until about 6400 cal years BP, suggesting that terrestrial climate, sea-surface and bottom water conditions in the LSLE area all shifted in response to the same large-scale forcing.Given the regional context, we may hypothesize that it corresponds to atmosphere-ocean reorganization after the final ice retreat over eastern Canada (e.g.Dyke, 2004).After about 6400 cal years BP, the PC1 of foraminifer assemblages remained relatively stable while the PC1 of pollen and dinocyst assemblages continued to shift, suggesting a possible decoupling between bottom water and land/sea surface conditions.After the Middle-Late-Holocene transition about 4200 years ago, the PC1 of pollen and dinocyst assemblages displayed distinct variations, pointing to different mechanisms and response times, as discussed in the following sections.

Postglacial vegetation history and climate
The postglacial pollen sequences from southern Quebec are generally divided into three main stages corresponding to the postglacial development of plant communities after glacial retreat: the non-arboreal stage, the afforestation stage and finally the forest stage (e.g.Fréchette et al., 2021;Richard, 1994).The dominance of tree taxa over the entire P42 pollen record may suggest forest vegetation in surrounding lands from at least 8000 years ago, which is consistent with the regional palynostratigraphy of coastal sites around the GSL indicating forested landscape from ~8000 cal years BP (Fréchette et al., 2021).The key taxa of the preceding afforestation stage (Alnus crispa, Populus and Juniperus) do not record comparable or significant percentages in core P42 (cf.Fréchette et al., 2021).These features indicate that our marine record encompasses only the forest stage.
During the early forest stage (>7360 cal years BP; zone PI; Figure 4), pollen and spore concentrations rapidly increased together with percentages of Picea, suggesting the development of dense forests on adjacent lands.An interesting feature in our record is the decline of Tsuga relative to Abies at around 5500-5000 cal years BP, which likely corresponds to the Tsuga regional decline (Bennett and Fuller, 2002;Fuller, 1998).After 3500 cal years BP, Picea dominated the assemblages and Abies slightly increased, highlighting the beginning of the closed boreal forest phase.The last 1260 years were marked by a drastic decrease in Abies percentages (9 to 0%) associated with a gradual increase in Alnus and Salix percentages.An increase of Alnus 1000 years ago is also recorded at two sites from the north shore of the Gaspé Peninsula (Lac J'Arrive, Lac à Euloge; Marcoux and Richard, 1995).In our data set, the particularly pronounced increase of Salix and Alnus around 560 cal years BP probably reflects an opening of the forest cover in the source area, which might be the regional signature of the Little Ice Age cooling.This event coincides with low SST in our record, which may be the result of enhanced atmospheric circulation from the North.Finally, Ambrosia pollen occurrence after ca.270 cal years BP likely relates to the spreading of ruderal plants associated with anthropic activities.
The timing of the major vegetation changes is asynchronous across the GSL territory, and the vegetation composition also differed between sectors (Fréchette et al., 2021).Our pollen record seems to display a largely mixed signal across the area.Compared to pollen spectra from lakes and peats, our marine record as others from the GSL (e.g.Jessen et al., 2011;Levac, 2003;Wu et al., 2022) shows an over-representation of bisaccate pollen grains, especially Pinus.In the topmost sediment of lacustrine cores from terrestrial sites around the GSL, Pinus pollen records less than 10% (Fréchette et al., 2021) whereas it is about 30% in core P42.Pinus is a thermophilic species related to atmospheric long-distance transport from southern mixed forests to the GSL area (Jessen et al., 2011;Wu et al., 2022).Among bisaccate pollen, that of Picea is also overrepresented in marine sediments but to a lower extent.Picea percentages in recent lake sediments are ~35% on the St. Lawrence North Shore and 15%-30% on the Gaspé Peninsula (Fréchette et al., 2021), whereas the topmost assemblage of core P42 records ~50% of Picea.In contrast, Betula pollen exceeds 40% in recent lake sediments from the GSL area while it averages less than 20% in core P42.Because of such discrepancies between the studied marine core and terrestrial sites around the GSL, the pollen record from core P42 is not used to reconstruct quantitatively the regional climate.
Although we cannot reconstruct climate parameters on a quantitative basis, the gradual increase of Picea relative to Pinus in our record indicates increasing proportions of inputs from the boreal forests during the last 8000 years.This could be explained by the expansion of Picea-dominated boreal forest on the North Shore and the Gaspé Peninsula Plateau (e.g.Fréchette et al., 2021), or the decreasing intensity of southerly winds bringing pollen from mixed and temperate forests (Jessen et al., 2011;Wu et al., 2022).Hence, the trend of increased Picea pollen abundance relative to Pinus suggests a cooling on adjacent lands, which is consistent with the pollen-based reconstructions of the mean temperature of the warmest month (MTWA) from nine terrestrial sites around the GSL (Fréchette et al., 2021; see Figure 9).

Palaeoceanography and regional climate change during the Middle Holocene
The dinocyst assemblages of core P42 and others from the GSL (e.g. de Vernal et al., 1993) suggest generally cool-temperate surface waters (e.g.Marret and Zonneveld, 2003;Zonneveld et al., 2013).Prior to ca. 7300 cal years BP, the assemblages feature the lowest percentages of N. labyrinthus and I. minutum, highest percentages of S. quanta and P. americanum, relatively high percentages of Brigantedinium spp., as well as rapidly increasing percentages of cyst of P. dalei.This translates into a possibly eutrophic estuarine environment with low salinities due to river discharge and a short sea-ice season (e.g.Head et al., 2001;Marret and Zonneveld, 2003;Potvin et al., 2013;Zonneveld et al., 2013).The dinocyst-based reconstruction indicates low temperatures (<10°C) and salinities (<28 psu), which we associate with early postglacial environments marked by high freshwater/meltwater discharges.This is also consistent with the decrease of δ 13 C org in the sediment that indicates an increase of terrestrial input of organic matter.
From about 7300 to 4200 cal years BP, the dinocyst assemblages record a transition in the percentages of several heterotrophic taxa.It is marked by the decrease of Brigantedinium spp., a taxon that often characterizes upwelling systems and high productivity in all thermic domains (e.g. de Vernal and Marret, 2007), the decrease of S. quanta, a taxon able to inhabit a wide range of environments from warm to cool, and the increase of I. minutum, which is typical of cold waters and sea-ice cover (de Vernal et al., 2005(de Vernal et al., , 2013;;Head et al., 2001;Potvin et al., 2013).The variations  (Fréchette et al., 2021).(e) Water table depth (WTD) reconstruction from two maritime peatlands (PLU and PTE) on Anticosti Island (Perrier et al., 2022).(f ) δ 18 O calcite (‰ VPDB) of B. marginata (this study).(g) Downcore variation of the PC1 of the assemblages of pollen (presented with reversed sign), dinocyst and foraminifer in core P42 (this study).LIS: Laurentide Ice Sheet.DO: dissolved oxygen.LCW: Labrador Current Water.
in dinocyst assemblages reflect a change in trophic conditions as illustrated by estimated productivity showing a decrease from 1037 to 513 g C m −2 year −1 .The reconstruction also demonstrates a rapid warming in surface water during summer, from 9.6 to 11.8°C between 7300 to 6600 cal years BP, and summer SST oscillating around 11°C afterwards.
The SST increase after ~7200 cal years BP occurred simultaneously with the increase of Picea percentages in core P42, but more than 2000 years after the air temperature increase (MTWA) as reconstructed from pollen data in the GSL region (Fréchette et al., 2021; see Figure 9).The onset of summer warming in surface waters even follows the regional Holocene thermal optimum marked by a maximum MTWA peak at ~7200 cal years BP (Fréchette et al., 2021; see Figure 9).This demonstrates a lag between regional climate over mainland areas and sea-surface conditions.A similar lag was recorded in maritime peatlands from Anticosti Island located ~400 km downstream from our core site indicating optimal conditions were reached as late as 5000 cal years BP (Perrier et al., 2022; see Figure 9).
From the base of the record (~8200 cal years BP) to about 6600 cal years BP, the salinity estimates show an increase parallel to SST, which suggests a close relationship between freshwater discharge and SST.Variations in freshwater discharge could correspond to a shift in climate regime with a change in precipitation over the St. Lawrence watershed.However, variation in salinity could also relate to changes in meltwater discharge rates from the residual Laurentide ice dome that finally subsided about 6000 years ago (Dyke, 2004).Because the pollen-based reconstructions of annual precipitation (PANN) from terrestrial sites indicate very little variation (Fréchette et al., 2021; see Figure 9), we discard the possibility of major changes in precipitation.Hence, our results rather indicate that ice meltwater discharge has been the main driver of low salinities in the early part of the Middle Holocene until ~6500 cal years BP.On this basis, taking into consideration the high insolation in the Early-Middle Holocene and the general thermal optimum peaking at ~7600-5000 cal years BP globally (Kaufman et al., 2020) and about 7200 cal years BP regionally (Fréchette et al., 2021; Figure 9), we hypothesize that warm conditions leading to an acceleration of ice melt during the final phase of the deglaciation resulted in enhanced cold freshwater discharges, which in turn prevented warming in the surface waters.After ca.6600 cal years BP, the summer salinity stabilized around 28.5 psu, which suggests limited meltwater input after the final ice sheet withdrawal.
The interval prior to ~6500 cal years BP is also marked by low concentrations of foraminifer tests despite abundant organic linings and high percentages of Quinqueloculina seminulum and Pyrgo williamsoni, both characterized by a thick and resistant test (Kaiho, 1994).Hence, dissolution probably has occurred in the sediment, perhaps as a response to the contribution of corrosive meltwaters in the bottom water masses.High primary productivity and C org fluxes to bottom waters resulting in high CO 2 release in the sediment, as already documented from δ 13 C records (Figure 8), may also explain the poor CaCO 3 preservation.

The Middle to Late-Holocene transition at ~4200 cal years BP
An important transition marks the dinocyst assemblages of core P42 at ~4200 cal years BP, which corresponds to the boundary between the Middle and Late-Holocene.This transition is highlighted by an abrupt drop in PC2 scores, which mainly corresponds to the almost disappearance of P. americanum.Together with the long-term decline of S. quanta and Brigantedinium spp.across the transition, the dinocyst assemblages could be translated into an important shift in trophic conditions (e.g.Bringué et al., 2014;Marret and Zonneveld, 2003;Zonneveld et al., 2013).This transition is also marked by estimates of sea-surface conditions.Around 4200 cal years BP, SSS increased by 2.5 psu and productivity decreased by 486 g C m −2 year −1 .Such variations reflect a drastic decrease in freshwater and nutrient inputs, corresponding to a change from estuarine to more oceanic conditions.The antiphase variation of reconstructed SSS and primary productivity at this transition, as well as over the entire record, suggests that freshwater and nutrient input could have been a main limiting factor of primary production in the LSLE.
At the transition from the Middle to Late-Holocene, there is also an increase in the δ 18 O of B. marginata followed by an important change in benthic foraminifer assemblages.There is an increase in the Labrador Current indicator N. labradorica and the cold-water species E. exavatum, while there is a decrease in the North Atlantic water indicator O. umbonatus and a warm water species E. exilis.These changes indicate a shift towards cool bottom waters, which could be due to increased contribution of the LCW relative to NACW and/or cooling of the NACW, a feature also documented by δ 18 O values in B. marginata (Figure 8).This points to a decoupling between surface water conditions, which are controlled by regional hydroclimate, and bottom water properties, which are driven by the large-scale ocean circulation patterns (Gilbert et al., 2005).Hence, the Middle to Late-Holocene bottom water cooling of the Laurentian Channel is consistent with the large-scale climate trends, notably in the northwest North Atlantic (e.g.Kaufman et al., 2020;Marcott et al., 2013;Sachs, 2007) and/or with a gradual decrease in the strength of the North Atlantic Current relative to the Labrador Current along the eastern Canadian margins.After ~4200 cal years BP, the relatively high SSS and SST in the LSLE reflect the establishment of modern-like hydroclimatic conditions in southeastern Canada.

The Late-Holocene
The Late-Holocene part of the record is marked by subtle changes in surface water conditions as suggested by the dinocyst assemblages.Most taxa experienced fluctuations without any noticeable long-term trend, except for P. americanum that increased markedly between ca.1200 and 500 cal years BP, which corresponds to an increase in reconstructed annual primary productivity.Estimated SST generally fluctuated between 10 and 12 °C.From ca. 4100 to 1500 cal years BP, SSS fluctuated inversely to SST, suggesting a change in the hydroclimatic regime that controls surface hydrography in the LSLE.According to our reconstruction of SST and SSS, oscillations between colder and drier (low runoff) versus warmer and wetter (high runoff) climates were recorded regionally.
The last millennium is characterized by a cool interval between 950 and 490 cal years BP accompanied by low SSS, suggesting possibly higher precipitation at the scale of the drainage basin, in agreement with a slight increase in the PANN reconstructed by Fréchette et al. (2021) (see Figure 9).This pulse corresponds to a maximum in productivity, which is consistent with dinocyst assemblages and the δ 13 C of the organic matter.After 430 cal years BP, a shift towards decreased bottom water oxygenation is shown by the occurrence of Brizalina subaenariensis, a species related to low dissolved oxygen concentration (Kaiho, 1994;Sen Gupta, 2003).This trend is associated with warming in bottom waters, highlighted by the increase of E. exilis, responsible for higher respiration rates and oxygen consumption together with the enhanced contribution of the oxygen poor Atlantic waters relative to the oxygen rich Labrador waters (Genovesi et al., 2011;Gilbert et al., 2005;Thibodeau et al., 2013).The isotopic record of benthic foraminifera in the study core, with decreasing δ 18 O and δ 13 C values, indicating a recent warming as well as reduced ventilation, provides further evidence for a significant change in the properties of bottom waters at the scale of the Laurentian Channel.

Conclusions
The paleoceanographic and palaeoclimatic data presented here show important changes during the Middle and Late-Holocene.They illustrate both large-scale ocean circulation patterns and sea-surface conditions related to the regional hydroclimatic regime, as well as the succession of postglacial forest on adjacent lands.The pollen data from core P42 record a mixed regional signal of long-term changes with gradual increase of Picea relative to Pinus, which would reflect the transition towards boreal forest and/or change in the atmospheric circulation with increased eolian inputs from the North relative to the South, both being related to a large-scale cooling trend.
Benthic foraminiferal assemblages and the isotopic composition of their tests permit us to infer a long-term trend of decreasing temperatures in bottom waters, which we tentatively associate with enhanced proportions of LCW entering the GSL.Superimposed on these trends that reflect large scale processes, higher frequency variations are also recorded, probably reflecting regional climate or ocean oscillations.Among these, the most spectacular event is recorded by dinocyst assemblages that show an abrupt transition from estuarine to oceanic conditions in surface waters at 4200 cal years BP.This transition is marked by a significant increase in salinity concomitant with a decrease in primary productivity.It possibly represents a tipping point in regional climate with decreased precipitation at the scale of the St. Lawrence drainage basin and/or a reorganization of regional ocean circulation.The fact that the transition coincides with the Middle to Late-Holocene boundary, which marks the end of the Holocene Thermal Optimum (e.g.Kaufman et al., 2020;Marcott et al., 2013), the full establishment of Labrador Sea Water formation (Hillaire-Marcel et al., 2001) and a shift towards aridification at low latitudes (e.g.DeMenocal et al., 2000) is interesting and suggests a major reorganization in hydroclimatic conditions related to ocean circulation.
Another important feature in core P42 deserves attention.Benthic foraminiferal assemblages and the isotopic composition of benthic foraminifer tests indicate a warming accompanied by decreasing concentration of dissolved oxygen in bottom waters over the last centuries, which appear independent from the surface variations.Hence, the data from the LSLE area illustrate a complex response of ecosystems to climate changes driven by global forcing and regional mechanisms.

Figure 1 .
Figure1.The St. Lawrence marine system and location of piston core COR0602-42PC (49°7.18′N -67°16.80′W).Other study sites mentioned in this article are marked by white squares (see Table1for coordinates and references).

Figure 2 .
Figure2.Left: age versus depth relationship in core P42 using the Bacon software ofBlaauw and Christen (2011) and the 'Marine20' calibration curve(Heaton et al., 2020), with a regional offset (ΔR) of 25 ± 72 years (see text).Semi-transparent blue zones represent calibrated age ranges as in Table1.The red dashed line represents the median values of the age model, and the grey dashed lines represent 95% confidence ranges.Right: dry bulk density (g cm −3 ) versus depth, calculated as the ratio of the mass of dried sediment to its total wet volume.
Note.Please refer to the online version of the article to view this figure in color.

Figure 4 .
Figure 4. (a) Relative abundance (%) of the main taxa (>1%) in pollen and spore assemblages, total concentrations (g −1 ) of pollen grains and spores, and curves of principal components PC1 (bold) and PC2 (dashed) in core P42.Pollen percentages and spore percentages are calculated on the total pollen sum.(b) Ordination diagram of the pollen taxa.

Figure 6 .
Figure 6.Sea-surface condition estimates from the Modern Analogue Technique (MAT) applied to dinocyst assemblages in core P42 (grey curves).The bold curves correspond to 5-point running means.The short horizontal lines represent the maximum and minimum distance of the set of five analogues for each sample.Modern averages of sea-surface parameters are indicated by vertical grey lines and horizontal grey zones on the upper axis.

Figure 7 .
Figure 7. (a) Relative abundance (%) of the main benthic foraminifer taxa (⩾3% in at least one sample), total concentration (tests g −1 ) of calcareous taxa, foraminifer lining to test ratio, and curves of PC1 (bold) and PC2 (dashed) in core P42.LC and NAC stand for Labrador Current and North Atlantic Current respectively.White zones represent samples with too few foraminifer counts (n < 20) to be statistically interpreted.(b) Ordination diagram of benthic foraminifer taxa (n = 15; ⩾3% in at least one sample).

Figure 8 .
Figure 8. Top: carbon isotope composition of Bulimina marginata tests versus that of the organic carbon in core P42; Bottom: comparison of oxygen isotope composition of benthic foraminifera at the study site (P42), versus those from site MD99-2220 (from de Vernal et al., 2023), upstream the St. Lawrence Estuary (SLE), and from MSM45-31-1 (from Lochte et al., 2020), on the Labrador Shelf; the dashed vertical line illustrates the changing trend towards a recent warming of the deep water mass filling the Laurentian Channel (see text).

Figure 9 .
Figure 9. (a) SSS and (b) SST reconstruction from core P42 presented as anomalies against the modern average value (this study).(c) Annual precipitation (PANN) and (d) mean temperature of the warmest month (MTWA) reconstruction from nine terrestrial sites around the GSL presented as anomalies against the modern average value(Fréchette et al., 2021).(e) Water table depth (WTD) reconstruction from two maritime peatlands (PLU and PTE) on Anticosti Island(Perrier et al., 2022).(f ) δ 18 O calcite (‰ VPDB) of B. marginata (this study).(g) Downcore variation of the PC1 of the assemblages of pollen (presented with reversed sign), dinocyst and foraminifer in core P42 (this study).LIS: Laurentide Ice Sheet.DO: dissolved oxygen.LCW: Labrador Current Water.

Table 1 .
Van Nieuwenhove et al. (2020)al.(2020).Results are presented in percentage diagrams, with certain taxa combined.For example, Spiniferites spp.includes Spiniferites Summary of regional palaeoclimatic studies cited in this paper. .