The central Peruvian upper margin presents some geomorphological differences: the shelf (200-m isobath) north of 13°S is relatively wide (∼50 km) with a pronounced break to the slope, while south of 14°S is narrow (∼15 km) and with a soft transition to the slope (Figure 1B). Sediments in the central Peruvian upper margin contain fine grains, are rich in organic carbon (Gutiérrez et al., 2006) and are highly reduced (Suits and Arthur, 2000). Sediment records show different structures, but laminations are common (Salvatteci et al., 2014a). Muddy laminated areas are found between Huacho and Callao (11 – 12°S; Reinhardt et al., 2002), but sedimentary records retrieved in the OMZ core off Pisco (14°S) are better laminated (Gutiérrez et al., 2006; Salvatteci et al., 2014a). ²¹⁰Pb sedimentation rates vary from 0.4 mm/y to 4 mm/y in the shelf and slope (Levin et al., 2002; Gutiérrez et al., 2009; Salvatteci et al., 2018). Figure 1B shows the main three equatorial currents arriving the Eastern Tropical Pacific. The Equatorial Undercurrent (EUC) and the primary and secondary Southern Subsurface Countercurrents feed the Peru Undercurrent (PUC) (Montes et al., 2010). Figure 1B also shows the two main subsurface currents off Chile and Peru, the PUC and the equatorward Chile-Peru Deep Coastal Current (CPDCC, Chaigneau et al., 2013) with an upper limit at ∼500 m is depicted as well.

To assess the oxygen variability in the last decades, we used instrumental data of dissolved oxygen by Winkler titration in situ analysis (Strickland and Parsons, 1968; Grasshoff et al., 1999) from multiple cruises carried out by the Peruvian Institute of Marine Research (IMARPE) and from the World Ocean Database for the 1960 – 2017 period (WOD 2018) (Figure 2A). In order to avoid bias we do not include information available in the last decades by other instrument (e.g., CTD, Argos). Quarterly data from 1960 to 2010 between 11 and 15° S for the oceanic area (30 – 100 nautical miles) were analyzed in different water column depths (60, 80, 100, 150, and 200 m). These depths were chosen after evaluating the dissolved oxygen (DO) mean profile for central Peru (Figure 2B). Data of DO were downloaded from the WOD 2018 for the 1960 – 2017 period and is available at https://www.nodc.noaa.gov/OC5/woa18/. For the depth time series, data with >3 replicates were used in the analysis. Afterward, measurements were integrated to show semi-annual variability. A regression analysis was performed for the time series to evaluate significant temporal trends in the time series.

Soutar box cores were collected off Callao (B0405-13, 12°00′S, 72°42′W, 184 m) and off Pisco (B0405-06, 14°07′S, 76°30′W, 299 m) in May 2004 on board of the R/V José Olaya Balandra (Figure 1B). Hereafter both cores will be referred as B-13 and B-6, respectively. The processing of the box cores is described in detail in Gutiérrez et al. (2006). The length of B-13 and B-6 is ∼78 and ∼74 cm, respectively. Both box cores were longitudinally sectioned in six slabs (I – VI) and each of those were used to determine physical characteristics (e.g., X radiography), chronology and geochemical and biogenic proxies. Subsampling for proxies was done following the laminae or bands in the record (Gutiérrez et al., 2006; Morales et al., 2006). Geochronology in both records was solved with radiocarbon dating of bulk organic sedimentary carbon and ²¹⁰Pb measurements (Gutiérrez et al., 2009; Supplementary Material therein). The last ∼33 cm (B-13) and ∼34 cm (B-6) represent the period from ∼1830 to 2004 (∼180 years). This time period consisted of a total of 37 and 41 subsamples (0.25 – 0.8 cm of thickness) in B-13 and B-6, respectively. A third box core nearby B-6 is considered in this study: B0506 – 14, 14°27′S, 76°43′W, 301 m (hereafter referred as B-14). The subsampling and processing of B-14 is described in detail in Salvatteci et al. (2014b). X-ray images of the three box cores are shown in Supplementary Figure S1A. The three box cores show laminated sections and the period from ∼1830 to 2004 starts after a shift in physical properties (e.g., dry bulk density and percentage of calcite, Supplementary Figure S1). The Pisco box cores (B-6 and B-14) are located in the core of the SE Pacific OMZ (Figure 1A,C) while the Callao box core (B-13) is closer to the upper limit of the OMZ.

Pre-treatment for box core subsamples for foraminifera is described in Salvatteci et al. (2012). These subsamples were treated with oxygen peroxide (30%) and sodium pyrophosphate under temperatures of 40–50°C for 5–10 minutes to disaggregate the organic matter. Residual material was later sieved through sieves of 355, 125, and 63 μm mesh size. The material retained in the 125 and 63 μm sieves was dried for about an hour at 40° C. Foraminifera were sorted and counted at the lowest taxonomical level, if possible, in the 63 – 125 μm size fraction under a stereoscopic microscope. Counts from this fraction were then integrated to the counts from the >125 μm size fraction published in Morales et al. (2006) and Gutiérrez et al. (2009). Splitting of large subsamples was performed using an Otto micro-splitter. A minimum of 300 tests was sorted in the splits. We calculated diversity indices on the raw abundance of foraminifera. The Shannon-Wiener index (H’) was calculated as an indicator of diversity; while the Pielou index (J’) was used to show how similar the abundances of the species are in a sample. Foraminiferal standing stocks were expressed in number of individuals per gram of dry sediment (#ind. g^-1). Species with percentages > 0.5% in at least three samples (5% of the total number of samples) were included in subsequent statistical analyses. In addition to counts of benthic foraminifera in the box cores, counts of planktonic foraminifera from the >125 μm fraction (published in Morales et al., 2006) were used to calculate the planktonic to benthic ratio (P/B) to assess calcite preservation in the sediment. We only applied the P/B in the >125 μm fraction as planktonic foraminifera are well preserved here.

In order to understand the foraminiferal assemblages in the cores, we calibrated the living foraminiferal stocks with oxygen and sedimentary redox gradients. For this purpose, we compiled abundance data of the main calcareous benthic foraminifera off Callao and Pisco. We complemented the dataset from Cardich et al. (2015) with new samples of living calcareous foraminiferal stocks from 2012 to 2016. These samples were taken during IMARPE surveys and during the Meteor 92 cruise. We also include other studies from the upper margin of central Peru: Pérez et al. (2002) and Mallon et al. (2012) to better represent the fauna zonation of benthic foraminifera.

We compiled published geochemical information for B-13, B-6, and B-14 to assess the temporal variability of benthic redox conditions (trace metals), water column denitrification (δ¹⁵N) and export production (TOC and biogenic silica). The methodology and data for the geochemical (Mo and δ¹⁵N) and organic concentrations [total organic carbon (TOC), biogenic silica] for B-13 and B-6 are shown in Gutiérrez et al. (2009). Data for B-14 (TOC, biogenic silica, Mo, Re and δ¹⁵N) were published in Salvatteci et al. (2014b, 2018). Subsampling for geochemical proxies in B-13 and B-6 was done every centimeter for trace elements and every two centimeters for TOC, biogenic silica and δ¹⁵N (Gutiérrez et al., 2009). Subsampling resolution was higher for B-14 (Salvatteci et al., 2014b, 2018). In cores B-13 and B-6, TOC was determined from total carbon measurements with a Thermo Electron CNS elemental analyzer, corrected for carbonate content (Gutiérrez et al., 2009) while in B-14 it was determined by RockEval (Salvatteci et al., 2014b). δ¹⁵N in B-6 and B-13 was determined by mass spectrometry after acidification at the Department of Geosciences, University of Arizona, United States. δ¹⁵N analyses in B-14 were measured on a continuous flow gas-ratio mass spectrometer at ALYSES laboratory (Bondy, France). Trace metals’ concentrations were analyzed by ICP-MS after hot-plate acid digestion in Polytetrafluoroethylene vessels and elimination of organic matter and removal of silica using acid treatments (HF, HNO3− and HClO4−). Besides Molybdenum, Uranium concentrations were measured in B-13 and B-6 and were not published before. Rhenium was only measured in B-14. The detrital (background, usually lithogenic) and authigenic (the excess part) fractions of all three metals were estimated. Authigenic Mo is a proxy for temporal to permanent sulfidic conditions (Scholz et al., 2011). Re and U are enriched under less reducing conditions (Algeo and Tribovillard, 2009), and can be used as a proxy for anoxic, but non-sulfidic, conditions. The enrichment factor (EF) of Mo, Re, and U was calculated using background values of andesite as an adequate representation of detrital material in the Peruvian margin. Mo, Re, and U to aluminum mass ratios in andesite are: 0.25 × 10^-4, 1.93 × 10^-9, and 0.34 × 10^-4, respectively (GEOROC database). The EFs were calculated with the formula: EF_{element x} = (X/Al)_sample/(X/Al)_andesite (Tribovillard et al., 2006). Additionally, Mo/U (for B-13 and B-6) and Re/Mo (for B-14) ratios were used as indicators of redox states.

We performed two principal component analysis (PCA) in order to detect structure in the relationships between variables. The first PCA was performed, in both cores, on a matrix of species data with a prior Hellinger transformation (squared root of relative abundances). Foraminiferal principal components were then compared with living assemblages to reconstruct past redox conditions. The second PCA was performed on a matrix of both foraminiferal and geochemical data. This PCA was used to understand the relationship between the foraminiferal assemblages and the different proxies for environmental changes. A varimax rotation was applied to the loadings. Finally, a non-parametric r Spearman correlation was applied to both foraminiferal and geochemical data. For multiple comparisons, the probability level was corrected by dividing the probability level α (p < 0.05) by the number of tests performed (Glantz, 2002). All statistical analyses were conducted using R statistical software version 3.5.1.

Water column mean DO is highly variable down to ∼125 m depth (Figure 2B). Below this depth, DO values are homogenous, more stable and always < 0.5 ml l^-1. We chose data from 60, 150, and 200 m (Figure 2C–E) as they are in the depth gradient from the oxycline to the OMZ core (Figure 2B). The trend of DO values from 1960 to 2010 at 60 m water depth (Figure 2B) differed from the trends at 150 and 200 m (Figure 2C,D). At 60 m, oxygen values were usually > 1.0 ml l^-1, showed high variation and amplitude, and there was not a clear trend in the time series. On the other hand, at 150 and 200 m, oxygen values were usually below 0.5 ml l^-1 and presented little variation. Both 150 and 200 m time series presented a slight positive trend from 1960 to 2010 (p < 0.01). Both time series had an estimated increase of ∼0.02 ml l^-1 per decade.

The abundance of benthic foraminifera was higher in Callao compared with Pisco, but both sites display similar temporal trends (Figure 3A). Densities of benthic foraminifera ranged from 2,504 to 73,660 ind. g^-1 in B-13 and from 117 to 24,498 ind. g^-1 in B-6. Four samples of B-6 were discarded from the analysis as they presented low unrealistic values (1 – 15 ind. g^-1). Both records presented similar patterns until ca. 1960, with three high peaks in B-13 (1834, 1847, and 1861) and two peaks in B-6 (1843 and 1856) and an increasing trend during the 20th century to reach the maximum abundance in ca. 1950. In the last decades of the records, densities in B-6 decreased, whereas densities increased in B-13 until 1990.

The diversity indices and the P/B ratio presented temporal and spatial changes (Table 1). A total of 68 species was identified in both records and a higher species richness was observed in Callao than in Pisco (Table 1). The species richness in B-13 (Callao) ranged from 14 to 35 species, while in B-6 (Pisco) it ranged from 5 to 21. The lowest diversity was found between ca. 1865 and ca. 1885 in B-13 (Table 1). The diversity then increased from 1910 to 2004 in B-13. The number of species did not change downcore in B-6. Species diversity (H_S) ranged from 1.3 to 2.1 in B-13 and from 0.3 to 1.9 in B-6. In both records, the Pielou’s evenness index (J’) increased from ca. 1835 to ca. 1870, then decreased until 1950 – 1960. In the last period, H_S and J’ increased again in B-6 and remained homogenous in B-13. Finally, P/B ratio values were higher in B-13 than in B-6 from the beginning of the records to ca. 1910 (Figure 3B), while from 1910 to 2004, the P/B was similar in both records. After a rise in both records from 1920 to 1945, the values rapidly decreased and remained low in B-6, while those in B-13 did it a decade later.

Temporal changes in the relative abundance of benthic foraminifera can help to determine the state of sediment redox conditions given that each species have a preferred habitat. The foraminiferal fauna in both records was greatly dominated by Bolivina seminuda and Buliminella tenuata (>50% of total abundance, Figure 3C,D). Both species presented a decadal to multidecadal negative co-variation (B-13: Spearman’s ρ = -0.79, p < 0.001; B-6: Spearman’s ρ = -0.93, p < 0.001), especially after ca. 1875 A.D. (Figure 3C). Bolivina costata (Figure 3E) and Nonionella auris (Figure 3F) presented high densities between ca. 1820 and ca. 1865 in B-13, but were not abundant in B-6. Bolivina pacifica showed an increasing trend from ca. 1865 to ca. 1960 in both records (Figure 3H). Other species presented different patterns between records (Figure 3I–M). The relative abundance of Suggrunda eckisi and Bolivina plicata was higher in the XX century compared with the XIX century. Buliminella elegantissima showed high relative abundances in some time periods in both records (Figure 3I). Finally, the presence of Virgulinella fragilis was rare and low abundances were recorded only in the last 40 years of the records (Figure 3G).

The mean spatial distribution of living benthic foraminifera in sediments between 45 and 300 m depth is shown in Figure 4A. The selected species presented a gradual zonation across the upper margin: some species were more abundant near the coast (B. costata and N. auris), some in the outer shelf (B. tenuata, B. elegantissima, and Epistominella sp.), and others thriving in the upper slope (B. seminuda, B. pacifica, and others).

The PCA results, based on the different species of benthic foraminifera, indicate that most of the variance of the data set can be explained by two principal components (Figure 4B,C). In Callao (B-13) and Pisco (B-6), 56.9 and 66.8% of the total variance of the downcore records was explained by two principal components. The PC1 in both B-13 (36.8%) and in B-6 (58.8%) was explained mainly by Bolivina seminuda and Buliminella tenuata (Figure 4B,C). PC1 resembles the negative covariation between B. seminuda and B. tenuata as observed in Figure 3B,C. In B-13, Epistominella sp. also contributed to PC1 (loading factor = 0.45; Figure 4B). PC2 in B-13 (20.1%) and in B-6 (8%) reflected the contrast of B. pacifica to different assemblages of species (Figure 4B,C). For B-13, B. pacifica and B. tenuata showed opposite loadings in comparison to B. costata. For B-6, loadings of B. pacifica were opposed to those of B. elegantissima and Epistominella sp. These observations are in agreement with the relative position of the sediment cores and the zonation in the recent samples (Figure 4A). In B-6, B. costata is absent given that it is a coastal species, while B. elegantissima and Epistominella sp. are present as they thrive in the shelf.

The results of the Varimax rotated PCA using geochemical and foraminifera-based proxies in both B-13 and B-6 are shown in Figure 5. Values for TOC, Si, δ¹⁵N, Mo EF, and Mo/U were used along the selected species in Figure 4 and the P/B ratio. In B-13 (Figure 5A), TOC and B. pacifica were opposed to Mo/U, Mo EF, Si, B. costata and N. auris (Varimax rotated PC1 – RC1). On the other axis, B. tenuata was contrasted to B. seminuda and the foraminiferal assemblage of deep sites (Varimax rotated PC2 – RC2). A similar pattern is observed in B-6, but with an inversion of the components: B. seminuda was contrasted to B. tenuata and other outer shelf species (RC1); while Mo/U, Mo EF, and Si were distributed in an opposite direction of B. pacifica (RC2). Additionally, TOC was distributed closer to δ¹⁵N in Pisco (B-6) than in Callao (B-13).

The new observations of living foraminiferal assemblages in addition to the findings in Cardich et al. (2015) helped to determine the temporal changes in sedimentary redox conditions. The distribution of the main species across the upper margin of the central Peru region presented a bathymetrical zonation responding to gradients of OM quality and sulfide concentration. A first group of coastal species (B. costata and N. auris) is more abundant toward the inner shelf (Figure 4A), where sediments are typically sulfidic and rich in fresh phytodetritus (Cardich et al., 2015) but are also exposed to interannual oxygenation events (Gutiérrez et al., 2008). A second group (Buliminella tenuata, B. elegantissima, and Epistominella sp.) is present across the margin, but concentrates in the outer shelf, which is less exposed to oxygenation events, but typically presents sediments with less labile OM, low sulfide tenors and Thioploca spp. mats. In third place, there is a group of species more abundant in the upper slope where the low oxygen values are not affected by high frequency temporal variability (OMZ core). In this area, the bottom waters are suboxic (dissolved nitrate is present and nitrite is high, Sommer et al., 2016) and the sediments are postoxic (anoxic, non-sulfidic) and rich in preserved OM (Dale et al., 2015). Within this last group, Bolivina seminuda shows a steep growing dominance from the inner shelf toward the upper slope, along the increase of nitrate and nitrite in the overlying waters. This pattern is consistent with the nitrate or nitrite respiration of B. seminuda (Piña-Ochoa et al., 2010; Glock et al., 2019). Bolivina pacifica is in general less abundant and follows B. seminuda’s distribution, but is absent in the inner shelf sulfidic sediments. Finally, B. plicata, Cancris spp. and P. subperuviana are only present in the upper slope.

Based on the observations from the foraminiferal ecological zonation (Cardich et al., 2015), the species composing the first two principal components of both B-13 and B-6 are part of the main living assemblages of the reduced shelf sediments off Callao and Pisco (Figure 4). Moreover, the foraminiferal PCs are in closely agreement with the faunal zonation. This assures the use of the PCs as proxies for biogeochemical gradients (e.g., redox) in our cores. Given that the variance of the PCs can be explained by no more than two or three dominant species (Figure 4B,C), we use these species ratios in the sediment cores to reconstruct the redox conditions. The use of ratios simplifies this interpretation of the PCs which is in agreement with the abiotic – biotic relationships (Figure 5). Thus, PC1 is expressed as the B. tenuata to (B. tenuata + B. seminuda) ratio [Bt/(Bt+Bs)] for both box cores. On the other hand, PC2 can be expressed as the Psp to (Aspx + Psp) ratio, Asp representing the anoxia assemblages: B. costata for the inner shelf (B-13) and Bu. elegantissima and Epistominella sp. for the outer shelf (B-6); and Psp representing the postoxia species B. pacifica. We did not use B. tenuata in the Psp/(Aspx + Psp) because this species was already present in the first ratio (Figure 4).

Principal component analysis results on geochemical and foraminiferal data (Figure 5) and Spearman correlations (Table 2) in the box cores confirm the above ecological interpretation. Moreover, there is a clear relationship between the different foraminiferal assemblages with reducing conditions in RC2 of both box cores (Figure 5A,B). This indicates that the sedimentary geochemical zonation in the margin is consistent and that more reducing conditions (e.g., sulfidic) are typically found near the coast. From the faunal distribution (Figure 4A) and the PCA results (Figure 5), an association between B. tenuata and mild reducing conditions is detected.

The co-variation of B. tenuata and B. seminuda may be associated with the water column denitrification (WCD). When comparing the B. tenuata to (B. tenuata + B. seminuda) ratio in B-13 and B-6 (Figure 6E,L) to the δ¹⁵N records (Figure 6B,I), a slight co-variation between the dominance of B. tenuata (B. seminuda) and stronger (weaker) WCD is noticed. This connection might be explained by changes in local export production to the sediments. A greater organic input is translated in prevailing anoxic sediments with no nitrate, low sulfide tenors and more labile OM. The thriving of B. tenuata is evident in this type of condition, which is typical of outer shelf sediments (Figure 4A; Cardich et al., 2015). On the other hand, postoxic sediments with preserved OM and available dissolved nitrate are found under lower organic input. B. seminuda and associated species thriving in the upper slope appear to be associated with availability of dissolved nitrate in the bottom waters (Glock et al., 2019; Cardich et al., unpublished results). The association of B. tenuata with δ¹⁵N exists in B-13 (Figure 5A) but not in B-6 (Figure 5B), giving more evidence that the intensity of WCD over the shelf and its effect on reducing conditions is mediated by local export production.

On the other hand, the ratio of shelf versus slope species is associated with redox extremes and environmental stability. Thus, lower (higher) values of the Psp/(Psp + Asp) ratio for Callao indicates sulfidic (mild postoxic) sediments with high labile (preserved) OM. Lower (higher) values also indicate exposition to coastal oxygenation events (stable bottom water suboxia). Likewise, lower (higher) values of the Psp/(Psp + Asp) for Pisco indicates light sulfidic (postoxic) sediments with labile (preserved) OM and less (more) stable conditions. These local disparities are a product of the bathymetric difference between the study sites. B-13 (180 m, outer shelf) is more exposed to coastal oxygenation processes than B-6 (300 m, upper slope), reflecting a biogeochemical gradient for coastal species. Meanwhile B-6 is more prone to register changes in the OMZ core of central Peru.

From the patterns of geochemical proxies along the records (Gutiérrez et al., 2009; Salvatteci et al., 2014b; Supplementary Figure S2), we differentiated three distinct time periods after ca. 1835:

Benthic foraminiferal ratios (Figure 6E,F) and rotated principal components (Figure 6G,N) for B-13 and B-6 also indicated the variability of redox and OM quality. The first time period (ca. 1835 to ca. 1875) showed a dominance of anoxia-tolerant foraminiferal assemblages inhabiting reduced sediments off Callao and Pisco, which evolved to sulfidic conditions off Callao. Average reducing conditions off Pisco (B-6) in the first period were as intense as in Callao (B-13), however, strong sulfidic conditions are observed in B-14. This suggests that the cause of reducing sediments come from the neritic water column (i.e., local productivity). The presence of diatom bands in B-13 and B-14 in this period (not visible in B-6 because of a hiatus; Salvatteci et al., 2014a) and peaks of biogenic silica content indicate that massive algal blooms were frequent. Algal blooms in this period might be enhanced by a prevailing water column stratification, inferred from the warm sea surface temperatures (Gutiérrez et al., 2011) and weak surface winds (Briceño Zuluaga et al., 2016). The second time period was represented by a (multi)decadal oscillation of the WCD and organic input (co-variation of B. seminuda and B. tenuata) and a sustained relaxation of reducing conditions to non-sulfidic sediments off Callao and Pisco (Figure 6). Weakening of sulfidic conditions (settling of postoxic sediments) continues during the last time period. However, the WCD off Pisco shows a strong weakening trend. Being situated in the shelf, the Callao site (B-13) is influenced mainly by coastal productivity fluctuations and the Pisco site (B-6) is more exposed to the oceanic variability. This implies that the variations of water column δ¹⁵N is affected by both local and regional processes. Thus, WCD multidecadal changes off Pisco can be attributed to variations of the OMZ intensity, since this signature is more marked off Pisco than it is off Callao. From the latter, we can attribute the variability of the OMZ off Pisco to the variability in subsurface circulation.

We interpret the OMZ weakening (from ca. 1875 to 2004) as a result of ventilation arriving either from the equator or from the south, or as an interplay of both sources. Equatorial oceanic circulation linked to the Peruvian upwelling system is complex. The EUC fuels the PCUC with relatively oxygen-richer waters in comparison to the primary and secondary Southern Subsurface Countercurrents (p and s in Figure 1C; Montes et al., 2014). A major and sustained contribution of the EUC to the PCUC would ventilate the Peruvian OMZ. However, ultimately “the water mass properties within the PCUC core result from a delicate balance between different sources in terms of oxygen” (Montes et al., 2014). Overall, the oxygen supply is low to the PCUC, but it contributes to the ventilation of the OMZ (Karstensen et al., 2008). An intensification of the EUC took place since 1900 according to a SODAS reanalysis in the central Equatorial Pacific (Figure 7A; Drenkard and Karnauskas, 2014). The velocity of the EUC increased in 16 – 47% per century since the half of the 19th century (Drenkard and Karnauskas, 2014). Nevertheless, a reconstruction of the variability of the Tsuchiya Jets in the same time period is needed. It is known that the PCUC is the main nutrient source for the coastal upwelling off Peru and northern Chile, favoring primary productivity and subsequent export production. As the latter processes presented an enhancement since ca. 1900 (Gutiérrez et al., 2011; Salvatteci et al., 2018), dissolved oxygen arriving in the PCUC might have been rapidly consumed. The equatorward subsurface current CPDCC may be another source of ventilation. At about 500 m depth (Chaigneau et al., 2013), the CPDCC transports relatively fresh and cold Antartic Intermediate Water (AAIW) northward (Pietri et al., 2014). Pietri et al. (2014) suggested that deep equatorward flow off Peru may be modulated by coastal trapped waves at intraseasonal time scales, enhancing the transport of AAIW to the tropics. This suggests that Equatorial circulation might drive intermediate circulation. However, this process needs more research at (multi)decadal time scales.

At this point, our observations suggest that the variability of subsurface oxygenation over the Peruvian margin at (multi)decadal to centennial scales might be a result of an interplay of subsurface ventilation and local productivity. The PCA performed on the geochemical and foraminiferal data (Figure 5) showed a close association between the proxies of anoxia (Mo EF and Mo/U) and Si contents. This correlation is stronger in B-13 (outer shelf), indicating that local siliceous productivity drives the bottom oxygenation variability. Besides, IMARPE instrumental data shows evidence of a slightly subsurface oxygenation at 200 m depth since 1960 (Figure 2; red lines in Figure 7C). In the same time period, a positive trend in chlorophyll-a contents has been reported within 100 km off the coast (13.5 – 14.5°S) (Gutiérrez et al., 2011). Oxygenation might have been more intense at intermediate waters, which could explain the reduction of the OMZ strength.

We compared our records with the main regional climatic indices to assess their influence on the OMZ variability. In that way, we used the Interdecadal Pacific Oscillation (IPO; Henley et al., 2015), the ENSO 1+2 index and the core velocity of the EUC from the SODA reanalysis (Drenkard and Karnauskas, 2014). The IPO integrates sea surface temperature anomalies of the Pacific Ocean with weak trade winds during its warm (positive) phase and strong trade winds during its cool (negative) phase. The IPO resembles a multidecadal “El Niño-like” pattern of climate variability (Henley et al., 2015). The ENSO 1+2 index is used to show how ocean-atmospheric circulation in the northern part of the Peruvian upwelling system impacts the study area. The velocity of the EUC is also shown as a possible consequent ventilation process driven by both IPO and ENSO 1+2 indices.

Deutsch et al. (2011) and Ito and Deutsch (2013) related the (multi)decadal variability of the Tropical Pacific OMZ to the natural variation of biogeochemical and physical processes. The authors explained that the suboxic zone of this region expands during La Niña events or during the cold phase of the Pacific Decadal Oscillation (PDO) by means of a shallower thermocline that promotes nutrient supply to the surface, a higher organic flux to the seafloor and an elevated respiration. The opposite occurs during EN events or the warm phase of PDO. These environmental changes might be associated with the δ¹⁵N records and the foraminiferal ratios [Bt/(Bt+Bs)] off Callao and Pisco in some degree. However, because of the difference in resolution between the climate indices and the proxies in the box cores, possible correlations might be masked. We suggest that additional time series analysis, such as cross-correlations, is needed to better infer cause-effect relationships at a finer scale.

As discussed above (see Understanding the Benthic Foraminiferal Assemblages in the Cores), a possible indirect link between the co-variation of B. tenuata and B. seminuda and WCD exists. An intensified WCD/OMZ over the margin during periods of frequent and strong LN events is associated with a high input of silica to the sediments. On the other hand, when strong EN events are frequent, the OMZ is deepened, WCD decreases, dissolved nitrate is more available in the bottom water and the upper margin sediments have less labile organic carbon. These oceanographical changes and their impact in the sedimentary geochemistry determine the dominance of either B. tenuata or B. seminuda. This suggests that there might be a link between species dominance and ENSO and IPO, possibly mediated by changes in the exported organic carbon and/or in oxygenation/denitrification.

Records of oxygenation/redox-proxies from other East Pacific upwelling zones are included in Figure 7: California (NE Pacific; Deutsch et al., 2014) and Concepción (SE Pacific; Srain et al., 2015). The weakening trend of the OMZ from a multidecadal to centennial scale is a feature present in both the NE Pacific and the SE Pacific (Figure 7B–E). Deutsch et al. (2014) conclude that the main factor controlling OMZ intensity in the NE Pacific is Walker Circulation through variation in thermocline depth as a natural variability not related to global warming. These records show a recent OMZ expansion since the mid-1990s (Figure 7B) because of the strengthening of the Walker Circulation (England et al., 2014). Deutsch et al. (2014) also suggest that when a new declining trend of trade winds settle, the NE Pacific OMZ will contract again. The reactivation of the OMZ expansion in the last decade is apparent off Concepción, Chile (Figure 7E), but not that evident as it is in California (Figure 7B). The sediments in Concepción shelf display a more reducing condition after 2000 (Srain et al., 2015). Despite the Concepción site being more representative of the variations of the oxycline and local productivity because of its location (∼100 m depth), it shows the general picture of the regional OMZ variations. The OMZ off Peru (based on the Pisco records), during the last decades, does not present the same pattern as Concepción and California given that suboxic waters are maintained in the core of the OMZ (Figure 7C,D).

The slight oxygenation trend near the equator (central Peru) as observed in the oxygen data and the paleo data (Figure 7C,D), and the OMZ expansion in the subtropics during the last 10 – 15 years appear to be regulated in some degree by EUC intensification. Subsurface equatorial source waters of the EUC are oxygen-rich for Peru but simultaneously result in oxygen-poor for higher latitudes (24.2 – 33.5°S, California and 36°S, Concepción). However, a shoaling trend of the oxycline off central Peru during the last decades has been determined indicating a vertical OMZ expansion (Bertrand et al., 2011; Espinoza-Morriberón, 2018). An alternate explanation for our results is deeper ventilation over the margin driven by intermediate circulation (e.g., CPDCC) occurring at the same as oxycline shoaling. The expansion of the volume of hypoxic waters off California during the last decades (since 1960s; Bograd et al., 2008) has been associated with changes in the source of water masses (EUC and northern subtropics) which affected ventilation in the basin (Bograd et al., 2015). This scenario might be the same at high subtropical latitudes in the SE Pacific. It must be noticed also that the Concepción record is shallower and might show a signal of local processes. Other coastal records from Callao (68 m, unpublished) and from Chile [100 m deep, Mejillones Bay (23°S); Díaz-Ochoa et al., 2011] exhibit a strong trend toward reducing conditions in the last decades. These results indicate that the progressive local productivity increase (among other factors) controls coastal hypoxia.