The sediment samples used in this study come from the piston core SIS-249, which measures 1.94 metres in length. This core was retrieved from the lower continental slope of the southern Brazilian continental margin at 2,091 metres below sea level (30°05′ S; 47°05′ W, Figure 1 ). Core SIS-249 was obtained during an oceanographic campaign in the austral spring-summer of 2007 by Fugro Brasil Ltda for the Brazilian National Agency of Petroleum, Natural Gas and Biofuels. Due to the presence of shallower allochthonous sands in the uppermost 48 cm of the core, this study focuses on the carbonate-rich pelagic mud and sandy mud, which lie between 48 and 194 cm. Within this interval, we collected 45 samples with a sampling spacing of 2 to 4 cm to analyse the planktonic foraminifera fossil assemblages and stable isotopes.

To assess planktonic foraminifera assemblages, each sample (~9 cc) was sequentially weighed, washed over a 63 µm sieve, dried at 55°C, and weighed again. To avoid juvenile specimens, which would induce taxonomic biases, planktonic foraminifera were only picked from the >150 µm size fraction (CLIMAP Project Members, 1976; Peeters et al., 1999). The processed samples were then divided with a micro splitter to recover at least 300 non-fragmented planktonic foraminifera tests per sample, as to support significant statistical variations of around 10% for assemblage analyses (Patterson and Fishbein, 1989). The taxonomic classification on the species level followed Schiebel and Hemleben (2017).

The chronology of core SIS-249 was first published by Rodrigues et al. (2018) based on one single Accelerator Mass Spectrometry (AMS) radiocarbon age combined with benthic oxygen stable isotope (δ¹⁸O) stratigraphy. The AMS radiocarbon date from Rodrigues et al. (2018) was measured on tests of Globigerinoides ruber at 58 cm core-depth and the age model developed in the software AnalySeries 1.1 (Paillard et al., 1996).

In the present study, we improved the age model by using six δ¹⁸O tie-points (two points from the new δ¹⁸O _G.rub record and four from δ¹⁸O _Uvig ), which are now correlated with the records from core GL-1090 (Santos et al., 2017a), a nearby core with high-resolution and well-calibrated age model (based on 14 AMS ¹⁴C and 13 stable oxygen isotope correlation points to two reference curves: Lisiecki and Raymo, 2005 and Govin et al., 2014). The new refined age model was created in the R-package “Bacon” v. 2.5.3, which implements Bayesian statistics (Blaauw and Christeny, 2011). We considered a propagated error of 2.5 kyr, conservatively estimated based on the mean accumulation rates of cores SIS-249 and GL-1090 (ca. 1.87 and 0.29 cm/kyr, respectively) and the <2 kyr age error from the GL-1090 reference curve (Santos et al., 2017a). The raw AMS radiocarbon age date from Rodrigues et al. (2018) was calibrated within the R-package “bacon”, using the Marine20 curve (Heaton et al., 2020) and applying a regional reservoir effect (ΔR) of −85 ± 40 years ( Tables S1 and S2 ). This estimate follows the Marine Reservoir Correction Database (http://calib.org/marine/), considering the ages of Nadal De Masi (1999), Angulo et al. (2005), and Alves et al. (2015).

We reconstruct past sea surface productivity using the ratio between the species Globigerina bulloides and Globigerinoides ruber (G.bull/G.rub), the relative abundance (%) of Globigerinita glutinata, and the benthic foraminifera accumulation rate (BFAR).

The G.bull/G.rub ratio is used to reconstruct upwelling events (Conan et al., 2002; Toledo et al., 2008) based on the contrasting ecological preferences of both species. The opportunistic species G. bulloides is associated with eutrophic waters in upwelling zones (Sautter and Thunell, 1991; Peeters et al., 2002; Zaric et al., 2005; Mohtadi et al., 2007; Lessa et al., 2014), while G. ruber, a symbiont-bearing shallow water-dwelling species, is abundant in tropical/subtropical planktonic foraminiferal provinces (Bé and Hutson, 1977; Kučera, 2007; Schiebel and Hemleben, 2017). Moreover, the G. glutinata abundance is typically higher in phytoplankton-rich waters and is used as a proxy for paleoproductivity (Conan and Brummer, 2000; Souto et al., 2011; Pereira et al., 2018).

We also applied the BFAR index, calculated here as the total number of benthic foraminifera multiplied by the sediment accumulation rate as it has been shown to be a reliable proxy for the organic carbon flux to the seafloor in the Brazilian margin (Dias et al., 2021). The BFAR index represents the increments of primary productivity export to the seafloor, in which the increases on benthic biomass are associated to the increasing food availability to the benthic community (Herguera and Berger, 1991; Guichard et al., 1997; Jorissen et al., 2007). Thus, these planktonic and benthic foraminiferal proxies together can indicate changes in the fertilisation mechanisms affecting the photic zone (i.e., biological pump, terrestrial nutrient input) and posterior organic matter export.

To assess changes in sea surface stratification, we studied both the stable isotopic composition of selected planktonic foraminifera species and assemblage counts. As planktonic foraminifera calcitic tests record a mean value corresponding to that of the local water mass properties (i.e., temperature, salinity) at the different depths where they live, their δ¹⁸O signal allows the reconstruction of the seawater conditions at different depth layers (Emiliani, 1954; Ravelo and Hillaire-Marcel, 2007). The species G. ruber dwells at surface-shallow depths, while the species Globorotalia inflata is a subsurface-thermocline dweller (Chiessi et al., 2007; Groeneveld and Chiessi, 2011; Schiebel and Hemleben, 2017; Lessa et al., 2020). Therefore, the gradient between the δ¹⁸O values of these two species (Δδ¹⁸O _{G. inf-G.rub} ) indicates a decreased (increased) upper water column stratification according to the lower (higher) Δδ¹⁸O _{G. inf-G.rub} values (Santos et al., 2017b and references therein).

Around 10 specimens of G. ruber and G. inflata were collected per sample for the δ¹⁸O analyses from the size fraction >250 µm to limit ontogenetic effects ​​(Elderfield et al., 2002). The tests were cleaned with distilled water using an ultrasonic bath to remove any contamination by external particles. The isotopic measurements were performed on a ThermoFisher Scientific MAT253 gas Isotope Ratio Mass Spectrometer coupled to a Kiel IV automated carbonate device at the Research Center for Geochronology and Isotopic Geochemistry (CPGeo) from the University of São Paulo. Isotopic data used the Vienna Pee-Dee-Belemnite (VPDB) reference standard. Here we report the standard deviation of the laboratory reference material used for normalisation (SHP2L; Crivellari et al., 2021), being 0.07‰ (n=20 standards) over the measurement period. The standard deviation of the VPDB δ¹⁸O values of the measured samples did not exceed 0.1‰. As no replicated measurements were carried out, we confirmed the consistency of our δ¹⁸O values and value offset between G. ruber and G. inflata with the expected values for our study site. We verify this by comparing our records with published Holocene (Chiessi et al., 2007) and last glacial and interglacial stage (Santos et al., 2017b) datasets.

In addition, we use the abundances of the Globorotalia truncatulinoides right coiling form, as it has been documented to represent a well-mixed upper water column (reduced stratification), since this morphotype migrates to relatively shallower depths to complete its reproductive cycle (e.g., Lohmann and Schweitzer, 1990 ; Renaud and Schmidt, 2003; Feldmeijer et al., 2015; Billups et al., 2016).

We quantify the effect of dissolution on foraminiferal assemblages by using: i) the ratio between the benthic and planktonic foraminifera (B/P ratio, Arrhenius, 1952; Parker and Berger, 1971; Kučera, 2007), ii) the number of whole planktonic foraminifera tests per gram of dry sediment (PF/g, Le and Shackleton, 1992; Suárez-Ibarra et al., 2021), iii) the abundance (%) of the coarse fraction i.e. larger than 63 µm (Berger et al., 1982; Gonzales et al., 2017; Suárez-Ibarra et al., 2021) and, iv) the CaCO₃ content (%) (Berger et al., 1982; Gonzales et al., 2017; Suárez-Ibarra et al., 2021). The CaCO₃ (%) was previously calculated by Rodrigues et al. (2018). The fraction >63 µm was estimated using a laser diffraction particle size analyzer Horiba Partica-LA-950X, which determined the grain size of the bulk sediment samples at the Centro de Estudo de Geologia Costeira e Oceânica (CECO) of the Universidade Federal do Rio Grande do Sul (UFRGS).

To divide the time series record into distinct cluster intervals (periods characterised by similar conditions) we carry out a clustering analysis. Posteriorly, to distinguish the dependencies of the grouping process we utilize an ordinate analysis. Both clustering and ordinate analyses are applied on the planktonic foraminifera species with relative abundances >1% (See Supplementary material ).

Additionally, since all proxies are inevitably affected by different environmental processes other than the targeted parameters, we decrease the bias by also synthesising the variation through time of i) sea surface productivity, ii) upper water column stratification, and iii) carbonate dissolution. To do so, we run principal component analyses (PCA) on the above proxies, based on the correlation matrix. The data were centralised and standardised by dividing the difference between the dataset mean and the sample value by the dataset standard deviation. The synthesised productivity, stratification and dissolution proxies were extracted from the first axes of the PCAs as PC1_P (productivity), PC1_S (stratification) and PC1_D (dissolution). Correlations between (and within) the first axes of the PCAs and other proxies were calculated using reduced major axis regressions. All PCA and correlation analyses were conducted using the software PAST (version 4.08; Hammer et al., 2001).

To test whether variations through time correspond to cyclic events, paced by orbital forcings, first we run a Multi-Taper-Method (MTM) test using the software “Acycle” (Li et al., 2019). The PC1_P curve was detrended using a locally weighted scatterplot smoothing (LOWESS) and then analysed with the MTM. We set the smoother with a Time-bandwidth product of “2”, and calculated the red noise following the Classical autoregressive (AR) model AR(1) (Husson, 2014). Second, we conducted a REDFIT spectral analysis utilising the software PAST (version 4.08; Hammer et al., 2001; Schulz and Mudelsee, 2002). The spectral analysis used the “Welch” window, a configuration of “4” for the oversampling and “2” as number of segments.

The relatively high G.bull/G.rub ratios, G. glutinata (%) and BFAR values ( Figures 2A–C ) imply a period of enhanced sea surface productivity during MIS 4 and, to a lesser extent, MIS 3. The high glacial productivity recorded by core SIS-249 agrees with previous studies for the southern Brazilian continental margin (SBCM), that also document enhanced productivity during the last glacial stage (e.g., Gu et al., 2017; Pereira et al., 2018; Portilho-Ramos et al., 2019; Frozza et al., 2020; Suárez-Ibarra et al., 2022). The increase of around 40% in G.bull/G.rub mean ratios from ~0.27 (110–83 ka) to ~0.39 (83–30 ka) suggests an increase in the upwelling of subsurface more nutrient-rich South Atlantic Central Water (Venancio et al., 2016; Lessa et al., 2017; Lessa et al., 2019; Portilho-Ramos et al., 2019).

As the upwelling of subsurface waters necessarily implies a break in the upper water column stratification, we should also expect a similar behaviour in the stratification proxies. This can be primarily confirmed by the high relative abundances of G. truncatulinoides right coiling form, which also point to a less stratified, well-mixed upper water column during MIS 4 ( Figure 2D , Lohmann and Schweitzer, 1990 ; Renaud and Schmidt, 2003; Feldmeijer et al., 2015; Billups et al., 2016; Pinho et al., 2021). Furthermore, G. truncatulinoides types II and V (right coiling forms) have also been associated with a shallow thermocline and eutrophic conditions (de Vargas et al., 2001; Ujiié and Lipps, 2009; Ujiié et al., 2010; Quillévéré et al., 2013). Finally, the less stratified upper water column during MIS 4 (and to a lesser extent MIS 3) is also supported by the lower Δδ¹⁸O _{G. inf-G.rub} values ( Figure 2E ), which implies a reduced stratification between the mixed layer and the thermocline (Santos et al., 2017b).

Considering that the utilised foraminiferal proxies also depend on other environmental parameters (for instance, changes in the relative abundances of one species can be due to variations in the abundances of other species), the use of the principal component analysis helps to synthesise the variation of the proxies through time. The PC1_P (the synthesised first principal component of the G.bull/G.rub ratios, G. glutinata (%), and BFAR proxies) from our core SIS-249 (in the SBCM) suggests a transition, from MIS 5 to MIS 4, to more eutrophic conditions in the upper water column ( Figure 3B ). This shift in nutrient availability is also evident in the clustering and PCoA analyses ( Figures S3 and S4 ), where the dominant/significant species shift from warm and oligotrophic conditions during MIS 5 (G. ruber albus, O. universa and G. menardii) to species associated with high productivity during MIS 4 (G. glutinata, G. bulloides and G. truncatulinoides dextral coiling). However, studies for the southeastern Brazilian continental margin (SEBCM) document a shift to less eutrophic conditions during the same time interval (e.g., Portilho-Ramos et al., 2015; Lessa et al., 2017; Lessa et al., 2019), associated with the eccentricity cycle. Our study suggests two different mechanisms fertilising the SBCM and SEBCM regions i) enhanced Fe-fertilisation through dust delivery and riverine input due to the strengthening of southwesterly winds during glacial stadials and ii) boosted upwelling delivering subsurface nutrient rich waters to the surface during interstadials, respectively.

Notably, our PC1_P (along with the PC1_S and PC1_D) varies alongside the obliquity cycle, as shown by the spectral analyses ( Figure 3 ). The Classic AR(1) ( Figure S5 ) and REDFIT ( Figure S6 ) give a >99% confidence for our PC1_P to be orbitally paced by the obliquity cycle (43 kyr). Our results support the hypothesis that the obliquity cycle modulates a dust delivered Fe-fertilising mechanism at the studied site (Lopes et al., 2021). Briefly, under low-obliquity values, annual average insolation decreases at the poles (Paillard, 2021), which would favour the expansion of Antarctica’s ice sheets (Doughty et al., 2021), intensifying the southwesterly winds (SWW, Toggweiler et al., 2006) and the wind-driven dust delivery. Then, during low-obliquity intervals, the north limit of the SWW belt migrates north (from 40°S to 30°S latitude, Lamy et al., 2015; Gili et al., 2017) close to our coring site. This would allow the northward/offshore transport of more nutrient-rich, fresher, and cooler waters of the RdlP closer to our study area, also fertilising the photic zone (e.g., Gu et al., 2017; Pereira et al., 2018; Portilho-Ramos et al., 2019; Bottezini et al., 2021; Bottezini et al., 2022).

In fact, a northern/offshore influence of RdlP outflows during periods of low obliquity values would be supported by the population collapses of the dinoflagellate species Operculodinium centrocarpum, associated with the Brazil Current, documented by Gu et al. (2017) (core GeoB2107-3) for the SBCM during MIS 4 and the MIS 3/2 boundary ( Figure 3D ). Also, the decreased influence of the core of the Brazil Current (which transports warm oligotrophic water) in the study area would lead to a decrease in upper water column stratification, evidenced by our PC1_S ( Figure 3B ). In addition, a northern/offshore intrusion of the nutrient-rich RdlP outflows (or freshwater bodies formed along the continental shelf) can also be inferred from the presence of freshwater diatoms during MIS 3–2 from core SIS-188 (Bottezini et al., 2021; Bottezini et al., 2022).

Moreover, strengthened SWW would also transport wind-driven particles from southern South America to the study area, increasing the Fe-fertilisation from terrestrial sources (Lopes et al., 2021). During MIS 5/4 and 3/2 boundaries, the elevated pollen concentrations (e.g., Gu et al., 2017, SBCM, core GeoB2107-3) can be interpreted as enhanced aeolian transport ( Figure 3 ), and thus, terrestrial dust fertilisation. Additionally, the presence of the Andean pollen Nothofagus, recorded by Gu et al. (2017) during the late MIS 3 and MIS 2, corroborates the idea of enhanced windiness, strengthened SWW and aeolian Fe-fertilisation from South American sources. On top of it, Lopes et al. (2021) compiled chemical studies analysing the provenance of glacial terrigenous dust in the southwest Atlantic, pointing southern South America (central-western Argentina and Patagonia) as the main source (e.g., Delmonte et al., 2010; Weber et al., 2012; Gili et al., 2017).

In summary, our record suggests two mechanisms to have fertilised the study area: i) upwelling, enhanced by the interplay of lower glacial sea level and the local bathymetry allowing the shoaling of nutrient-rich South Atlantic Central Water (also richer in silicic acid, Portilho-Ramos et al., 2019) to the sub-surface (e.g., Mahiques et al., 2007; Lessa et al., 2016) and; ii) obliquity-paced fertilisation from southern South America (wind-driven dust and the RdlP outflows) through enhanced aeolian transport, associated with expanded SWW. Both processes lead to a decrease in the upper water column stratification, first, by the ascension of cooler and less salty subsurface South Atlantic Central Water and, second, by the offshore displacement of the core of the Brazil Current. Ultimately, the correlation value between PC1_P and PC1_S (ρ<0.05, r=0.51) corroborates the synergy between the boosted biological pump and less stratified conditions.

Enhanced primary productivity (i.e., G.bull/G.rub and G. glutinata %) and decreased upper water column stratification (PC1_S) are mirrored by high accumulation rates of total organic carbon during MIS 4 (TOC*AR, Figure 3 , TOC data from Rodrigues et al., 2018, and accumulation rates from this study). In addition, the BFAR proxy, commonly employed to reconstruct past surface productivity by examining organic matter export to the seafloor (Herguera and Berger, 1991), demonstrates significant co-variation (ρ<0.05) with both G.bull/G.rub (r=0.5) and G. glutinata relative abundances (r=0.55). Although the BFAR proxy has been widely used to quantitatively estimate past productivity changes, this relationship is not always straightforward (e.g., Naidu and Malmgren, 1995; Den Dulk et al., 2000; Jorissen et al., 2007). This limitation arises from the lack of calibration across various productivity settings, among other factors (for an in-depth discussion, please refer to Jorissen et al., 2007). However, Dias et al. (2021) demonstrated that the BFAR can indicate variations in the Cabo Frio Upwelling System, located in the SEBCM, where higher BFAR values correspond to increased export of fresh marine organic carbon.

The above-mentioned coupling suggests efficient glacial carbon sequestration via biological pump in the study area. A schematic representation of coupled benthic-pelagic changes observed in mid-latitude southwest Atlantic during the 107–30 ka time interval is presented in Figure 4 .

The here evidenced increase in the export of organic matter (OM) to the seafloor from MIS 5 to MIS 4 (and to a lesser extent in MIS 3) aligns with the findings of Rodrigues et al. (2018). They observe a decrease in the infaunal/epifaunal benthic foraminifera ratio in core SIS-249, indicating a shift in the bioavailability of OM from refractory to labile forms (Gooday, 1993; Smart et al., 1994; Gooday, 2003; Garcia-Chapori et al., 2014; de Almeida et al., 2022). This change has significant implications for calcium carbonate preservation. As highlighted in our study, during phytoplankton blooms in glacial periods, higher amounts of OM (more labile, more soluble, and more easily remineralised) are exported to the seafloor, resulting in increased CO₂ release, lower pH, and the dissolution of CaCO₃ ( Figure 4 ). In addition, despite the relatively brief sinking time of planktonic foraminifera tests (Takahashi and Bé, 1984; Schiebel and Hemleben, 2017), they could experience minor dissolution during transit through the CO₃ ²⁻ under-saturated intermediate SCWs (Frenz et al., 2003). These findings shed light on the important role of sea surface productivity and water mass properties in influencing the dynamics of carbonate preservation during different climatic periods.

Despite the sparser sampling in PC1_D compared to PC1_P and PC1_S, our study found a strong correlation between PC1_P and PC1_D in core SIS-249 (ρ<0.05, r=0.80, Table 1 ). This correlation supports the idea that enhanced dissolution is triggered by glacial high productivity, as it is also documented during the last deglaciation (Suárez-Ibarra et al., 2022, Figure 3C ), even in deeper cores under the presence of modern non-corrosive Northern Component water. While inferring changes in paleoproductivity based on compositional assemblages from partially dissolved samples can be imprecise, our study addresses this issue by complementing the data with geochemical and sedimentological analyses. These additional analyses, like the δ¹⁸O signal from the species G. inflata and G. ruber (both dissolution resistant according to experimental studies from Petró et al., 2018), as well as the BFAR and the TOC*AR, are less biased by dissolution.

Due to the OM enrichment in ¹²C from photosynthesis (O'Leary, 1988; Ravelo and Hillaire-Marcel, 2007), remineralisation at the seafloor is expected to drive negative excursions in the endobenthic foraminiferal δ¹³C values, as seen in MIS 5c and 4 ( Figure 3E ). Nevertheless, although δ¹³C _Uvig is significantly anti correlated with our PC1_P (ρ<0.05), the correlation value is weak (r=-0.31) and is likely impacted by additional changes in the dissolved inorganic carbon δ¹³C of pore waters and/or the local bottom water mass geometry (Ravelo and Hillaire-Marcel, 2007; Hesse et al., 2014). Lopes et al. (2021) suggested an increased influence of more corrosive SCW at the core site during MIS 4 and 3, based on an apparent phasing between the δ¹³C _Uvig record from core SIS-249 and the dust flux record from Dome C ice core (Lambert et al., 2008), as the dust-fertilisation in the Southern Ocean would be recorded in the δ¹³C of the bottom water masses ( Figure 3E ). Although we do not exclude the possible influence of a higher proportion of SCW on carbonate preservation, our data show a relationship between the changes in the organic carbon cycle (periods of enhanced sea surface productivity associated with higher exportation of labile OM) and the deep-sea carbonate system. This highlights the need for further investigations applying a conservative or quasi-conservative water mass circulation tracer to elucidate the water masses mixing proportions at our site during the last glacial inception.

Furthermore, a third process that could have impacted the calcium carbonate preservation is the intensity of the Atlantic Meridional Overturning Circulation (AMOC). In a scenario with a sluggish AMOC, the reduced water flow at the seafloor leads to an accumulation of respired CO₂, resulting in decreased water pH and subsequent corrosion of CaCO₃ (McManus et al., 2004; Thornalley et al., 2013; Howe et al., 2016b). Despite the lower temporal resolution of our PC1_D record, it is evident that the positive excursions of ²³¹Pa/²³⁰Th ( Figure 3F ), indicative of a sluggish AMOC (McManus et al., 2004; Böhm et al., 2015), do not align with the longer-term events of intense dissolution recorded by core SIS-249, which also seem to follow the 43 kyr-cycle. Furthermore, our findings are consistent with the conclusions of Suárez-Ibarra et al. (2022), who, using a dataset with higher temporal resolution, found that surface productivity was the primary factor driving calcium carbonate dissolution close to our study site.

Additionally, we would like to highlight that these processes, i) the biological pump and ii) water mass configuration, are not mutually exclusive and both may respond to the obliquity cycle. Under low obliquity values and decreased insolation at high latitudes, the (sub)polar fronts migrate equatorward: (i) displacing the north limit of the SSW belt closer to the study site, enhancing the wind-driven dust and continental river fertilisation, and boosting the biological pump; (ii) expanding the southern ice coverage, leading to an increased brine production and enhanced SCW formation (Govin et al., 2009). This would result in the volumetric increase in the deep Atlantic of more corrosive SCWs which also act in part as a CO₂ storage reservoir.

Finally, although the enhanced remineralization causes part of the OM and biogenic carbonate to be recirculated back in the system instead of being exported to the sediments, we suggest an efficient biological pump in the mixed-layer, which removes inorganic carbon and points out the importance of the SBCM as a glacial carbon capture and sink of atmospheric CO₂ and thus part of the puzzle of glacial-interglacial CO₂ changes.