Multidecadal Changes in Marine Subsurface Oxygenation Off Central Peru During the Last ca. 170 Years

Introduction

Massive oceanic midwaters with permanent oxygen-deficiency are found in the open ocean and in regions influenced by coastal upwelling (e.g., East Pacific, Arabian Sea). These areas are known as oxygen minimum zones (OMZ) and are produced by the combination of high oxygen demand, sluggish circulation, and low-oxygen source waters (Helly and Levin, 2004; Karstensen et al., 2008). The OMZs show strong variability at multiple time scales, ranging from interannual (Cabré et al., 2015; Graco et al., 2017) to geological scales (Jaccard and Galbraith, 2012). OMZs play a relevant role in the global cycles of carbon and nitrogen, and related greenhouse gasses (Paulmier and Ruiz-Pino, 2008) and constrain productivity and biodiversity (Levin et al., 2009).

Over the last decade, several studies have reported global ocean deoxygenation trends since the middle of the 20^th century (Bograd et al., 2008; Stramma et al., 2008; Keeling et al., 2010; Horak et al., 2016; Schmidtko et al., 2017). The main drivers of deoxygenation in the upper water column are increased stratification, reduction of oxygen solubility and acceleration of oxygen consumption, all due to global warming (Keeling and Garcia, 2002; Breitburg et al., 2018). As a result, an expansion of OMZs have been documented (Stramma et al., 2008). Another type of oxygen loss occurs in coastal waters, where dead zones are caused by human-caused eutrophication. These zones have increased in number (Breitburg et al., 2018) and their interaction with expanding OMZs represent a global threat for ecosystems and coastal fisheries (Stramma et al., 2010). Global deoxygenation is expected to continue in the future with climate change, as most of numerical models predict (Bopp et al., 2013). However, several global models fail to reproduce recent regional deoxygenation trends and differ in the prediction of future regional trends, particularly in the tropics where most of OMZs are hosted (Schmidtko et al., 2017).

The Peruvian upwelling system is one of the most productive areas of the world ocean (Chávez and Messié, 2009) and is associated with a permanent, intense and shallow OMZ (Gutiérrez et al., 2008; Graco et al., 2017). The source of low oxygen waters and nutrients that feeds coastal upwelling is the Peru–Chile Undercurrent (PCUC; Montes et al., 2014). The OMZ is wider off central Peru and impinges the continental margin, generating strong biogeochemical gradients in the surface sediments (Gutiérrez et al., 2006). The Peruvian upwelling system is also subjected to a significant interannual variability linked to equatorial Kelvin waves and the El Niño-Southern Oscillation (ENSO) that affects upwelling, nutrient availability and coastal ventilation (Chavez et al., 2008; Espinoza-Morriberón et al., 2017). The oxycline and the coastal subsurface oxygenation off Peru are modulated by this interannual variability and temporal regimes can be defined by the occurrence, frequency and intensity of El Niño (EN) events and Kelvin waves (Gutiérrez et al., 2008; Graco et al., 2017). The seasonal cycle is driven by local respiration, vertical mixing/stratification and mesoscale circulation as eddies and filaments (Thomsen et al., 2016; Vergara et al., 2016; Graco et al., 2017).

Several geochemical and biogenic proxies are used to infer geochemical conditions. Redox-sensitive metals are usually used to infer paleo-redox conditions. Poorly oxygenated overlying waters and the flux of organic matter control sedimentary reducing conditions (McManus et al., 2006). Metal enrichment (excess part from a background) occurs during or after deposition, and each metal presents different sensitivities in the redox gradient (Tribovillard et al., 2006). In general, higher metal enrichment in the sediments suggest more reducing conditions (Tribovillard et al., 2006). Scholz et al. (2011) summarized the behavior of main redox-sensitive metals in the Peruvian margin. Among them, molybdenum and uranium show different sensitivities to redox conditions that can be used to depict sulfidic (Mo) to suboxic (U) conditions. Rhenium is another metal that behaves as U (Colodner et al., 1995; Crusius et al., 1996). Benthic foraminiferal assemblages represent another reliable proxy for paleo-redox reconstructions. The living community is abundant in marine sediments, is highly sensitive to environmental changes and their calcareous tests are preserved in the sediment record (Gooday, 2003). Calcareous benthic foraminifera are abundant in OMZ sediments such as off California (Phleger and Soutar, 1973), Peru (Mallon et al., 2012; Cardich et al., 2015), and in the Arabian Sea (Caulle et al., 2014, 2015). Calcareous tests of dead foraminifera are typically well preserved in low-oxygen sediments, but oxygenation events might promote carbonate diagenesis (Jahnke et al., 1997). In the Peruvian margin, the foraminiferal assemblages are associated with particular biogeochemical states defined by the interplay of porewater sulfide concentration and the quality of sedimentary organic matter (Cardich et al., 2015).

Large changes in OMZ intensity have been reconstructed in the Eastern Tropical South Pacific during the Late Quaternary (Gutiérrez et al., 2009; Scholz et al., 2014; Salvatteci et al., 2016). Based on multiproxy records (e.g., trace metals, foraminifera, organic carbon, δ¹⁵N) from two sites located in the current OMZ core (dissolved oxygen ∼0.2 ml/l), centennial-scale biogeochemical regimes were evidenced (Gutiérrez et al., 2009). For instance, a major shift from low productivity and suboxic bottom waters to nutrient-rich, oxygen-depleted waters is evidenced toward the ending of the Little Ice Age period (early 19th century). The current low oxygen levels in the OMZ off Peru were established after a large climatic reorganization in the tropical Pacific involving the Intertropical Convergence Zone, the South Pacific Subtropical High and the Walker circulation (Sifeddine et al., 2008; Gutiérrez et al., 2009; Salvatteci et al., 2014b). Changes in this ocean-atmosphere circulation controls the Humboldt system variability at longer time-scales (Salvatteci et al., 2016, 2019).

Here, we studied the decadal to multidecadal subsurface oxygenation variability from the central Peruvian continental margin using a multiproxy approach combining multiple proxies from sediment records and oxygen measurements. First, we assessed the oxygenation trends at 60, 150, and 200 m depth based on instrumental data for the 1960 – 2010 period. Second, we assessed benthic foraminiferal records in order to depict sedimentary redox changes associated with variations in OMZ intensity throughout the 180-year record. The interpretation of benthic foraminiferal indicators is achieved through a calibration with ecological distribution data shown by Cardich et al. (2015) and complemented with new data presented in this work. Our goal is to determine the relative importance of large-scale (i.e., ocean-atmosphere) versus local processes (i.e., export production) driving sub-surface oxygenation off Peru. To do so, we compared our results to other geochemical proxies analyzed in the same box cores (described in detail in Sifeddine et al., 2008; Gutiérrez et al., 2009, 2011; Salvatteci et al., 2014b, 2016; Briceño Zuluaga et al., 2016) and to regional paleo-records.

Materials and Methods

Study Area

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.

Figure 1. Site location and description of box cores used in this study. (A) Mean dissolved oxygen concentration (ml l^-1) at 400 m depth at the Tropical East Pacific (data available in World Ocean Database 2013). (B) Map of central Peru indicating the location of B0405-13 (B-13; 12°00’ S, 77°42′ W, 184 m), B0405-06 (B-6; 14°07′ S, 76°30′ W, 299 m) and B0506-14 (B-14; 14°27′S, 76°43′W, 301 m). (C) Dissolved oxygen concentration in an oceanographic section along the coast of America (black line in inset in A), showing the position of the cores in relation to the OMZ and the principal subsurface currents. Triangles are our box cores (B-13, B-6 and B-14) and black circles are other sites considered in the discussion: Soledad and Pescadero (Deutsch et al., 2014) and ST-18 (Srain et al., 2015).

Dissolved Oxygen Concentration

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.

Figure 2. Time series (1960 – 2010) of dissolved oxygen off central Peru. (A) Map showing the IMARPE oceanographical stations between 30 and 100 NM and 11 and 15° S. (B) Mean dissolved oxygen profile with error values in the water column for the area depicted in (A) (data available in World Ocean Database 2018). Time series of dissolved oxygen at 60 m (C), 150 m (D) and 200 m (E) depth. Red areas indicate error values. Trend lines with confidence intervals (orange) for each depth are indicated.

Core Sampling

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 (Figures 1A,C) while the Callao box core (B-13) is closer to the upper limit of the OMZ.

Benthic Foraminifera

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.

Geochemical Proxies

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, ${\mathrm{HNO}}_3^{-}$ and ${\mathrm{HClO}}_4^{-}$). 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.

Statistical Analyses

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.

Results

Dissolved Oxygen Off Central Peru in 1960–2010

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 (Figures 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 (Figures 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.

Benthic Foraminifera in the Box Cores (∼1835–2005 Period)

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.

Figure 3. Records of foraminifera in in Callao (B-13, black circles) and Pisco (B-6, white circles) box cores: (A) Total densities; (B) P/B ratio and (C–M) relative abundances (%) of main benthic foraminiferal species.

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.

Table 1. Diversity indices of benthic foraminifera for box cores B-13 and B-6.

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, Figures 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 (Figures 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).

Figure 4. Interpretation of benthic foraminiferal assemblages. (A) Distribution of living species in central Peru from several databases: Callao and Pisco upper margin for the 2009 – 2011 period (Cardich et al., 2015) and for the 2012 – 2017 period (Benthos Lab of IMARPE); Callao and Huacho upper margin (M77-2 cruise, 2008; Mallon et al., 2012); Callao upper slope (Panorama Leg3a cruise, 1997; Pérez et al., 2002) and (M92-2 cruise, 2013; Cardich et al., un published). Width of bars indicate mean relative abundance for each species. (B,C) Loadings of principal component 1 (black bars) and principal component 2 (red bars) for B-13 (Callao) and B-6 (Pisco). Opposite values of B-6’s PC2 scores are shown to maintain coherence in the axes.

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 (Figures 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 (Figures 4B,C). PC1 resembles the negative covariation between B. seminuda and B. tenuata as observed in Figures 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 (Figures 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.

Principal Components of Geochemical and Biogenic Proxies in the Box Cores (∼1835–2005 Period)

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).

Figure 5. Principal component analysis based on both geochemical and foraminiferal matrices for B-13 (A) (Callao; left) and B-6 (B) (Pisco; right). Geochemical proxies are: Total organic carbon (TOC), δ¹⁵N, biogenic silica (Si), Mo/U, Mo enrichment factor (Mo EF). Foraminiferal proxies are the same benthic species in Figure 4, plus the planktonic to benthic ratio (P/B).

Discussion

Understanding the Benthic Foraminiferal Assemblages in the Cores

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 (Figures 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 (Figures 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.

Table 2. Spearman correlation matrix for biogenic and geochemical proxies from B-13 and B-06.

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 (Figures 6E,L) to the δ¹⁵N records (Figures 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.

Figure 6. Trends of geochemical and foraminiferal proxies in Callao (left panel) and Pisco (right panel) along to ENSO index (EN 1+2; A,H): δ¹⁵N anomalies (B,I), as a proxy of water-column denitrification; Mo EF (C,J) as proxy of intensity of anoxia; Mo/U ratio (D,K) as proxy of suboxic (low values) vs. sulfidic (high values) conditions; Bt/(Bt+Bs) index (E,L); Asp/(Asp+Psp) index (F,M); and PC scores from Figure 5 as indicators of redox conditions [RC1 of B-13 (G) and RC2 of B-6 (N)].

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.

(Multi)Decadal Sedimentary Patterns and Response to Climatic Variability

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:

(i) From ca. 1835 to ca. 1875, with high export productivity (large peaks of biogenic silica), intense water column denitrification (high values of δ¹⁵N) and sulfidic sediments (high Mo EF and low Mo/U values);

(ii) From ca. 1875 to ca. 1960, with increasing export productivity, a decreasing trend of water column denitrification and a relaxation of reducing sedimentary conditions; and

(iii) From ca. 1960 to ca. 2005, the increasing trend of export productivity continues, but with local differences of reducing sedimentary conditions and water column denitrification.

Benthic foraminiferal ratios (Figures 6E,F) and rotated principal components (Figures 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.

Figure 7. Ocean circulation indices and records for (de)oxygenation in the East Pacific upwelling systems: (A) Equatorial Undercurrent (EUC) zonal velocity (Drenkard and Karnauskas, 2014); (B) mean anomalies of δ¹⁵N in two locations off California (black dots in Figure 2B; Deutsch et al., 2014); (C) δ¹⁵N record in B-14 (gray circles, this study) and instrumental dissolved oxygen anomalies at 200 m depth off central Peru (red line, inverse axis, this study); (D) Psp/(Psp+Asp) index anomalies (this study); and (E) Mo excess content anomalies off Concepción (Srain et al., 2015).

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.

Regional Trend of Oxygenation in the Eastern Pacific

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 (Figures 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 (Figures 7C,D).

The slight oxygenation trend near the equator (central Peru) as observed in the oxygen data and the paleo data (Figures 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.

Conclusion

We show evidence for a subtle oxygenation trend from 1960 to 2010 in the core of the OMZ off central Peru. This trend is part of the OMZ variability at (multi)decadal to centennial time scales since the second half of the 19th century evidenced by a multiproxy approach. Near the coast, an intensification of the subsurface deoxygenation during the mid-nineteenth century resulted from the interplay of local and regional factors. Decadal to multidecadal variability of redox conditions might be controlled by local productivity changes affecting the OMZ intensity. At a larger scale, subsurface ventilation through the EUC and/or CPDCC is the ultimate factor controlling OMZ variability at multidecadal to centennial scale off Peru and probably in the Tropical South Eastern Pacific region. To be fully understood, the mechanisms behind the subsurface (de)oxygenation trends in Eastern Boundary Upwelling Ecosystems during past periods need the use of subsurface circulation proxies together with earth models at regional scale in addition to geochemical/biogenic proxies.

Author Contributions

JC and DG designed the manuscript and landed the original idea of the manuscript. AS, RS, and FB-Z participated in data interpretation and discussion since the conception of the idea. RS analyzed the geochemical proxies. DR and CA analyzed the benthic foraminifera from the Pisco box core. JC analyzed benthic foraminifera in the Callao box core. MG and TA analyzed the oceanographic data from Peru and in its interpretation. All authors participated in the general discussion of the manuscript.

Funding

This project was funded by the MAGNET Program of CONCYTEC N°007-2017. This work was supported by the International Joint Laboratory “PALEOTRACES” (IRD, France; UPMC, France; UFF, Brazil; UA, Chile; UPCH, Peru), the Department of Geochemistry of the Universidade Federal Fluminense (UFF, Brazil), and the Peruvian Marine Research Institute (IMARPE). JC acknowledges the financial support from CNPq (Grant 142001/2013-9).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We give special thanks to Ioanna Bouloubassi and Myriam Khodri for their comments and suggestions. We also thank the crews aboard the RVs ‘SNP 2’ and ‘José Olaya Balandra’ during the sampling cruises.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2019.00270/full#supplementary-material

References