Benthic Foraminiferal Assemblages from the Laurentian Channel in the Lower Estuary and Gulf of ST. Lawrence, Eastern Canada: Tracers of Bottom-Water Hypoxia

Historical measurements of dissolved oxygen (DO) concentrations in the bottom waters of the Lower St. Lawrence Estuary (LSLE) reveal a decrease of more than 50%, from ∼125 μM L^–1 in the 1930s to an average of ∼60 (±6) μM L^–1 between 1985 and 2018 (Gilbert et al., 2005; Jutras et al., 2020a). This impacted benthic marine ecosystems (Belley et al., 2010) with deleterious effects on the marine biota, including reduced growth and reproduction of commercially exploited marine species such as cod and northern shrimp (D'Amours, 1993; Chabot & Dutil, 1999; Dupont-Prinet et al., 2013). Furthermore, at hypoxic levels (<62.5 µM L^–1), fish populations are increasingly vulnerable to additional stressors, natural and anthropogenic (Breitburg, 2002).

Oxygen depletion in the deep waters of the LSLE and the Gulf of St. Lawrence (GSL) has been attributed to various factors, including changes in the northwest Atlantic Ocean circulation pattern (Gilbert et al., 2005; Thibodeau et al., 2018; Jutras et al., 2020a), increased organic matter and nutrient fluxes (i.e., eutrophication) in the LSLE (Thibodeau et al., 2006, 2010; Jutras et al., 2020a, b), as well as increased bottom-water temperatures and the associated enhanced microbial respiration rates (Genovesi et al., 2011). In a recent study, Jutras et al. (2020a) determined that the dominant causes of oxygen depletion in the bottom waters of the Laurentian Channel (LC) may have changed over the last 50 years: the impact of enhanced microbial respiration in response to warmer bottom-water temperatures and eutrophication was likely determinant between the 1970s and the late 1990s, while a lower contribution of cold, oxygen-rich Labrador Current Water (LCW), relative to the warmer, oxygen-depleted North Atlantic Central Water (NACW) may be invoked between 2008 and 2018.

Regardless of the causes of oxygen depletion, the decrease in bottom-water oxygenation in the LC over the past century resulted in variations in the benthic foraminiferal assemblages preserved in the sediment (Thibodeau et al., 2006; Genovesi et al., 2011). As benthic foraminifera are sensitive to oxygenation (Sen Gupta & Machain-Castillo, 1993; Bernhard et al., 1997; Osterman, 2003; Platon et al., 2005), the abundance and diversity of species can be used to qualitatively and quantitatively (Ohkushi et al., 2013; Tetard et al., 2017, 2021a, b; Hoogakker et al., 2018; Erdem et al., 2020) assess the level of bottom-water DO concentrations. Temperature, primary production, organic carbon fluxes, and water mass distribution may also modulate the distribution of assemblages (e.g., Jorissen et al., 2007). Previous studies by Hooper (1975) and Rodrigues & Hooper (1982) provided a comprehensive overview of the benthic foraminiferal distribution in sediment collected in the LSLE and GSL in the 1970s. A few other studies documented recent changes in benthic foraminiferal assemblages at a few sites in the LSLE (Thibodeau et al., 2006, 2010) and GSL (Genovesi et al., 2011).

In this study, we analyzed the geochemical composition and micropaleontological content of four sediment cores collected in June 2018 in the LC of the LSLE and GSL. The records encompass decades to a few centuries. The two main objectives were: (1) to document how the recent and sub-recent benthic foraminiferal assemblages changed with the temporal and spatial evolution of the hypoxic zone across the LC, and (2) to semi-quantitatively evaluate the rates of changes in benthic foraminiferal assemblages in relation to bottom-water oxygenation.

The St. Lawrence Estuary system on the eastern coast of Canada is considered the largest semi-enclosed estuary in the world (Fig. 1). The freshwater discharge at Quebec City is estimated at 11,000 m³ s^–1, second only to the Mississippi in North America (Bourgault & Koutinonsky, 1999). Its main bathymetric feature is the Laurentian Channel (LC), a 300–500 m deep, U-shaped submarine valley that extends ∼1280 km landward from the continental shelf break to the head of the Lower St. Lawrence Estuary (LSLE) at Tadoussac (Figs. 1, 2). The LSLE and the Gulf of St. Lawrence (GSL) are characterized by estuarine circulation. During most of the sea ice-free season, the water column consists of a three-layer system based on its thermal stratification. The surface layer (0–30 m) flows seaward and results from the mixing of water from the northwest North Atlantic and freshwater from the St. Lawrence River and tributaries along its north shore (Dinauer & Mucci, 2018). It is characterized by relatively low salinities (S_P = ∼25–31, where S_P is the practical salinity) and seasonally variable temperatures. It overlies the Cold Intermediate Layer (CIL), a cold (T = –1 to –2°C) and saline (S_P = ∼32 to 33) water mass that forms in winter in the gulf and flows landward at depths between 50 and ∼150 m (Gilbert & Pettigrew, 1997; Galbraith, 2006). The deep waters below ∼150 m originate from a mixture of warm, oxygen-poor North Atlantic Central Water (NACW) and cold, oxygen-rich Labrador Current Water (LCW), that enter the GSL through Cabot Strait (Dickie & Trites, 1983; Gilbert et al., 2005; Fig. 2B). They are relatively warm (T = 4–7°C) and salty (S_P = ∼33–35). In winter, cooling and decreased freshwater input result in vertical mixing of the surface layer and CIL (Galbraith, 2006; Smith et al., 2006; Saucier et al., 2009). The surface waters are in contact with the atmosphere and, thus, well-oxygenated. However, the deep waters, which are isolated from the atmosphere by the persistent stratification, progressively lose DO through microbial respiration and remineralization of organic matter during their 4 to 7-year landward journey from the continental shelf break to the head of the LC (Bugden, 1991; Gilbert, 2004; Benoit et al., 2006). Consequently, highly oxygen-depleted bottom waters are observed throughout the LSLE and the western GSL (Gilbert et al., 2005; Jutras et al., 2020a, 2022). Any change in the ocean circulation that modifies the relative contributions of the LCW and NACW to the bottom waters of the LC, in addition to variations in organic matter fluxes and temperature, may play a role in modulating oxygenation conditions in the LSLE and GSL.

Four sediment cores were recovered in the LSLE and GSL in June 2018 aboard the R/V Coriolis II (Fig. 1; Table 1) using an Ocean Instrument Mark II box corer. Each box core (0.12 m² × 0.5 m high) was transferred to a specially designed stainless-steel glove box (Edenborn et al., 1987) for subsampling. Stations 23, 21, 18.5, and 16, at which the box cores were recovered, are part of a network of recurring stations along the LC that have been visited regularly over the past 40 years. Each core was subsampled at 0.5-cm increments for the first centimeter. The core from station 23 was subsampled at 1-cm increments from 1 to 5 cm, 2-cm increments from 5 to 15 cm, and 3- or 4-cm increments to the bottom of the core. Cores from stations 21, 18.5, and 16 were subsampled at 1-cm increments to the bottom of the cores. At each depth interval, sediment samples were transferred and stored in two pre-weighed plastic vials. One set of vials was later weighed, freeze-dried, and weighed again to determine the water content. Sediment porosity was calculated from the water content, the salinity of the overlying seawater, and a dry sediment density of 2.65 g cm^–3. The freeze-dried sediments were then ground and homogenized with an agate pestle and mortar in preparation for later analyses (e.g., organic carbon). A known volume (6 to 8 mL) of a Rose Bengal solution (2 g L^–1 in ethanol, density = 0.79 g mL^–1) was added to each of another set of vials to stain the living foraminifera. All subsamples we analyzed were stained, except for those below 25 cm from station 23. These vials were subsequently weighed, and the sediment porosity was used to determine the dry-sediment weight. The exact location, water depth, overlying water salinity, and DO concentration at each sampling station are reported in Table 1.

The chronological framework was established from ²¹⁰Pb measurements for all studied cores, and from both ²¹⁰Pb and ¹³⁷Cs measurements for station 23. The ²¹⁰Pb activity in sediment was obtained indirectly by measuring the decay rate of its daughter isotope ²¹⁰Pb (t_1/2 = 138.4 days; α = 5.30 MeV) by alpha spectrometry assuming secular equilibrium. Following the protocol proposed by Flynn (1968), ²⁰⁹Po was added to the freeze-dried samples as a spike to determine the extraction and counting efficiency and to measure the decay rate. The samples were submitted to a series of acid digestions (HCl, HNO₃, HF, and H₂O₂), and the purified polonium was deposited on a silver disk. The ^209–210Po activities were measured in a silicon surface-barrier multichannel α-spectrometer (EGG-ORTEC, type 576A) coupled to a MAESTRO Multichannel Analyzer Emulation Software. Uncertainties were estimated from counting statistics to be approximately 6% for two standard deviations (2σ).

Sedimentation rates were estimated using the ²¹⁰Pb-decay constant (λ), the linear least-squares regression slope of the logarithmic function of excess ²¹⁰Pb and the CF (CRS) model with Monte Carlo uncertainty under the assumption of a constant rate of supply. The Plum software (Aquino-López et al., 2018) under R (R Core Team, 2021) was used to define the age models for each core (Fig. 3), taking into consideration the total (²¹⁰Pb_tot) and supported (²¹⁰Pb_sup) ²¹⁰Pb activities, the density (g cm^–3) of samples, and the coring date.

The supported fraction of ²¹⁰Pb was determined in the bottom section of the cores using the activity of the parent isotope ²²⁶Ra. To confirm the validity of the chronological framework, we measured the vertical distribution of ¹³⁷Cs activity in the sediment core recovered from station 23. The ¹³⁷Cs activity was measured at 661.6 keV and the ²²⁶Ra activity was measured at 295.22, 352, and 609.32 keV using a Canberra low-background high-purity Ge well-detector. The reproducibility of ²²⁶Ra and ¹³⁷Cs activity measurements was estimated at ±1% (1σ) based on replicate analyses (n = 6) of the Columbia River Basalt reference material (BCR-2; USGS).

The standard deviation (1σ) on the C_tot and N_tot analyses was estimated at, respectively, 0.1% and 0.3%, based on replicate measurements of organic analytical standard substances (acetanilide, atropine, cyclohexanone-2,4-dinitrophenyl-hydrazone, and urea).

For isotopic analyses of C_org, freeze-dried, ground and homogenized sediment samples were fumed in a desiccator for 48 hours in the presence of concentrated hydrochloric acid to dissolve carbonate phases. The isotopic composition of C_org (δ¹³C_org) was measured with a Micromass Isoprime^TM ratio mass spectrometer in-line with an Elementar Vario MicroCube^TM elemental analyzer in continuous flow mode. Data are reported in the δ notation in ‰ with respect to V-PDB (Coplen, 1995; Fig. 4). As determined from replicate measurements of standard materials, the analytical uncertainty was estimated at ±0.1‰. To normalize the results on the NBS19-LSVEC scale, two internal reference materials (δ¹³C = –28.73 ± 0.06‰ and δ¹³C = –11.85 ± 0.04‰) were used, and a third reference material (δ¹³C = –17.04 ± 0.11‰) was analyzed as an unknown to assess the accuracy of the normalization.

Each sediment sample was treated according to the procedure described by de Vernal et al. (1999). Depending on the liquid content of the subsamples, 5 to 13 cm³ of wet sediment were weighed, dried at room temperature, weighed again, and wet-sieved through 63 μm and 106 μm mesh sieves. Both these fractions were examined under a binocular microscope at 40× and 100× magnification. Calcareous and agglutinated benthic foraminifera were hand-sorted, enumerated, identified, and counted according to the nomenclature of Rodrigues (1980) and Loeblich & Tappan (1987). Planktonic foraminifera were not counted, as they only occurred occasionally. Well-preserved specimens of calcareous benthic foraminifera, photographed using a Keyence^TM microscope, are presented in Figures A1 and A2. The list of all observed taxa and their relative abundance (e.g., rare, occasional, common) in the LSLE and GSL sediments are reported in Table A1.

In oxygen-depleted environments, the protoplasm of dead foraminifera may be preserved for several months (Bernhard, 1988; Corliss & Emerson, 1990), which could lead to an overestimation of the abundance of living individuals from Rose Bengal staining (Fontanier et al., 2002). Hence, we applied strict staining criteria, as described in Fontanier et al. (2002), to overcome the inaccuracies of the method, by counting individuals as living only if the Rose Bengal stained all chambers except the last one. Partly stained tests were still noted as an indicator of the presence of an undecayed protoplasm.

The relative abundance of various taxa, expressed as percentages, were calculated for all taxa but only common and occasional taxa are reported. The relative abundance of taxa that accounts for ≥3% of the assemblages in at least one sample per core are presented in Figure 5.

Multivariate analyses using the CANOCO 5 software (ter Braak & Šmilauer, 2012) were conducted with the percentages of calcareous taxa to help define zonation and assemblages. We performed principal component analyses (PCA) and illustrated the scores of the two main principal components (first axis PC1; and second axis, PC2) on ordination diagrams available in the Supplementary Materials. We also performed a constrained redundancy analysis (RDA) based on the calcareous assemblage at station 23 and temperature and DO as environmental variables. Finally, we calculated dissimilarity coefficients using the squared chord distance (SCD) to define the depth of the most pronounced changes in assemblages (Overpeck et al., 1985). The PCA, RDA, and multivariate SCD analyses were conducted on taxa that account for ≥3% of the assemblage in at least one sample to reduce crowding and the effect of rare taxa on the analyses.

Uniform values of the excess ²¹⁰Pb activity (²¹⁰Pb_ex) in the upper part of the cores are interpreted as delimiting the depth of the mixed layer (Fig. 3). At station 23, the ln(²¹⁰Pb_ex) data can be fit to a simple linear equation with a slope of −0.05 (R² = 0.98), corresponding to a sedimentation rate of ∼0.70 ± 0.03 cm a^–1. At station 21, the linear fit to the ln(²¹⁰Pb_ex) data yields a slope of −0.08 (R² = 0.96), corresponding to a sedimentation rate of ∼0.42 ± 0.03 cm a^–1 (Fig. 3B). At station 18.5, the linear fit to the ln(²¹⁰Pb_ex) data yields a slope of −0.12 (R² = 0.96), corresponding to a sedimentation rate of ∼0.27 ± 0.03 cm a^–1 (Fig. 3C). At station 16, the results suggest that there may have been a change in sedimentation rates. The first slope over the upper 5 cm is −0.10 (R² = 0.96), corresponding to a sedimentation rate of ∼0.31 ± 0.03 cm a^–1. The second slope over the 5−20 cm interval is −0.04 (R² = 0.97), corresponding to a sedimentation rate of ∼0.11 ± 0.03 cm a^–1 (Fig. 3D). Based on ²¹⁰Pb and ¹³⁷Cs measurements, Genovesi et al. (2011) also concluded that sedimentation rates in this area of the GSL may have varied through time over the previous two centuries. The sedimentation rates obtained in this study are consistent with estimates from Thibodeau et al. (2006) at station 23 and Genovesi et al. (2011) near station 16.

Since none of the ²¹⁰Pb_ex profiles follow a decreasing exponential curve under the mixed layers, we assume a constant flux of ²¹⁰Pb to set a chronology. Our chronological framework established with the Plum model was confirmed by measurements of ¹³⁷Cs activities throughout the core recovered at station 23. Assuming the location of the ¹³⁷Cs peak at a sub-bottom depth of ∼25 cm (not shown) corresponds to 1986 CE, the year of the Chernobyl incident, our ²¹⁰Pb-based age model at this station is accurate. Furthermore, the correlation of the benthic foraminiferal species Brizalina subaenariensis and Eubuliminella exilis between our station 23 core and core CR02-23 investigated by Thibodeau et al. (2006), taken at the same location, indicates that the bottom of our core (41 cm) is not older than ∼1960, as the two species first appeared around this time in CR020-23.

According to our age models and the above-mentioned correlation, we can assess that the cores we recovered from station 23 encompass the last ∼50 years of sedimentation (i.e., ∼1965–2018 CE). Thus, the geochemical and micropaleontological analyses performed at 0.5- and 1-cm intervals provide time series with sub-annual to decadal resolution. The Plum settings and output for each station are available in the Supplementary Materials (Fig. S1).

The organic carbon content of the sediment in the four study cores ranges from 1.26% to 1.91%. There is a slight increase in organic carbon content from the LSLE to the GSL, with an average of 1.55% at station 23, 1.53% at station 21, 1.61% at station 18.5, and 1.75% at station 16. However, these values do not necessarily translate into higher organic carbon fluxes and productivity as they may reflect a lower dilution by organic material inputs of terrestrial origin, as would suggest the δ¹³C_org data. The organic carbon fluxes can be computed based on the estimated sedimentation rates (Fig. 4). The average total mass accumulation rates are 1.54, 0.99, 0.68, and 0.25 g cm^–2 a^–1 at stations 23, 21, 18.5, and 16, respectively.

The organic carbon content of all four cores typically decreases slightly with depth, while the C_org:N_tot molar ratio increases slightly (Fig. 4). This may be attributed to the preferential degradation of nitrogen-rich organic matter during early diagenesis (Meyers, 1997; Muzuka & Hillaire-Marcel, 1999).

The δ¹³C_org data show distinct signatures at each study site, with lower values in the LSLE and higher values in the GSL. The mean δ¹³C_org are –24.2‰ at station 23, –23.2‰ at station 21, –22.6‰ at station 18.5, and –22.1‰ at station 16. They reflect the higher proportions of terrestrial versus marine organic matter at stations in the LSLE than in the GSL.

Specimens of calcareous and agglutinated taxa were well preserved, with few broken or altered tests. A total of 52 calcareous and four agglutinated taxa were identified in the >63 µm fraction of the four study cores. Among those, seven calcareous taxa are common to abundant at all sites, representing >10% of the assemblages in many samples: Brizalina subaenariensis, Eubuliminella exilis, Globobulimina auriculata, Bulimina marginata, Nonionellina labradorica, Pullenia osloensis, and Elphidium clavatum. Additionally, nearly 25 taxa are considered occasional, whereas 20 taxa are considered rare as they account for <3% of the assemblages in samples of each core.

We observed a loss in the number of species towards the core top at all stations, indicating a loss of diversity over time (Fig. 5). In the LSLE, at station 23, the number of species decreased by about 60% since the 1960s. In the GSL, at station 16, species diversity dropped by almost 65% over the same time period, while the entire sedimentary record shows a decrease of 75% over the last ∼160 years. The number of species decreased by about 50% over the last ∼30 years at station 21, and by about 60% over the last ∼50 years at station 18.5.

In the LSLE, the overall assemblage is dominated by the calcareous taxa B. subaenariensis, E. exilis, G. auriculata, B. marginata, N. labradorica, and P. osloensis (Fig. 5A, B), which together account for about 80% of the assemblages in cores from stations 23 and 21. Important secondary taxa in the LSLE include Cassidulina reniforme and E. clavatum.

In the GSL, we observe a higher diversity than in the LSLE with assemblages dominated by B. subaenariensis, E. exilis, P. osloensis, and Alabaminella weddellensis which together account for ∼67% of the assemblage at station 18.5. At station 16, they are dominated by E. exilis and P. osloensis, which account for ∼40% of the assemblage. Secondary species include C. reniforme, Elphidium spp. (mostly E. clavatum), G. auriculata, B. marginata, and N. labradorica.

Brizalina subaenariensis is the dominant species in all study cores, except at station 16 where E. exilis dominates. This suggests that B. subaenariensis has a higher tolerance for reduced oxygen levels than E. exilis, since the bottom-water DO concentration is much higher (∼117 µM) at station 16 than at sites in the western part of the LC (<80 µM; Fig. 6).

The PCA results from each core analyzed independently show different ordinations, but the opposition between B. subaenariensis and E. exilis versus most other taxa remains a characteristic feature (see Supplementary Materials, Fig. S2). Three main assemblage zones were defined based on these PCA results, and their boundaries were set using both PC1 and PC2, and maximum SCD coefficients (Fig. 5).

Zone 1 is generally characterized by the dominance of B. subaenariensis, E. exilis, G. auriculata, which are low-oxygen tolerant species (Kaiho, 1994; Sen Gupta, 1999; Brüchert et al., 2000), as well as by a low diversity of species. It encompasses approximately the last 20 years, from the early 2000s to 2018, as derived from the age model in the core from station 23.

Zone 2 is likely a transitional zone, as it corresponds to a decrease in species diversity and the successive decrease or disappearance of most occasional taxa, such as B. marginata, P. osloensis, E. clavatum, as well as an increased abundance of B. subaenariensis, E. exilis, and G. auriculata. It encompasses a period from about ∼1960 to ∼2000 in stations 23 and 16.

Zone 3 is observed in the lower part of the core from station 16. It is characterized by the highest taxonomical diversity, low percentages of B. subaenariensis and G. auriculata, and the common occurrence of E. exilis, P. osloensis, and E. clavatum, as well as the indicator of well-oxygenated waters Quinqueloculina seminulum (Kaiho, 1994).

The significant loss in species diversity and the change towards assemblages dominated by low oxygen-tolerant taxa between Zone 1 and 2 around the late 1990s and early 2000s in all four cores illustrate how rapidly benthic foraminifera responded to the changes in ocean circulation happening at the same time (Jutras et al., 2020a).

The data from station 23, where instrumental measurements of bottom-water conditions are available for an extended period of time (∼90 years), allowed us to perform an RDA to further evaluate the relationship between bottom-water temperatures, DO concentrations, and the downcore measurements of geochemical variables with benthic foraminiferal assemblages (ordination diagram available in the Supplementary Materials; Fig. S3; Table S1). The RDA reveals that oxygenation is the principal environmental parameter affecting the calcareous benthic foraminiferal assemblages in our study area, as it is highly correlated with the RDA axis 1 (R² = 0.9024), which explains 52% of the variance. The relationship with food quality, as illustrated by the C_org:N_tot and δ¹³C_org, also stands out as an important environmental parameter, as suggested by the RDA axis 2 (R² = 0.4225 and R² = 0.3984, respectively), which explains 11% of the variance.

With the DO concentrations highlighted as the principal environmental parameter affecting benthic foraminifera, the results of the PCA in all analyzed samples (see Fig. 7) combined with ecological notes from the literature led us to define three main groups of benthic foraminifera in relation to the oxygenation conditions of their habitat. Ordination diagram shows that PC1, which represents 47% of the variance, is weighted by a strong opposition between B. subaenariensis, E. exilis, and G. auriculata on the positive side and all other species on the other site.

The first group, characterized by the three co-dominant species B. subaenariensis, E. exilis, and G. auriculata, is related to hypoxic bottom-water conditions ([DO] ≤ 62.5 µM) in the LC. This assemblage is typically associated with Zone 1 assemblages at stations in the LSLE. At stations in the GSL, G. auriculata is not a co-dominant species.

The second group is characterized by the presence of Islandiella islandica, A. weddellensis, Gyroidina lamarckiana, P. osloensis, and Bolivinellina pseudopunctata. It corresponds to varying, but on average relatively low bottom-water DO concentrations. Based on the assemblages observed in the core from station 23, we estimate that bottom-water DO concentrations ranged between 62.5 and ∼100 µM.

The third group is characterized by the presence of Paracassidulina neocarinata, Buccella spp., Lagena/Oolina spp., Q. seminulum, Pyrgo williamsoni, Stainforthia concava, and Elphidium spp. It would be an indicator of relatively high bottom-water oxygen levels in the LC. Given the assemblages in the core from station 23, we estimate that bottom-water DO concentrations ranged from ∼100 to ∼140 µM when this assemblage was present at the sediment surface.

It should be noted that these groups are specific to the LC, LSLE, and GSL, and these proxies may not be applicable to the global ocean. Furthermore, the benthic environmental conditions defined by these groups are likely to have changed through time. Although bottom-water oxygenation conditions in the LSLE have been low over the past 40 years, there may have been episodic pulses of higher oxygen concentrations, as suggested by the variability of the areal extent of the hypoxic zone (Jutras et al., 2022) and the co-occurrence of both high and low oxygen-tolerant species in groups 2 and 3.

The foraminiferal assemblages of the LSLE and the GSL include a number of species that are today found in Arctic and subarctic environments (e.g., C. reniforme, I. islandica, Islandiella norcrossi, N. labradorica, S. concava, Valvulineria arctica; Osterman et al., 1999; Lloyd et al., 2007; Jennings et al., 2020; Cage et al., 2021; Ovsepyan et al., 2021). However, we also identified species that are today abundant in temperate to subtropical regions. Noteworthy are Sagrina subspinescens (commonly named Bolivina subspinescens), a relatively warm Atlantic Water species (Rasmussen & Thomsen, 2004), and P. neocarinata, a species found along the western Atlantic margin from the Gulf of Mexico to the coast of Nova Scotia (Cushman, 1922; Phleger & Parker, 1951; all as C. laevigata var. carinata Cushman; Parker, 1954; Vilks & Rashid, 1976: as C. laevigata; Sen Gupta & Aharon, 1994; Platon, 2001; Cage et al., 2021). Moreover, many of the species we observe are known to tolerate low bottom-water oxygenation (B. subaenariensis, E. exilis, G. auriculata, and B. marginata), including P. neocarinata that has also been reported to tolerate the hypoxic conditions that develop seasonally in the Gulf of Mexico (Tichenor, 2013). This complex combination of species may be explained by inputs from three different water masses into the GSL (Fig. 2): the LCW that brings oxygen-rich polar water, the NACW that supplies oxygen-poor warm water, and freshwater from St. Lawrence River that brings nutrients from a large watershed.

In response to the recent changes in the ocean circulation in the northwestern Atlantic Ocean, the consequent modification of benthic environmental conditions in the GSL (Genovesi et al., 2011; Thibodeau et al., 2018; Jutras et al., 2020a), and the development of hypoxic waters (Jutras et al., 2022), benthic foraminiferal assemblages swiftly shifted towards the quasi-exclusive dominance of low-oxygen tolerant or indicator species (i.e., B. subaenariensis, E. exilis, and G. auriculata; Figs. 6, 8) at the expense of the other taxa (Fig. 5). At station 23, where hypoxic conditions (<62.5 µM) have been persistent for nearly 40 years, these three taxa make up 90% of the calcareous assemblage. Although the geographical expansion of the hypoxic zone had only made occasional incursions into the western GSL until 2020 (Jutras et al., 2022), these incursions are reflected in the important decrease of species diversity throughout the LC.

At stations 18.5 and 16 in the GSL, the lack of a long-term instrumental record of DO prevents defining quantitative relationships with foraminiferal assemblages and the relative abundances of B. subaenariensis and E. exilis. Both species are present from the base of the study cores dating back to ∼1860s and their abundance varies little throughout the cores, especially at station 16 (Fig. 5D). The proliferation of E. exilis at station 16 seems to be fostered by temperature, as this species prefers warmer waters (Caralp, 1989; Kaiho, 1994; McKay et al., 2015, 2016; Tetard et al., 2021a). The abundance of B. marginata and P. neocarinata, two other warmer-water species (Murray, 1991; Sen Gupta & Aharon, 1994; Platon, 2001), also increases towards the top of the sequence before disappearing between the late 1990s and the early 2000s at station 16, while the abundance of the colder-water species N. labradorica, P. osloensis, E. clavatum, and V. arctica (Vilks, 1969; Schafer & Cole, 1986; Knudsen et al., 1996; Polyak et al., 2002; Rytter et al., 2002) decreases (Fig. 5D). The increase in warmer-water species at the expense of colder-water species, paired with the disappearance of N. labradorica, a species associated with the LCW (Bilodeau et al., 1994), in the early 1960s in the GSL reflects changes in the relative contribution of parental waters to the LC observed in recent decades, more specifically a lower contribution of the LCW relative to the NACW (Gilbert et al., 2005; Jutras et al., 2020a, 2022).

From a recent historical perspective, the work of Hooper (1975) and Rodrigues & Hooper (1982) is interesting as it documents the distribution of benthic foraminifera in the LC in grab samples collected in the 1960s. Hence, a comparison between their data and ours allows for a further assessment of the evolution of the benthic environment over several decades. The most notable difference between our results and those of Hooper (1975) is the absence of B. subaenariensis in the LSLE and the overall higher taxonomic diversity. At the time of sampling by Hooper (1975) in 1965, bottom-water DO concentrations at the western edge of the GSL were likely in excess of ∼95 µM (see Fig. 6).

In the assemblages from the GSL that were described by Rodrigues & Hooper (1982), the dominance of E. clavatum and S. concava (as Fursenkoina loeblichi; Rodrigues & Hooper, 1982) is a contrasting feature from what we observe in our cores, and in the GSL in general (Genovesi et al., 2011). Elphidium clavatum is an opportunistic species typically associated with low organic carbon fluxes, low salinity, and relatively cold and well-oxygenated waters (Knudsen et al., 1996; Sen Gupta, 1999; Karlsen et al., 2000), although some studies report its occurrence in oxygen-depleted waters (e.g., Sen Gupta et al., 1996). Likewise, Stainforthia has been described as an opportunistic genus associated with high seasonal productivity and cold waters (Alve, 1995; Gustafsson & Nordberg, 2001; Polyak et al., 2002), and S. concava has been associated with Arctic waters (Lloyd et al., 2007). Thus, changes from assemblages dominated by E. clavatum and S. concava to assemblages dominated by E. exilis, as those we observe in zone 1 of all our study cores, could result from the change in the relative contribution of the parental waters and/or the concomitant warming of the bottom waters, both of which impact the benthic ecology, including oxygenation on the sea floor. It is difficult to disentangle these interrelated parameters. Nevertheless, our data and RDA results imply that the DO content is ultimately the primary parameter affecting the foraminiferal assemblages.

As for the environmental interpretation of the three main groups described in the previous section, it should be noted that this relationship is specific to the LC in the LSLE. It might tentatively be applied in the GSL, but not in the global ocean. We thus tentatively applied the relationship to a long-time series from the Esquiman Channel in the eastern GSL (Core CR06-TCE, Fig. 1) to reconstruct the low oxygen event reported to have occurred between 4,000 and 6,000 years ago by Thibodeau et al. (2013). The tentative DO reconstruction presented in Figure 9 highlights hypoxic levels between 4,000 and 6,000 years ago.

The micropaleontological analyses we conducted on four sediment cores from the LC document the multidecadal succession of benthic foraminiferal assemblages and their rapid spatial and temporal evolution in response to changes in benthic environmental conditions in recent decades. The overall foraminiferal assemblages that include both cold- and warm-water species are complex. They likely relate to variations in the contribution of both oxygen-rich polar water from the LWC and oxygen-poor warm water from the NACW.

Whereas the four study sites record a decrease in taxonomic diversity of ∼60–65% across the LC since the 1960s, important changes in benthic foraminiferal assemblages are observed at approximately the same time at all studied sites, around the late 1990s and the early 2000s, coincident with the areal expansion of the hypoxic zone. We observe the decrease or disappearance of most taxa in favor of hypoxia-tolerant indicator species Brizalina subaenariensis, Eubuliminella exilis, and Globobulimina auriculata. These changes reflect how fast benthic foraminiferal assemblages respond to variations of bottom-water DO concentrations across the LC.

Whereas other environmental parameters such as temperature could explain multidecadal changes in the benthic microfauna, DO concentrations stand out as the major determinant parameter. Multivariate analyses highlight a close relationship between the benthic foraminiferal assemblages and bottom-water oxygen concentrations in the LSLE, from which we derived a regional proxy that allows us to estimate changes in bottom-water oxygenation prior to the period of instrumental measurements. This proxy deserves further calibration but appears promising to reconstruct the evolution of bottom-water oxygenation in relationship with large-scale changes in ocean circulation and their impact on the LSLE and the GSL during the present interglacial.

We thank the Fonds de recherche du Quebec Nature et technologies (FRQNT) for their financial support to TA through scholarships and a Strategic Network grant to the Geotop Research Center on the dynamics of Earth system. This project was funded by an NSERC Ship-Time grant to YG and AM as well as NSERC Discovery grants to AM, AdV, CHM and YG. We would also like to acknowledge funding to MSS from the Independent Research Fund Denmark [Grant No. 0135–00165B (GreenShelf)], and from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 869383 (ECOTIP). We acknowledge the National Council of Science and Technology (CONACYT) postdoc grants from Mexico to VCB [CVU no. 174856]. We would like to give special thanks to Gilles Desmeules as well as the captain and crew of the R/V Coriolis II for their dedication and help at sea. Additional thanks to many other colleagues who participated in this study: A. Adamowicz and J-F. Hélie for stable isotope measurements, N. Sanderson for her support with the Plum software, and M. Jutras for her assistance with maps and transects figures. We also thank JFR Editor M. Robinson, Associate Editor B. Wilson, and anonymous reviewers for their constructive comments. Supplementary Materials can be found linked to the online version of this article.

Table S1. RDA axes 1 and 2 scores of calcareous benthic foraminiferal taxa at station 23 with dissolved oxygen (DO), temperature, ratio of organic carbon to total nitrogen (C_org:N_tot), and δ¹³C_org, as environmental variables. RDA axes 1 and 2, respectively, explain 52% and 11% of the variance. Species with the highest and lowest scores in bold. The multivariate analyses were done using the CANOCO 5 software (Ter Braak & Šmilauer, 2012).

Table S2. Record of dissolved oxygen (DO) concentrations at 250 to 355 m of water depth in the LSLE based on data extracted in 2018 from the BioChem database made available by the Department of Fisheries and Oceans Canada, and a compilation of data acquired on the R/V Coriolis II (extended from Gilbert et al., 2005; Jutras et al., 2020a).

Table S3. A comparison of taxa identified by Hooper (1975), Rodrigues & Hooper (1982) and our study in the LC sediments. Large circles correspond to frequent occurrences, small circles to presence, whereas “x” indicate an absence in the sediments. If Hooper (1975), Rodrigues & Hooper (1982) used a different species name than the name used in this study, this former name is given in brackets below the name used by us.

Table S4. Raw counts of all foraminiferal taxa. Counts of fully and partially Rose Bengal-stained individuals are in parentheses besides the total number.

Table S5. Raw counts of individuals sized between 63 and 106 µm.

Figure S1. Plum age-depth models for A station 23 B station 21 C station 18.5 D station 16. The blue rectangles represent the ²¹⁰Pb activity (Bq kg⁻¹). The red line represents the mean model. The grey dashed lines are the 95% confidence intervals. The purple rectangles represent the ²²⁶Ra activity (Bq kg⁻¹). The priors (green lines) and posteriors (grey plots) for each model are shown in the mini plots above the age-depth models. Based on estimated sedimentation rates, the prior parameters corresponding to accumulation rates (acc.mean) were set to 1.4 a cm⁻¹, 2 a cm⁻¹, 3 a cm⁻¹, and 10 a cm⁻¹ for stations 23, 21, 18.5 and 16, respectively.

Figure S2. Ordination diagrams of calcareous benthic foraminiferal taxa (3% in at least one sample) from A Station 23 (n = 13). PCA axes 1 and 2, respectively, explain 61% and 30% of the variance. B Station 21 (n = 11). PCA axes 1 and 2, respectively, explain 60% and 20% of the variance. C Station 18.5 (n = 13). PCA axes 1 and 2, respectively, explain 81% and 7% of the variance. D Station 16 (n = 24). PCA axes 1 and 2, respectively explain 66% and 10% of the variance.

Figure S3. Ordination diagram of RDA of calcareous benthic foraminiferal taxa (n = 19; 3% of the calcareous assemblage in at least one sample), with dissolved oxygen (DO), temperature, ratio of organic carbon to total nitrogen (C_org:N_tot), and δ¹³C_org, as environmental variables at station 23 (see Table S1). RDA axes 1 and 2, respectively, explain 52% and 11% of the variance.