Bottom-water deoxygenation at the Peruvian margin during the last deglaciation recorded by benthic foraminifera

Deciphering the dynamics of dissolved oxygen in the mid-depth ocean during the last deglaciation is essential to understand the influence of climate change on modern oxygen minimum zones (OMZs). Many paleo-proxy records from the eastern Pacific Ocean indicate an extension of oxygen-depleted conditions during the deglaciation, but the degree of deoxygenation has not been quantified to date. The Peruvian OMZ, one of the largest OMZs in the world, is a key area to monitor such changes in near-bottom-water oxygenation in relation to changing climatic conditions. Here, we analysed the potential to use the composition of foraminiferal assemblages from the Peruvian OMZ as a quantitative redox proxy. A multiple regression analysis was applied to a joint dataset of living (rose-bengal-stained, fossilizable calcareous species) benthic foraminiferal distributions from the Peruvian continental margin. Bottom-water oxygen concentrations ([O₂]_BW) during sampling were used as the dependant variable. The correlation was significant (R²=0.82; p<0.05), indicating that the foraminiferal assemblages are rather governed by oxygen availability than by the deposition of particulate organic matter (R²=0.53; p=0.31). We applied the regression formula to three sediment cores from the northern part of the Peruvian OMZ between 3 and 8^∘ S and 997 and 1250 m water depth, thereby recording oxygenation changes at the lower boundary of the Peruvian OMZ. Each core displayed a similar trend of decreasing oxygen levels since the Last Glacial Maximum (LGM). The overall [O₂]_BW change from the LGM and the Holocene was constrained to 30 µmol kg⁻¹ at the lower boundary of the OMZ.

1 Introduction

Oxygen minimum zones (OMZs) occur where intense upwelling and high primary productivity result in elevated oxygen consumption within the water column in combination with sluggish ventilation (Wyrtki, 1962; Helly and Levin, 2004; Fuenzalida et al., 2009). In today's world oceans, the most pronounced OMZs with oxygen concentrations <20 µmol kg⁻¹ are observed offshore of northwest and southwest Africa, in the Arabian Sea and Bay of Bengal in the Indian Ocean, and along the continental margin of the eastern Pacific at low latitudes (Helly and Levin, 2004; Paulmier and Ruiz-Pino, 2009). Warmer conditions have contributed to expansion of OMZs during the last decades (Stramma et al., 2008; Schmidtko et al., 2017; Levin, 2018; Oschlies et al., 2018). Paleoceanographic reconstructions of bottom-water oxygenation during the last deglaciation are a valuable approach to understand the dynamics of OMZs during changing climatic conditions (e.g. Jaccard and Galbraith, 2012; Moffitt et al., 2015; Praetorius et al., 2015). The eastern equatorial Pacific (EEP) has been the focus of various paleoceanographic studies to unravel the dynamics of surface productivity and bottom-water oxygenation (Oberhänsli et al., 1990; Heinze and Wefer, 1992; Cannariato and Kennett, 1999; Loubere, 1999; Hendy and Pedersen, 2006; Martinez and Robinson, 2010; Moffitt et al., 2014; Scholz et al., 2014; Salvatteci et al., 2016; Tetard et al., 2017; Balestra et al., 2018; Hoogakker et al., 2018; Cardich et al., 2019). The region is characterized by a strong and shallow OMZ maintained as a result of persistent upwelling (Pennington et al., 2006). Previous studies in the region used established paleo-proxies such as sedimentary textures (laminations), productivity indicators (C_org, δ¹⁵N, biogenic opal), redox-sensitive elements (e.g. U, Mo, Cd, V, Mn) and benthic foraminiferal distributions (e.g. Jaccard et al., 2014; Moffitt et al., 2015). Most of the studies reported that cold glacial periods were generally associated with contracted OMZ, whereas warm or interglacial periods were associated with an expansion of the OMZ at intermediate depths. However, only a few studies attempted bottom-water oxygen reconstructions of the Peruvian margin, and they used records from sediment cores recovered from depths shallower than 400 m (Oberhänsli et al., 1990; Heinze and Wefer, 1992; Scholz et al., 2014; Moffitt et al., 2015; Salvatteci et al., 2016). In a review, Schönfeld et al. (2015) demonstrated that the deposition of laminated sediments indicated oxygen concentrations <7 µmol kg⁻¹ and that accumulation rates of sedimentary organic carbon could be used to quantify oxygen concentrations of >10 µmol kg⁻¹. A geochemical approach focusing on the Peruvian margin used redox-sensitive elements (Fe, Mo, U) and found a 5 to 10 µmol kg⁻¹ decrease from glacial to interglacial periods in the centre of the OMZ at depths between 100 and 500 m (Scholz et al., 2014). Because of the lack of complete records from the continental slope off Peru (Reimers and Suess, 1983; Erdem et al., 2016) and the limited applicability of some of the redox proxies (laminated sediments, Mo, U) at higher oxygen levels, paleo-oxygen reconstructions were not possible at the lower, dysoxic–oxic boundary of the Peruvian OMZ. However, benthic foraminiferal faunas were markedly structured with oxygenation at these depths (Mallon et al., 2012). Therefore, the present study aimed to reconstruct paleo-oxygen conditions since the Last Glacial Maximum (LGM) by using benthic foraminiferal records from sediment cores from the Peruvian OMZ between 900 and 1250 m water depth. We compiled all available information on living (rose-bengal-stained) benthic foraminiferal faunas from the Peruvian margin as a calibration dataset to investigate the following questions. (1) Did the Peruvian OMZ structure show differences in terms of vertical and horizontal extension since the LGM? (2) If there are such differences, can we actually quantify these changes in bottom-water oxygen concentrations ([O₂]_BW) by reverting to the oxygen demands of today's living faunas? (3) If so, how much did the [O₂]_BW levels change since the LGM?

1.1 Benthic foraminifera as an oxygen proxy

Certain benthic foraminiferal species and assemblages have been suggested as proxies for bottom-water oxygenation, in particular for low-oxic to anoxic conditions (e.g. Sen Gupta and Machain-Castillo, 1993; Kaiho, 1994; Alve and Bernhard, 1995; Bernhard et al., 1997; Baas et al., 1998; Nordberg et al., 2000; Leiter and Altenbach, 2010). Since a high flux of particulate organic matter to the sea floor prevails in OMZs, these species also flourish under elevated food availability. Applications of certain species to reconstruct ancient bottom-water oxygen concentrations were often undermined by the TROX model (Barmawidjaja et al., 1992; Jorissen et al., 1995). This conceptual model explains the microhabitat structure of benthic foraminifera in the sediments as driven by the availability of both the organic matter and dissolved oxygen (Van der Zwaan, 1999). A growing number of publications reporting the living (rose-bengal-stained) benthic foraminiferal distributions and their ambient environmental conditions (Phleger and Soutar, 1973; Mackensen and Douglas, 1989; Sen Gupta and Machain-Castillo, 1993; Bernhard et al., 1997; den Dulk et al., 1998; Jannink et al., 1998; Levin et al., 2002; Schumacher et al., 2007; Cardich et al., 2012, 2015; Mallon et al., 2012; Caulle et al., 2014) showed some features in common. First, benthic foraminiferal faunas generally show a low diversity and high population density in oxygen-depleted environments; (2) most but not all living specimens in OMZs were found dwelling in the first 1 or 2 cm of the surface sediments here, even at moderate flux rates of particulate organic matter; (3) species with a thin, porous test wall (e.g. Bolivinids) always outnumber the agglutinated and porcelaneous species at low oxygen levels. Pore densities of the tests have been recognized as indicators of bottom- and pore-water redox conditions (Kaiho, 1994; Glock et al., 2011; Kuhnt et al., 2013; Rathburn et al., 2018). A comparison of different OMZ settings showed that benthic foraminiferal assemblages and distributions could be used to identify spatial changes of the OMZ provided certain threshold values and ranges are considered. However, the terminology, threshold values and units used in previous benthic foraminiferal studies were markedly different (Table 1). Here, we consider the following classification: microxic conditions <5 µmol kg⁻¹, dysoxic conditions 5–45 µmol kg⁻¹ and oxic conditions >45 µmol kg⁻¹. Overall, an extreme low-oxygen, even anoxia-tolerant association was found within the OMZ core, a transitional species group was recorded around the lower boundary of the OMZ (>20 µmol kg⁻¹), and a cosmopolitan and much more diverse fauna was observed outside the OMZ (e.g. Table 2; Schumacher et al., 2007; Mallon et al., 2012; Caulle et al., 2014).

Table 1 Classification of different environments and thresholds of bottom-water dissolved oxygen ([O₂]_BW) for benthic foraminifera and benthic biota in general.

Table 2 Common species considered in this study and reported from different oxygen-depleted environments.

1.2 Regional setting

The Peruvian OMZ is one of the most pronounced OMZs in the world (Fig. 1; Paulmier and Ruiz-Pino, 2009), covering the Peruvian continental shelf and upper slope, with its thickest part between 5 and 15^∘ S and 50 and 750 m water depth (Fig. 2; Fuenzalida et al., 2009). The intensity of the OMZ is dependent on the low ventilation of advected intermediate waters, diapycnal mixing and the extremely high primary productivity in the surface waters (Karstensen et al., 2008; Brandt et al., 2015). The productivity is maintained by the wind-driven upwelling of cold, nutrient-rich and oxygen-poor waters from intermediate depths (e.g. Pennington et al., 2006). The main source of these upwelled waters is the Peru–Chile Undercurrent (PCUC). It originates around 3–5^∘ S and flows southward between 50 and 300 m water depth (Montes et al., 2010; Chaigneau et al., 2013). The PCUC is fed by the equatorial undercurrent (EUC) and southern subsurface countercurrents (SSCCs; Montes et al., 2010). Below the PCUC, the northward-flowing Chile–Peru Deep Coastal Current (CPDCC) carries cold Antarctic Intermediate Waters as a thin layer (AAIW; Chaigneau et al., 2013).

Figure 1 Location map of the study area (square in a) in the eastern equatorial Pacific and structure of the OMZ at 400 m water depths. The purple area and contour line indicate dissolved oxygen values <0.5 mL L⁻¹ (ca. <20 µmol kg⁻¹); World Ocean Database 2013 (Boyer et al., 2013). (b) Detail map of the study area showing the locations of the surface samples (x) and the sediment cores (red circles). See Fig. 2 for the sample locations in relation to water depth and the OMZ intensity.

2 Material and methods

2.1 Sediment cores

Three sediment cores were considered in this study (M77/2-50-4, 52-2 and 59-1). They were collected during expedition M77 leg 2 with R/V Meteor in 2008 from the continental slope between 3 and 9^∘ S and water depths of 997 m and 1250 m around the lower boundary of today's OMZ (Figs. 1 and 2, Table 3). The age models of all cores were based on radiocarbon dating from the planktonic foraminiferal species Neogloboquadrina dutertrei which is supplemented with benthic oxygen isotope curves that were correlated to stacked standard records (Fig. 3; Stern and Lisiecki, 2014) and Antarctic ice cores (Members et al., 2006). The radiocarbon datings were performed at the Leibniz Laboratory for Radiometric Dating and Stable Isotope Research, University of Kiel (CAU) and by Beta Analytic Inc., Florida. For core M77/2-59-1, the data were reported by Mollier-Vogel et al. (2013) and for the other cores by Erdem et al. (2016). The radiocarbon ages were later calibrated applying the Marine13 marine calibration set (Reimer et al., 2013). Reservoir age corrections were done according to the marine database (http://calib.org/marine/, last access: 12 June 2020) ranging from 89 to 338 years (Erdem et al., 2016). The closest ¹⁴C reservoir age records are from the equatorial Pacific (Zhao and Keigwin, 2018) and off Chile (Siani et al., 2013; Sarnthein et al., 2019), and they significantly differ from each other. The reservoir ages at our three study sites likely varied over the last deglaciation, and yet due to their regional proximity, they are not expected to be significantly different from each other. All ages are expressed in thousands of years (ka) before 1950 CE (abbreviated as cal ka). We focused on the following time intervals with 300- to 500-year resolution at each core: the late Holocene (LH; 3–5 cal ka), the early Holocene (EH; 8–10 cal ka), the Bølling Allerød/Antarctic Cold Reversal (BA/ACR; 13–14.5 cal ka), the Heinrich Stadial-1 (HS1; 15–17.5 cal ka), and the Last Glacial Maximum (LGM; 20–22 cal ka). For benthic foraminiferal analyses, 10–20 cm³ sediment samples were wet sieved on a 63 µm screen immediately after they were taken, and the residues were dried at 40 ^∘C. They were later split with an Otto microsplitter when needed, in order to attain similar total numbers of specimens, around 300 per sample (Murray, 2006). The foraminifera were dry picked, collected in Plummer cell slides, sorted by species, fixed with glue and counted. Benthic foraminiferal assemblage compositions and taxonomic references were previously reported elsewhere (Erdem and Schönfeld, 2017). Bulk sedimentary TOC (%) and δ¹⁵N data of core M77/2-59-1 were taken from Mollier-Vogel et al. (2019), whereas TOC and δ¹⁵N of core M77/2-52-2 were taken from Doering et al. (2016) and Glock et al. (2018), respectively (Table 3).

Table 3 Metadata of sediment cores used in this study.

NA: not available.

Figure 2 (a) Depth vs. latitude profile of the dissolved oxygen concentrations measured during the M77 legs 1 & 2 expeditions (red line in Fig. 1, b show the position of this profile; CTD data compilation after Schönfeld et al., 2015), together with the location of the surface samples (x) and the sediment cores (circles and core names). (b) Particulate organic carbon rain rates (RRPOC) taken from Dale et al. (2015), otherwise calculated following Martin et al. (1987). Rain rates are grouped indicating different values (<5, 5–10, 10–20 and >20 mmol m⁻² d); see text for details. Contour lines are the same as in (a).

2.2 Surface samples and living benthic foraminifera

Information on the living (rose-bengal-stained) benthic foraminifera was compiled from four independent datasets. These studies did not report the dead (not stained) benthic foraminiferal assemblages. This is mainly due to the fact that empty (dead) foraminiferal tests are very rare and outnumbered by living specimens in surface samples from the upper Peruvian OMZ. Since the majority of the specimens in surface sediments were stained, it is not likely that there would be a substantial difference in the statistical analysis if the dead assemblage were used. The compiled dataset comprises 53 samples from the Peruvian continental shelf and slope from water depths of 48 to 2092 m between 1^∘45^′ and 17^∘28^′ S (Fig. 2, Table 4). Four of the samples from a transect around 12^∘30^′ S were collected in December–January 1998 during the Panorama Expedition, leg 3a, with R/V Melville. Eight surface sediment samples and environmental parameters were collected from the continental shelf and uppermost slope around 12 and 14^∘ S during different monitoring cruises with R/V SNP 2 and José Olaya Balandra in August and April 2009, 2010 and 2011. For the present study, averaged bottom-water oxygen values are used for these stations (for details see Cardich et al., 2015). The largest dataset was gathered from 33 surface samples collected in October to December in 2008 during R/V Meteor expedition M77 legs 1 and 2 (Mallon, 2012; Mallon et al., 2012). The eight most recent samples and additional data were collected in May 2017 during R/V Meteor expedition M137 (Sommer, 2017). The surface sediment samples in all studies comprise the topmost 10 or 30 mm. The faunal census is based on the >63 µm size fraction. In the case of fractionated subsamples, we combined the values of different grain size fractions considering the volumes and splits reported for each subsample.

Table 4 Metadata of surface samples from the living benthic foraminifera dataset. Particulate organic carbon rain rates (RRPOC) of stations indicated in bold were taken from Dale et al. (2015); others were calculated (see text for details). Bottom-water oxygen concentrations were taken from reference publications.

2.2.1 Data reduction: consideration of taphonomy

The inventory of living benthic foraminiferal species was compared with that from the sediment cores after compilation of the joint dataset. As expected, only single specimens of three agglutinated species were found in some samples from cores 50-4 and 52-2. Agglutinated species have a lower preservation potential after burial in the sediment because their organic cement is decomposed during early diagenesis or changing redox conditions (Schröder, 1988; Mackensen et al., 1990). Consequently, only species with calcareous tests were considered for further analyses. Even though agglutinated species were not used for our downcore application, they are well known for their low tolerance to oxygen minimum conditions (Bernhard and Bowser, 1999; Gooday and Rathburn, 1999; Gooday et al., 2000; Levin et al., 2002; Mallon, 2012). We therefore considered their abundances separately, and the proportions of agglutinated species in living faunas were used for a comparison of the results.

Another exception was made for Hoeglundina elegans. The aragonitic test of this species has a low preservation potential (Gonzales et al., 2017). Hoeglundina elegans was observed in four surface samples with more than 5 % abundance, but it was almost absent from the foraminiferal assemblages of sediment cores. Therefore, we excluded H. elegans from the census data as well. The faunal data from surface sediment and core samples were reduced accordingly.

Figure 3 Benthic δ¹⁸O isotope curves of sediment cores with age control points. Black lines are reference records: EPICA ice core (Members, 2006) and intermediate water Pacific stack (IP stack; Stern and Lisiecki, 2014). Symbols indicate the tie points used for age models. Each core is shown with different colours throughout this article.

Considering the water depths at which the sediment cores were taken (>900 m) and the depth range of our surface samples (Fig. 2), we further reduced the dataset by taking only surface samples collected from water depths >300 m into account. Species showing at least three occurrences with percentages of more than 5 % from these samples were listed and considered for further analyses. The final reduced dataset included in total 16 species from 35 samples.

2.2.2 Environmental data

Bottom-water oxygen concentrations ([O₂]_BW) were measured at the time of sampling and reported to vary between 0.0 and 100.4 µmol kg⁻¹ (Table 4). Dissolved oxygen concentrations at the R/V Melville stations around 12^∘30^′ S were previously reported (Levin et al., 2002). For the eight R/V SNP 2 and José Olaya Balandra stations, [O₂]_BW data were measured on each cruise (Cardich et al., 2015). We thus considered averaged oxygen concentrations for the R/V SNP 2 and José Olaya Balandra stations. The [O₂]_BW of the R/V Meteor M77 stations was taken from Mallon et al. (2012) and from a synoptic compilation of all conductivity–temperature–depth (CTD) hydrocasts from this area (Schönfeld et al., 2015). The [O₂]_BW of the R/V Meteor M137 stations was extracted from the expedition dataset (Gerd Krahmann, personal communication, 2018). Bottom-water oxygen concentrations around the lower boundary in the region are relatively stable throughout the year, whereas the upper boundary is quite dynamic (Paulmier and Rioz-Pino, 2009). We compared the [O₂]_BW from stations close to each other, and those sampling campaigns took place in different years and seasons. The only difference was observed at two 820 m stations (January 1998 vs. May 2017; ca. 10 µmol kg⁻¹); the other stations did not show significant contrast between different sampling periods (Table 4). Therefore, we consider the reduced living assemblage data to be least affected by seasonality of [O₂]_BW.

Rain rates of particulate organic carbon (RRPOC; mmol m⁻² d) for six of the stations were taken directly from Dale et al. (2015; Table 4). The rain rates for the other stations were estimated by using equations provided for water depths between 100 and 1000 m (Martin et al., 1987; Dale et al., 2015). Different primary productivity values were used for the RRPOC calculations. For the region around 11–12^∘ S, values reported in the Supplement of Dale et al. (2015) were used since the measurements were done close to sampling during the M77 expedition (Mallon et al., 2012). For the region around and south of 15^∘ S, values reported by Martin et al. (1987) were considered. For the northern part of the study area, estimates from Pennington et al. (2006) were used.

2.3 Statistical analyses

We standardized the reduced data matrix by calculating the proportions of the involved species referring to the total calcareous (calcitic) species of each sample as described above. Relative abundances (percentages) were preferred instead of absolute abundances (individuals per cubic centimetre) since this information would create an inconsistency when applied to surface sediment samples and to fossil sediments in the same manner. Diversity and dominance values were calculated for the reduced calcareous (calcitic) assemblages applying centred bootstrapping, together with Q-mode hierarchical cluster analysis and canonical correspondence analysis (CCA), including comparison with the environmental variables [O₂]_BW and RRPOC.

All statistical and diversity analyses were performed with the PAleontological STatistics (PAST) software, version 3.11 (Hammer et al., 2001). Q-mode hierarchical cluster analysis was applied using the unweighted pair group method (UPGMA) based on a Bray–Curtis similarity matrix. Canonical correspondence analysis (CCA) was performed on the same datasets to see the relation between the species, stations and environmental variables (Fig. S3 in the Supplement). Additionally, multiple regression analysis was applied, and the significance was assessed in order to evaluate the reliability for downcore applications. Coefficients and intercept value from the multiple regression analysis were later used as data entry for a polynomial transfer function (Supplement) to calculate past bottom-water oxygen concentrations ([O₂]_BW) from foraminiferal assemblages of sediment core samples.

3 Results

3.1 Living benthic foraminiferal distributions

The entire dataset of living calcareous benthic foraminifera from the Peruvian margin comprised 53 surface sediment samples and 127 different calcareous (calcitic) species. When the oxic (>45 µmol kg⁻¹), dysoxic (5–45 µmol kg⁻¹) and microxic (<5 µmol kg⁻¹) classification (Table 1) was applied to this dataset, 27 samples were classified as microxic, 20 samples were classified as dysoxic and the remaining 6 samples were classified as oxic (Fig. 4). The diversity of calcareous species and the relative abundance of the total of agglutinated species increased with bottom-water oxygen. In particular, a marked increase in the proportion of agglutinated species was observed at stations with [O₂]_BW >15 µmol kg⁻¹, and a further increase was recorded around 40 µmol kg⁻¹. The proportion of the agglutinated taxa was more than 50 % at sampling stations under oxic conditions. The assemblages were, however, less diverse compared to dysoxic samples, as depicted by low Fisher alpha indices and the dominance of single species. This pattern implies that these six oxic stations did not represent the total calcareous taxa very well; neither do the stations in the reduced dataset of 16 species and 35 samples that we statistically analysed. The microxic samples showed less diverse assemblages with higher dominance, whereas the dysoxic samples showed higher diversities and a lower dominance. The Fisher alpha diversity indices were lower and dominance was higher at stations under dysoxic conditions where the living faunas were dominated by Bolivinids (Fig. 4).

Figure 4 Relative abundances of Bolivinids and agglutinated species at each station in relation with the bottom-water oxygen measured during sampling. Oxic (>45 µmol kg⁻¹), dysoxic (5–45 µmol kg⁻¹) and microxic (<5 µmol kg⁻¹) classification was applied to these samples, and the colours green, orange and black were used respectively as indicators of the thresholds (see also Table 1). Dominance and Fisher alpha diversity indices were also calculated for the same samples.

In the centre of the OMZ, at microxic stations, the most abundant species was Bolivina seminuda (Fig. S1), followed by B. costata, B. interjuncta, B. spissa and Bolivinita minuta with increasing oxygen concentrations. Uvigerina peregrina became the most abundant species outside the OMZ core, showing similar trends as agglutinated species and other calcareous species which increased markedly as well. Pyrgo murrhyna and Melonis barleeanum were observed only at oxic stations. Q-mode hierarchical cluster analysis on the entire census data (53 samples and 127 calcareous species) indicated three clusters (Fig. S2). The census dataset was later reduced to the most abundant 16 species from stations deeper than 300 m. All of the species from this subset were grouped in two clusters (Fig. S2a, b). Application of CCA to these abundant species showed similar results (Figs. S3 and S4). Most of the samples and species grouped together as indicated by their positive loading of axis 2. Overall, lower [O₂]_BW and higher RRPOC did not show an eye-catching relation as expected. When multiple regression results were used to compare these observations, the fitted [O₂]_BW showed a significant correlation with measured values at the sampling stations whereas the RRPOC estimates did not: R²=0.82, p<0.05 and R²=0.53, p=0.31, respectively (Fig. 5 and Table 5).

Table 5 List of the living benthic foraminifera species which were considered in the transfer function, regression coefficients and 1σ errors that were calculated with multiple regression analyses.

Figure 5 Graphs showing the regression results applied to the living benthic foraminifera dataset for bottom-water oxygen ([O₂]_BW) and POC rain rates (RRPOC).

3.2 Application of [O₂]_BW estimates to sediment cores

Erosion, reworking and high energetic bottom conditions prevail at the continental slope of the Peruvian margin. The Holocene from core 50-4 was missing (Erdem et al., 2016), and due to the high sedimentation rates, core 59-1 covers only the late Holocene (LH), early Holocene (EH), Bølling Allerød/Antarctic Cold Reversal (BA/ACR), and Heinrich Stadial-1 (HS1). Consequently, it was only possible to compare these time intervals in all of the sediment cores. In the following, we describe the results of [O₂]_BW quantification for each core separately, from south to north. We abstained from applying the same approach to reconstruct past RRPOC for the downcore records because the regression analyses did not show a significant correlation. Additionally, these cores were recovered around or deeper than 1000 m, which was the maximum depth for any reliable RRPOC calculations.

Core M77/2-50-4 was collected from 1013 m water depth, and the [O₂]_BW was 52 µmol kg⁻¹ during the time of sampling. In total, 20 samples were analysed covering the time intervals BA/ACR, HS1 and LGM; the Holocene was missing in this core. Among the 138 benthic foraminiferal species identified, only one was agglutinated (Dorothia goesi) and occurred as a single specimen in samples from the early HS1 and LGM. The estimated paleo-oxygen concentrations ranged from 35 to 43 µmol kg⁻¹ during the LGM, varied between 21 and 40 µmol kg⁻¹ during the HS1, and varied between 9 and 15 µmol kg⁻¹ during the BA/ACR. Deviations ranged from 9 to 17 µmol kg⁻¹. Core M77/2-52-2 was collected from 1249 m water depth. The [O₂]_BW was 74 µmol kg⁻¹ at the core location during sampling. In total, 27 samples were analysed, and 170 species were identified of which three were agglutinated. The estimated [O₂]_BW indicated stable conditions during the LGM with values ranging from 52 to 61 µmol kg⁻¹, which was followed by a decrease from 56 to 46 µmol kg⁻¹ from mid-HS1 to mid-BA/ACR and a much more distinct decrease at the end of BA/ACR from 53 to 10 µmol kg⁻¹. After the deglaciation, the fluctuations were weaker, indicating rather stable conditions with values of 13 to 30 µmol kg⁻¹ during the EH and 23 to 32 µmol kg⁻¹ during the LH. The standard deviations (1SD) ranged from 8 to 20 µmol kg⁻¹. Core M77/2-59-1 was recovered from the northernmost part of the study area, from 997 m water depths, and the [O₂]_BW was 54 µmol kg⁻¹ during sampling. The core location has been under the influence of strong riverine input. The sedimentation rates throughout the core were 50 and 170 cm ka⁻¹, rather high compared to the other cores (Mollier-Vogel et al., 2013). Because of these high sedimentation rates, the LGM was not retrieved by this core. In total, 20 samples were analysed and 161 species were identified. Agglutinated species were scarce but present. They were more frequent in samples older than 9 cal ka. The estimated [O₂]_BW ranged from 34 to 60 µmol kg⁻¹ during HS1, from 63 to 18 µmol kg⁻¹ during the BA/ACR, and varied between 11 and 30 µmol kg⁻¹ during the EH and between 3 and 25 µmol kg⁻¹ during the LH. One sample showed a negative value, obviously an artefact of the method. Overall standard deviations were calculated as ranging from 9 to 28 µmol kg⁻¹. When the average values were considered in each time interval, the decrease from HS1 to the end of BA/ACR was 16 µmol kg⁻¹ in core 50-4, 20 µmol kg⁻¹ at core 52-2 and only 6 µmol kg⁻¹ in core 59-1. The difference between the BA/ACR and the early Holocene [O₂]_BW was around 11 µmol kg⁻¹ at core 52-2 and around 24 µmol kg⁻¹ at core 59-1.

4 Discussion

4.1 Living benthic foraminifera in relation to OMZ settings

Similar to previous observations from other modern OMZs, benthic foraminifera living in the Peruvian OMZ showed high population densities and low diversity in the centre. Certain species, predominantly Bolivinids, were observed to be the most abundant species here. They are known for their high tolerance to suboxic and even anoxic conditions (Table 2; Mullins et al., 1985; Schumacher et al., 2007; Piña-Ochoa et al., 2010; Glock et al., 2011; Mallon et al., 2012; Cardich et al., 2015). In the vicinity of the OMZ core, around the boundary, a more diverse assemblage (cluster A; Fig. S2) was observed. Considering the environmental information gathered on some of the species of this group, such as Bolivina spissa, Bolivinita minuta, Cassidulina delicata, Epistominella pacifica and Uvigerina peregrina (Table 2), the assemblage seemingly represented a transitional community which reacted to a broader range of environmental variables including increasing water depths, oxygen content and organic carbon content. The living taxa from stations below the lower boundary of the OMZ indicated a much more diverse assemblage (cluster B; Fig. S2). This assemblage involved mostly Miliolids which are known for their intolerance to low oxygen (e.g. Caulle et al., 2014). Agglutinated species dominated the whole assemblage from these samples, which again mirrors rather oxic conditions. The data are in good agreement with previous observations on the low tolerance of agglutinated species to oxygen-depleted conditions (<0.2 mL L${}^{-\mathrm{1}}=\sim \mathrm{9}$ µmol kg⁻¹; Gooday et al., 2000). In the CCA on the living fauna, this Miliolid-associated assemblage (cluster B), together with the agglutinated species, showed a positive relation with [O₂]_BW (Fig. S3). Similar trends of such different assemblages in relation with changing oxygen and organic matter were also observed in the Arabian Sea (Jannink et al., 1998; Caulle et al., 2014). Even though there are different species involved in the assemblages representing the OMZ core (microxic), around the boundary of the OMZ (dysoxic) and outside the OMZ (oxic), the transitional appearance of these assemblages from low diversity to high density to more diverse and cosmopolitan assemblages indicates strong similarities.

Comparison of living taxa in different OMZ settings revealed that each OMZ has its own genuine assemblages. The most abundant species were not observed at similar abundances in other OMZs, indicating that they are specifically adapted to the conditions in these regions. For example, Bolivina dilatata is dominant in the Arabian OMZ (Jannink et al., 1998), and Bolivina costata is frequent in the core of the Peruvian OMZ (this study; Cardich et al., 2012; Mallon et al., 2012). While Bolivina dilatata is widely distributed in the Atlantic Ocean too, including the Mediterranean, Bolivina costata is seemingly endemic to the western South American margin. This “trapped” occurrence might be due to the structure and shape of the OMZ, prevailing for more than 0.5 million of years, and the strong adaptation of these assemblages to conditions during a long time (e.g. Heinze and Wefer, 1992). Similar assumptions were made for the OMZ core assemblages in the northern Arabian Sea (Jannink et al., 1998). Therefore, expecting to find the same specific species within the same oxygen range in different OMZs is potentially misleading. Although there are not many frequent species observed in different OMZs, Bolivina seminuda and Bulimina exilis are within the few common species (Table 2; Bernhard et al., 1997; den Dulk et al., 1998; Cardich et al., 2015). Elevated proportions of these two species could be used as an extremely low-bottom-water-oxygen indicator in downcore records (e.g. McKay et al., 2015; Praetorius et al., 2015; Tetard et al., 2017). Additionally, total numbers of the Miliolids could be oxic condition indicators as previously suggested (den Dulk et al., 2000), but this approach would not produce sensible results outside OMZ settings. Accordingly, comparisons of relative abundances of these low-oxygen-tolerant and low-oxygen-intolerant species in the fossil record might be used in determining dysoxic–oxic transitions (e.g. Kaiho, 1994; Cannariato et al., 1999; Schmiedl et al., 2003; Tetard et al., 2017; Balestra et al., 2018).

4.2 Peruvian margin oxygen history since the LGM

The records revealed a distinct overall decrease in [O₂]_BW during deglaciation around the lower OMZ boundary. The spatial delineation of oxygenation changes between the LGM and the Holocene is limited since our records do not fully cover all the considered time intervals. When the records were stacked, the estimates showed a decreasing trend starting from the LGM, a distinct drop during the deglaciation with fluctuations from the HS1 to the BA/ACR and a slight increase followed by relatively stable concentrations during the Holocene. Nonetheless, Holocene [O₂]_BW values were still lower than those of the LGM. This trend is consistent with other results reported in reviews of bottom-water deoxygenation during the deglaciation in the eastern Pacific Ocean (Jaccard and Galbraith, 2012; Moffitt et al., 2015). Both reviews consider different proxies (e.g. lamination, δ¹⁵N, redox-sensitive elements) which are known to indicate oxygen-depleted conditions as recorded by sediment cores from above 500 m depth off Peru. When the records are evaluated individually, relatively stable LGM [O₂]_BW values at both available records were followed by a continuous decrease during the deglaciation in the southernmost core. The [O₂]_BW from the northern cores suggested that bottom waters experienced increased oxygen levels during the early BA/ACR. This is potentially due to weakening of surface productivity or, considering the positions of both cores on the continental slope (900–1300 m), due to intermediate water mass advection of oxygen-rich waters. Water mass contributions might have changed during the deglaciation at these cores sites as previously reported (e.g. Antarctic Intermediate Water (AAIW) from the south and North Pacific Intermediate Water (NPIW) from the north; Bova et al., 2018). Overall, the highest decrease was observed in the deepest core, implying that the absolute changes at the lower boundary were larger than in the shallower depths. Another quantification approach compared Fe concentrations and Mo∕U ratios in core M77/2-24-5 from the upper slope off Peru at 11^∘ S and 210 m water depths (Scholz et al., 2014). They found a drop of 5 to 10 µmol L⁻¹ during the deglaciation, suggesting that when the OMZ intensifies (or diminishes) the change is profound around its borders and the conditions are rather stable close to its centre. More recently, a modelling approach quantified oxygen concentrations during the deglaciation and reported similar estimates with relatively stable conditions in shallower depths and a larger decrease in intermediate and deep waters (their Fig. 6; Yamamoto et al., 2019).

Figure 6 Comparison of the estimated bottom-water oxygen concentrations [O₂]_BW applied to three sediment cores from south to north. Stars indicate the modern [O₂]_BW measured during sampling in 2008. Error bars are calculated using 1σ errors using values shown in Table 5 and using equations described in the Supplement.

The downcore distribution of benthic foraminiferal species that were not included in the quantification approach provided further corroborating evidence. For instance, the disappearance of Prygo murrhyna and agglutinated species after HS1 at core 50-4 and the increasing abundances of B. costata during the same time interval suggested that bottom-water conditions successively became dysoxic (Erdem and Schönfeld, 2017). Similarly, the presence of P. murrhyna and agglutinated species throughout core 52-2 suggested that the prevailing oxygen levels have been moderate during the considered time intervals. Due to hiatus and sampling resolution, our assessments for the Holocene are limited. In addition, large standard deviations observed in the northernmost core 59-1 raised questions about the applicability of the method for the Holocene at this core (Fig. 6). Nevertheless, both Holocene records of cores 59-1 and 52-2 suggested a trend of increasing bottom-water oxygenation that is in accordance with recently published results from the region (Salvatteci et al., 2016, 2018; Mollier-Vogel et al., 2019). Moreover, it is possible that bottom waters became more oxic after the late Holocene as reported for the shelf during the last 100–150 years (Cardich et al., 2019). However, we cannot comment further for the rest of the Holocene [O₂]_BW trend on the basis of currently available information.

Finally, the estimated values together with standard deviations indicated negative values, which might be considered a challenge for the applicability of this approach. Once we looked closer, only one data point at core M77/2-59-1 indicated a negative value (Table S1 in the Supplement), which is not physically possible. The negative values shown on the scale bars (in Fig. 6) are inferred by the large standard deviations from a few samples. This is potentially an artefact of the transfer function where the estimates are biased towards lower values, because the majority of reference surface samples are from depths shallower than the sediment cores investigated, hence at lower oxygen concentrations. And yet, our estimations are within a realistic range of positive [O₂]_BW within the statistical uncertainty. Moreover, we are confident in the [O₂]_BW differences in each time interval considered, even though the absolute estimates for each sample might be biased because of the dominance of the low-oxygen samples in the reference dataset. As with any new approach, precision will improve as additional samples and more species information become available.

4.3 Comparison with other proxies and records from the region

Total organic carbon content (TOC (wt %)) from cores 52-2 and 59-1 showed an increasing trend with the onset of the deglaciation followed by relatively similar trends during the Holocene (Fig. 7; Doering et al., 2016; Mollier-Vogel et al., 2019). Higher TOC values during the Holocene suggested higher preservation and thus enhanced productivity and relatively low oxygen values prevailing at core locations since the EH in comparison to the LGM. Meanwhile, sedimentary δ¹⁵N values from both cores indicated a higher surface productivity and/or water column denitrification and hence a stronger OMZ during the deglaciation compared to the LGM and the Holocene (Glock et al., 2018; Mollier-Vogel et al., 2019). The opposite δ¹⁵N_sed trend observed at both cores during the BA/ACR is intriguing. In the case of more oxygenated bottom waters, as depicted by our [O₂]_BW estimates from cores 59-1 and 52-2, δ¹⁵N_sed is expected to be lower as observed at core 52-2. This opposing trend suggests that different regional dynamics had an interfering impact on N dynamics at core 59-1 during the early BA/ACR. On the other hand, TOC values from both records are in agreement with bottom-water oxygenation estimates from our [O₂]_BW approach during this period.

Figure 7 Comparison of [O₂]_BW estimated concentrations with bulk sedimentary δ¹⁵N and TOC (wt %) from the same cores (data are available for cores M77/2-59-1 and M77/2-52-2).

Later in the record both cores indicated a distinct drop in the bottom-water oxygen concentrations that were also mirrored by the increasing denitrification during the BA/ACR as displayed by high δ¹⁵N_sed values. Other redox and productivity proxies from the region revealed similar trends. Laminated sediments indicating [O₂]_BW of <7 µmol kg⁻¹ (e.g. Schönfeld et al., 2015) were observed to be widest in extent during HS1, in particular on the continental shelf and upper slope between 10 and 18^∘ S (Erdem et al., 2016). Stacked records from the Peruvian margin at 14^∘ S revealed a correlation between enhanced productivity and low bottom-water oxygen as indicated by increasing Mo∕U values within the laminated sediments during the HS1 period (Salvatteci et al., 2016). Enhanced surface productivity, increasing denitrification and bottom-water deoxygenation during the early deglaciation were observed in various records from the eastern equatorial Pacific (Pedersen, 1983; Pedersen et al., 1988; Hendy and Pedersen, 2006; Martinez and Robinson, 2010; Bova et al., 2018). The distinct increase in relative abundances of the phytodetritus bloom-feeding Epistominella exigua in cores 52-2 and 59-1 during the deglaciation supported these observations, indicating an enhanced organic matter flux to the sea floor (Erdem and Schönfeld, 2017). Nevertheless, the highest percentages of E. exigua at core 52-2 were recorded when the bottom-water oxygenation indicated more or less stable conditions. This decoupling of productivity and deoxygenation suggested that bottom-water deoxygenation did not always co-occur with enhanced surface productivity at the core location. The discrepancy might be a corroborating evidence that parameters other than primary productivity were influencing the oxygen dynamics in the Peruvian OMZ as they do today, for instance ventilation of intermediate waters or changes in the hydrodynamics (Karstensen et al., 2008; Brandt et al., 2015). However, it should be kept in mind the Peruvian margin potentially represents a much more complex structure in terms of productivity as it is mirrored in biogenic opal and δ³⁰Si records (Doering et al., 2016) as well as the impact of different intermediate water masses and stratification (Bova et al., 2015, 2018). The same caution should be applied to our E. exigua records. The species was observed in the northern cores but not in core 50-4 situated much closer to the main upwelling area. This might be explained by the upper food flux tolerance limit of ca. 200 g C m⁻² yr⁻¹ for this species (Altenbach et al., 1999: their Fig. 6). Differences observed in bottom-water oxygenation in our records might also be an indication of changing dynamics in the mixed layer (e.g. primary productivity and upwelling) instead of changes in intermediate water masses. However, further and broader investigation focusing on other records and proxies is needed to confirm this. Even though our estimated [O₂]_BW values were lower than anticipated, the quantification approach is consistent in terms of absolute changes and coherent with other proxies. Our data and comparisons with other proxies indicated an expansion of the northern part of the Peruvian OMZ in terms of thickness and wideness during the last deglaciation with a decrease in [O₂]_BW of 30 µmol kg⁻¹ at its lower boundary.

5 Conclusions

The use of benthic foraminiferal assemblages as a bottom-water oxygenation proxy has been under debate since the oxygen deficiency indicator species can be found in many other environments as well. When certain thresholds are applied, for instance microxic (<5 µmol kg⁻¹), dysoxic (5–45 µmol kg⁻¹) and oxic (>45 µmol kg⁻¹), benthic foraminiferal assemblages were observed, indicating similar transitional trends in different oxygen minimum zone settings world-wide, even though they are composed of different species. The present study reports an extensive dataset based on four independent studies of living (rose-bengal-stained) benthic foraminiferal distributions from the continental shelf and slope off Peru. The faunal distribution data were compared with bottom-water oxygen concentrations ([O₂]_BW) measured during the sampling periods. Certain species and assemblages showed a much better correlation with [O₂]_BW than with rain rates of particulate organic carbon (RRPOC). Application of a multiple regression analysis with [O₂]_BW as the dependant variable indicated that the foraminiferal assemblages along the Peruvian margin are rather governed by oxygen availability than by the deposition of particulate organic matter. The correlation with [O₂]_BW was significant (R²=0.82; p<0.05); therefore we applied the transfer function to three sediment cores taken from the lower boundary of the Peruvian Oxygen Minimum Zone (OMZ) in order to quantify the past [O₂]_BW. The data revealed a drop in [O₂]_BW of 30 µmol kg⁻¹ at the lower boundary of the OMZ during the last deglaciation. The overall bottom-water oxygen trend revealed a distinct drop during the deglaciation, with a slight increase around 15 cal ka observed in the northern cores. The largest decrease in [O₂]_BW was observed at the later part of the deglaciation, which was followed by a slight increase as recorded in the northern cores during the Holocene. This general trend is in line with previous paleo-oxygenation proxy records at intermediate depths from the eastern Pacific Ocean, supporting the viability of the benthic foraminiferal approach.