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 15 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 20 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 2 ] BW ) during sampling were used as dependant variable. The correlation was significant (R 2 = 0.82; p<0.05) indicating that the foraminiferal assemblages are rather governed by oxygen availability than by the deposition of particulate organic indicator found in many other When certain thresholds are applied, for instance microxic (<5 μmol/kg), dysoxic (5 - 45 μmol/kg), oxic (>45 μmol/kg), benthic foraminiferal assemblages were observed indicating similar transitional trends in different oxygen minimum zone 5 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 2 ] BW ) measured during the sampling periods. Certain species and assemblages showed a much better correlation with [O 2 ] BW than with rain rates of particulate organic carbon (RRPOC). Application of a 10 multiple regression analysis with [O 2 ] BW as 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 2 ] BW was significant (R 2 = 0.82; p<0.05), therefore we applied the transfer function to four sediment cores taken from the lower boundary of the Peruvian Oxygen Minimum Zone (OMZ) in order to quantify the past [O 2 ] BW . The data revealed a drop in [O 2 ] BW of 30 μmol/kg at the lower boundary of the OMZ 15 during the last deglaciation. The decrease was largest at the deepest core site whereas differences were not significant closer to the centre of the OMZ. The overall deoxygenation trend started with the onset of Heinrich Stadial 1 and it was first observed at the southernmost core. It was later followed by a distinct drop during the deglaciation in all other cores. A slight increase was observed 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 20 Pacific Ocean,

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 ([O2]BW) during sampling were used as dependant variable. The correlation was significant (R 2 = 0.82; p<0.05) indicating that the foraminiferal assemblages are rather governed by oxygen availability than by the deposition of particulate organic matter 25 (R 2 =0.53; p=0.31). We applied the regression formula to four sediment cores from the northern part of the Peruvian OMZ between 3°S and 8°S and 600 m to 1250 m water depths; 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 [O2]BW change from the Last Glacial Maximum and the Holocene was constrained to 30 μmol/kg at the lower boundary of the OMZ, whereas at shallower depths [O2]BW 30 was relatively stable along the deglaciation. The deoxygenation trend was time-transgressive. It commenced at the southern core, and gradually spread to deeper waters and to the northernmost core location. This pattern indicates a gradual expansion of the OMZ during the last deglaciation, as a result of increasing surface productivity in the Eastern Equatorial Pacific and decreasing advective oxygen supply to intermediate waters off Peru.
1 Introduction 5 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 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 10 Levin, 2004;Paulmier and Ruiz-Pino, 2009). Warmer conditions have contributed to expansion of OMZs during the last decades 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 15 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). 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 20 region used established paleo-proxies such as sedimentary textures (laminations), productivity indicators (Corg, δ 15 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). The 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 25 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 m 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 5 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 . 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 600 and 1250 m water depth. We compiled all available information on living (rose Bengal 10 stained) benthic foraminiferal faunas from the Peruvian margin as 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

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 20 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 both, the availability of 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 25 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; 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 30 dwelling in the first one or two 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 recognised 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 5 could be used to identify spatial changes of the OMZ provided certain threshold values and ranges are considered (Table 1). Here, we consider the following classification: microxic conditions <5 μmol/kg, dysoxic conditions 5-45 μmol/kg, 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; 10 (Schumacher et al., 2007;.

Regional setting
The Peruvian OMZ is one of the most pronounced OMZs in the world (Figure 1; (Paulmier and Ruiz-Pino, 2009), covering the Peruvian continental shelf and upper slope, with its thickest part between 5°S and 15°S and 50 to 750 m water depths (Figure 2; (Fuenzalida et al., 2009). The intensity of the OMZ is dependent on the low 15 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 depths (Montes et al., 2010;Chaigneau et al., 2013). The 20 PCUC is fed by the Equatorial Undercurrent (EUC) and Southern Subsurface Countercurrents (SSCCs; (Montes et al., 2010). Below the PCUC, northward flowing Chile-Peru Deep Coastal Current (CPDCC) carries cold Antarctic Intermediate Waters as a thin layer (AAIW; (Chaigneau et al., 2013). 25 Four sediment cores were considered in this study (M77/2-47-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°S and 9°S and water depths of 600 m and 1250 m around the lower boundary of today's OMZ ( Figure 1 and Figure 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 correlated to stacked standard records (e.g., (Lisiecki and Raymo, 2005) or Antarctic ice cores (EPICA 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 M77/2-59-1: Mollier-Vogel et al. 5 (2013), and for others Erdem et al. (2016), respectively). The results 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.qub.ac.uk/marine/) ranging from 89 to 338 years (Erdem et al., 2016). All ages are expressed in thousands of years (ka) before 1950 AD (abbreviated as cal ka). For this study, 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 10 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 to 20 cc 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 15 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) (Table 3). 20

Surface samples and living benthic foraminifera
Information on the living (rose Bengal stained) benthic foraminifera was compiled from four independent datasets. They comprise 53 samples from the Peruvian continental shelf and slope from water depths of 48 to 2092 m between 1°45'S and 17°28'S ( Figure 2, Table 4). Four of the samples from a transect around 12°30'S were collected in December -January 1998 during Panorama Expedition, Leg 3a, with R/V Melville. Eight

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 core 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 10 (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 15 the results.
Another exception was made for Hoeglundina elegans. The aragonitic test of this species has a low preservation potential . 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 20 reduced accordingly.
Considering the water depths at which the sediment cores were taken, (>600 m) and the depth range of our surface samples (Figure 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 25 species from 35 samples.

Environmental data
Bottom water oxygen concentrations ([O2]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, [O2]BW data were measured on each cruise . We thus considered averaged oxygen concentrations. The [O2]BW of the R/V Meteor M77 stations were taken from  and from a synoptic compilation of all CTD hydrocasts from this area (Schönfeld et al., 2015). The [O2]BW of the R/V Meteor M137 stations were extracted from the expedition dataset (Gerd Krahmann, GEOMAR, pers. comm.). Rain rates 5 of particulate organic carbon (RRPOC; mmol/m 2 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. 10 (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.

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 15 preferred instead of absolute abundances (individuals per cm 3 ) 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 [O2]BW and RRPOC. 20 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 were applied using the Unweighted Pair Group Method (UPGMA) based on a Bray-Curtis similarity matrix. Canonical Correspondence Analysis (CCA) was performed to the same datasets to see the relation between the species, stations and environmental variables. Additionally, multiple regression analysis was applied and their significance was assessed in order to 25 evaluate their reliability for downcore applications. Coefficients and intercept value from the multiple regression analysis was later used as data entry for a polynomial transfer function to calculate past bottom-water oxygen concentrations from foraminiferal assemblages of sediment core samples.

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  μmol/kg), and microxic (<5 μmol/kg) classification was applied to this dataset, 27 samples were classified as 5 microxic, 20 samples were classified as dysoxic and the remaining 6 samples were classified as oxic (Figure 3).
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 [O2]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, 10 however, less diverse as 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 the stations do in the reduced dataset of 16 species and 35 samples that we statistically analysed. The microxic samples showed less diverse assemblages with higher dominances, whereas the dysoxic samples showed higher diversities and a lower dominance. The Fisher alpha diversity indices were lower and 15 dominance was higher at stations under dysoxic conditions where the living faunas were dominated by Bolivinids

Application of [O2]BW estimates to sediment cores
Erosion, reworking and high energetic bottom conditions prevail at the continental slope of the Peruvian margin.
The Holocene from cores 47-2 and 50-4 were 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 5 intervals in all of the sediment cores. In the following, we describe the results of [O2]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, three sediment cores were recovered around or deeper than 1000 m which was the maximum depth for any reliable RRPOC calculations. 10 Core M77/2-47-2 was retrieved from 626 m water depths, and the prevailing oxygen at this location was around Deviations ranged from 9 to 17 μmol/kg. Core M77/2-52-2 was collected from 1249 m water depth. The [O2]BW was 74 μmol/kg at the core location during sampling. M77/2-52-2 is the only core which spans all time intervals considered in this study. In total, 27 samples were analysed, and 170 species were identified of which three were agglutinated. The estimated [O2]BW indicated stable condition during the LGM with values ranging from 52 to 61 25 μmol/kg, which was followed by a decrease from 56 to 46 μmol/kg during HS1, and a further much more distinct decrease during the 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 [O2]BW was 54 μmol/kg during sampling. 30 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 as 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 [O2]BW ranged from 34 to 60 μmol/kg during HS1, from 63 to 18 μmol/kg 5 during the BA/ACR, 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 BA/ACR was 16 μmol/kg in core 20 μmol/kg at core 52-2 and only 6 μmol/kg in core 59-1. Whereas the shallowest core 47-2, indicated that [O2]BW values slightly increased (5 μmol/kg) from 10 HS1 to BA/ACR. The difference between the BA/ACR and the early Holocene [O2]BW was around 11 μmol/kg at core 52-2 and around 24 μmol/kg at core 59-1.

Living benthic foraminifera in relation with OMZ settings
Similar to previous observations from other modern OMZs, benthic foraminifera living in the Peruvian OMZ 15 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 even anoxic conditions (Table 2; Mullins et al., 1985;Schumacher et al., 2007;Piña-Ochoa et al., 2010;Glock et al., 2011;. In the vicinity of the OMZ core, around the boundary, a more diverse assemblage (cluster A; supp. Figure (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 much more diverse assemblage (cluster B; supp. Figure S.2). This assemblage involved mostly Miliolids which are 25 known for their intolerance to low-oxygen (e.g., . Agglutinated species dominated the whole assemblage from these samples which again mirrors rather oxic conditions. The data is in good agreement with previous observations on the low tolerance of agglutinated species to oxygen depleted conditions (<0-20 ml/l = ~9 μmol/kg; Gooday et al., 2000). In the CCA on the living fauna, this Miliolids associated assemblage (cluster B), together with the agglutinated species, showed a positive relation with [O2]BW (supp. Figure S.3). Similar trends of such different assemblages in relation with changing oxygen and organic matter were also observed in the Arabian Sea . 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 -high density to more diverse 5 and cosmopolitan assemblages indicate 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  and Bolivina costata is frequent in the core of the Peruvian OMZ (This study; 10 (Cardich et al., 2012;Glock et al., 2019). 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 since 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 15 Northern Arabian Sea . 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;. Elevated proportions of these two species could be used as extreme-low bottom-water oxygen indicator in downcore records (e.g., McKay et al., 2015;20 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 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). 25

Peruvian margin oxygen history since the LGM
The records revealed a distinct decrease in oxygen during deglaciation in all of the cores from the lower OMZ boundary. Three cores did not cover all the considered time intervals, limiting the spatial delineation of oxygenation changes between the LGM and the Holocene. 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 30 to the BA/ACR, and a slight increase followed by relatively stable concentrations during the Holocene. None-theless, Holocene [O2]BW 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, δ 15 N, redox sensitive elements) which are known to indicate oxygen depleted conditions as recorded by sediment cores from 5 above 500 m depths off Peru. Comparing the differences between our cores, the highest overall decrease was observed in the deepest core, implying that the absolute changes at the lower boundary were larger than in the core of the OMZ. This might be a reason why profound oxygenation changes were not recorded in core 47-2.
Moreover, this core location was at least 100 m shallower during the LGM hence was potentially within the OMZ core. This record also shows that when the OMZ intensifies (or diminishes) the change is profound around its 10 borders and the conditions are rather stable close to its centre. This regional dynamics has also been disclosed with other approaches. For instance; 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). Furthermore, differences and delays were observed in the timing of the decreases among the different cores. In the southern core 50-4, the most distinct decrease started with the onset of the HS1, whereas the decreasing trend did not commence before the end of this period in the slightly deeper core 52-2 ( Figure 6). The amount of the 20 oxygen decrease from HS1 to the BA/ACR was around 15 to 25 μmol/kg in the cores 50-4 and 52-2, whereas it was only 5 μmol/kg in the northern core. At this core, the values raised later by ca. 10 μmol/kg in the BA/ACR to the Holocene transition whereas the decrease was only around 6 μmol/kg in the core 52-2. This indicated a south to north time-transgressive and decreasing trend in [O2]BW. The downcore distribution of benthic foraminiferal species that were not included in the quantification approach, provided further, corroborating evidence. For 25 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 is limited. 30 Besides, large standard deviations observed in the northernmost core 59-1 raised questions about the applicability of the method for the Holocene at this core ( Figure 6). Never-the-less; both Holocene records of core 59-1 and 52-2, suggested a slight recovery of the OMZ that is in accordance with recently published results from the region (Salvatteci et al., 2016;Salvatteci et al., 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 . However, we cannot comment further for the rest of the Holocene [O2]BW trend on the 5 basis of currently available information.  15 Total organic carbon content (TOC (w.%)) from core 52-2 and 59-1 showed an increasing trend with the onset of  (Figure 7; Glock et al., 2018;Mollier-Vogel et al., 2019). At core 59-1, the increase started not earlier than at the onset of the BA/ACR, whereas core 52-2 showed an increase at the onset of the HS1 already. 25 These results are in agreement with our oxygen reconstructions for the same cores. In core 52-2 for instance, the deoxygenation at the end of HS1 followed relatively stable conditions while δ 15 Nsed showed lower values and increased again not earlier than at the beginning of the BA/ACR. This coupling of δ 15 Nsed and bottom-water oxygen suggested a recovery of oxygen depletion in the lower boundary of the OMZ, which is in accordance with the pore-water nitrate reconstruction from the same core . Later in the record both cores 30 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 δ 15 Nsed values.

Comparison with other proxies and records from the region
Other redox and productivity proxies from the region revealed similar trends. Laminated sediments indicating [O2]BW of <7 μmol/kg (e.g., Schönfeld et al., 2015) were observed 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 5 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 was 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 10 relative abundances of the phytodetritus, bloom-feeding Epistominella exigua in core 52-2 during the HS1 and in core 59-1 during the entire deglaciation supported these observations indicating an enhanced organic matter flux to the sea floor (Erdem and Schönfeld, 2017). Never-the-less, 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 15 enhanced surface productivity at the core location. The discrepancy might be a corroborating evidence that other parameters 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 δ 30 Si records (Doering et al.,20 2016) as well as impact of different intermediate water masses and stratification (Bova et al., 2015;Bova et al., 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 -2 yr -1 for this species

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)