Oceanic nutrient rise and the late Miocene inception of Pacific oxygen-deficient zones

The oxygen (O₂) content of the ocean is critical for the health of marine ecosystems and productive fisheries. Observations show that the oceans have been losing O₂ over the past five decades in both open ocean and coastal zones, attributed in part to anthropogenic climate change and coastal nutrient discharge (1, 2). However, simulating the lowest O₂ regions of the ocean, the oxygen-deficient zones (ODZs), is a challenge for climate models, undercutting the predictions of future ocean deoxygenation (3).

In general, the O₂ content in the ocean interior is determined by the balance between O₂ supply through ocean circulation and O₂ consumption by the respiration of sinking organic matter imported from the surface ocean. In the ODZs, ocean circulation patterns minimize O₂ supply (4), and wind-driven upwelling drives high productivity in overlying surface waters and thus a rain of sinking organic matter, causing O₂ concentrations to fall below ∼5 μmol/L, at which point nitrate is used for microbial respiration (5). Three major ODZs are found in the modern ocean, two of which are located in the eastern Pacific Ocean (the Eastern Tropical South Pacific [ETSP] and the Eastern Tropical North Pacific Ocean) at depths of roughly 150 to 500 m (SI Appendix, Fig. S1). Observations have demonstrated the expansion of the Pacific ODZs over the past several decades (6). However, it remains unclear whether they will continue to expand under global warming, because the O₂ supply to the depths of the ODZs as well as the productivity in overlying surface waters may be affected by regional changes in atmospheric and oceanic conditions (3, 7). For example, some studies have suggested that, on longer timescales (i.e., more than a century), the ODZs may shrink after the initial expansion (8).

The history of the ODZs may provide important insights into their sensitivities to global warming and other human perturbations. However, reconstructions of variations in Pacific ODZs have thus far been restricted mainly to the Pleistocene (9, 10), with little known on longer timescales such as during past warm periods (e.g., the Pliocene, the Miocene). We measured foraminifera-bound δ¹⁵N (where δ¹⁵N = [(¹⁵N/¹⁴N_sample)/(¹⁵N/¹⁴N_air) − 1] × 1,000) to reconstruct a 12-Mya history of the South Pacific ODZs, which reveals a large, gradual expansion since the late Miocene. We also observe a concomitant rise in phosphate concentrations ([PO₄³⁻]) in the South Pacific reconstructed over the same interval, using the phosphorus (P) and iron (Fe) content of hydrothermal sediments. Taken together, these data illuminate nutrient dynamics and their impacts on the O₂ content of the ocean, the global carbon cycle, and climate change since the late Miocene.

Our study sites Deep Sea Drilling Project (DSDP) 598 (19.00 ^oS, 124.68 ^oW, 3,699 m) and OC-73–3-20 (19.25 °S, 113.58 °W, 3,081 m) are on the Southern Eastern Pacific Rise (SEPR; Fig. 1), a midocean ridge system with fast spreading rate and high hydrothermal activity (11). Today, OC-73–3-20 lies ∼8 km west of the SEPR axis, with DSDP 598 located ∼1,150 km to the west of OC-73–3-20. With an estimated spreading rate of 77 mm/y for the SEPR (12), DSDP 598 would have been only ∼200 km away from the SEPR axis 12 Mya. In this study, the 12-My records were generated from DSDP 598. Core OC-73–3-20, which is too short for long-term reconstructions, was used as a late Pleistocene comparison site to validate the proxies and correct for paleodepth variations.

These SEPR sites are ideal for the reconstruction of the history of Pacific ODZs for two reasons. First, these sites are close to the ETSP and the thermocline is directly affected by the ODZs through the transport of the South Equatorial Current, resulting in a high-nitrate δ¹⁵N of 12 to 13‰ from the water column denitrification in the ODZs (SI Appendix, Fig. S1) (13). Second, these sites are located within the oligotrophic zone of the South Pacific gyre, with complete nitrate consumption by phytoplankton in the surface ocean. The foraminifera-bound δ¹⁵N records from the SEPR sites, unlike the records from within the ODZs, are not complicated by fractionations due to incomplete nitrate consumption. Previous studies suggest that in ocean regions with complete nitrate consumption, foraminifera-bound δ¹⁵N is a proxy for the thermocline nitrate δ¹⁵N, because the thermocline nitrate is the dominant source of N to the euphotic zone, even in regions with relatively high rates of N₂ fixation (14). Hence, it is expected that the foraminifera-bound δ¹⁵N from the SEPR provides a direct measure of its thermocline nitrate δ¹⁵N and thus the history of ODZs in the ETSP. Indeed, we observed very high Holocene-age δ¹⁵N values (14 to 16‰) in OC-73–3-20 (SI Appendix, Fig. S2), consistent with the high δ¹⁵N values seen in the modern thermocline nitrate. The Pleistocene foraminifera-bound δ¹⁵N in DSDP 598 matches the Holocene δ¹⁵N values observed in OC-73–3-20, supporting the use of DSDP 598 foraminifera-bound δ¹⁵N to reconstruct the long-term history of the ODZs in the ETSP.

In addition to the foraminifera-bound δ¹⁵N record, the SEPR sites provide a unique opportunity to reconstruct the nutrient content in the deep Pacific Ocean using the P and Fe content in these hydrothermal sediments. Phosphate removal by Fe oxyhydroxides in midocean ridges has long been recognized as a sink of P from seawater (15). The fast adsorption/scavenging of phosphate by abundant Fe oxyhydroxides is limited by the available phosphate, such that the Fe:P ratio in hydrothermal sediments is mainly influenced by the ambient deep-ocean [PO₄³⁻]. It has been observed that the Fe:P ratio in hydrothermal precipitates displays a tight correlation with the in situ [PO₄³⁻] (16). Furthermore, the hydrothermal phosphate sink has a minimal impact on the local phosphate concentration (i.e., ≤0.1 μM) (17), such that the processes that generate the Fe-P tracer do not self-limit the P scavenging process. These mechanics provided an approach to use measurements of the Fe and P contents of hydrothermal sediments as a proxy for seawater [PO₄³⁻] in the past.

Although the mechanics of the Fe-bound P proxy are well understood, this proxy has been little used in paleoceanographic reconstructions, mainly because of a lack of suitable sedimentary materials. The SEPR sites provide a unique opportunity to apply this proxy because (1) the sites are adjacent to an active hydrothermal vent with high Fe oxyhydroxides production; (2); Fe-bound P is the dominant P species in the sediments; (3) the sites are in an oligotrophic setting with very low organic matter content, minimizing postdepositional alterations (SI Appendix, SI Text); and (4) the sediment cores used in this study have good age models with no stratigraphic disturbance. To further validate the Fe-bound P proxy at the SEPR sites, we analyzed the Fe and P content in OC-73–3-20 spanning the past ∼40 ky. Independent reconstructions of deep-ocean phosphate concentrations from foraminifera Cd/Ca data suggest little change on this timescale (18, 19). Indeed, the Fe and P content display a strong correlation in OC-73–3-20 with only one slope, consistent with the Cd/Ca reconstructions and supporting the use of the Fe-bound P proxy on longer timescales.

In DSDP 598, the Fe and P contents display tight correlations during each epoch/stage (SI Appendix, Fig. S3, r² > 0.97). However, the slopes of this Fe-P correlation decreased dramatically over the past 10 My—from ∼94 in the Tortonian to ∼31 in the Messinian and ∼19 in the Pleistocene (SI Appendix, Fig. S3). The Pleistocene samples in DSDP 598 display a Fe-P slope only a little higher than that of nearby core OC-73–3-20 (19.4 ± 1.76 versus 14.87 ± 1.52), supporting the robustness of this proxy for deep-water phosphate concentrations. To quantify the deep-ocean [PO₄³⁻] for the past 12 My, we computed the running slopes of the Fe-P records for a window size of 1.5 My and converted them to the deep-ocean [PO₄³⁻] by calibrating our late Pleistocene observations with the modern [PO₄³⁻] value at the study site (see SI Appendix, Fig. S6 and SI Text). Over the entire interval, the correlation coefficients of Fe-P slopes remain higher than 0.85 (SI Appendix, Fig. S6), consistent with the proximity of this site to the ridge and signifying that hydrothermal processes have persistently been the dominant control on the Fe and P content of the sediments captured by this core. We applied a depth correction to the reconstructed [PO₄³⁻], because the DSDP 598 site has been migrating to deeper water as it moves further off axis. However, this correction causes negligible changes to the results (i.e., <0.05 μM) because [PO₄³⁻] is largely homogeneous below 3,000 m in the Pacific Ocean (Fig. 1).

In summary, we generated records of foraminifera-bound δ¹⁵N and deep-ocean [PO₄³⁻] in the South Pacific since 12 My (Fig. 2). The most prominent feature of the foraminifera-bound δ¹⁵N record is a gradual and large increase of 10 to 11‰ since the late Miocene, from average δ¹⁵N values of ∼6‰ before 9 My to 16 to 17‰ by 1 My. With an offset of 2 to 3‰ between Holocene foraminifera-bound δ¹⁵N and subsurface nitrate δ¹⁵N (SI Appendix, Fig. S2 and SI Text), the subsurface nitrate δ¹⁵N in the eastern Pacific before 9 My would have been 3 to 4‰, indicating much smaller ODZs at that time. Thus, the late Miocene marks a dramatic intensification of South Pacific ODZs. In parallel, the deep-ocean [PO₄³⁻] record displays a large increase since the late Miocene. Before 8 My, the [PO₄³⁻] in the deep Pacific was only ∼0.5 μM—a value that is less than 50% of the minimum [PO₄³⁻] in the modern open ocean. From 8 to 6 My, the deep-Pacific [PO₄³⁻] increased by 0.8 μM, representing the largest increase over the past 12 My. This increase is corroborated by a benthic foraminifera Cd/Ca record from the west Pacific (SI Appendix, Fig. S7) and coincides with similar changes in global benthic foraminifera δ¹³C records, as well as the expansion of C4 plants on land exemplified by the soil carbonate δ¹³C record from northern Pakistan (Fig. 2). The δ¹⁵N change appears to have preceded the increase in [PO₄³⁻] change by ∼1–2 My (Fig. 2). Because the δ¹⁵N and [PO₄³⁻] records were generated on the same samples, this 1–2 My lag is not due to age model uncertainties. Rather, the timing of the δ¹⁵N change is more consistent with the sea-surface temperature (SST) records (Fig. 2), which showed that cooling may have preceded the δ¹³C and [PO₄³⁻] changes. It is thus likely that the early change in the δ¹⁵N record before 8 My was related to global cooling, such as through high-latitude ocean stratification and its impact on deep-ocean O₂ concentrations (20). However, the effect of global cooling on the expansion of ODZs since the late Miocene is only secondary, as indicated by the model results below. Nevertheless, most of the δ¹⁵N increase occurred in concert with the [PO₄³⁻] increase.

Although our study site is well within the oligotrophic waters of the South Pacific today (Fig. 1), one may wonder whether the equatorial Pacific high-nutrient tongue extended further south in the past, causing incomplete surface nitrate consumption and lowering foraminifera-bound δ¹⁵N at our study sites in the early part of the record. However, the effect of the high-nutrient tongue would be to increase, not decrease, foraminifera-bound δ¹⁵N at the latitudes of our core sites. This is the case because the residual nitrate and particulate nitrogen at the margins of the high-nutrient tongue are elevated in δ¹⁵N due to prior nitrate consumption (with isotopic fractionation) closer to the equator, coupled with flow away from the equator (13, 21, 22). Thus, the observed low foraminifera-bound δ¹⁵N values before 8 My require a smaller ODZs and an equatorial nutrient tongue that was either similar to or smaller than that of today. Indeed, the lower whole ocean [PO₄³⁻] that we infer for the late Miocene may well lead to a contraction in the equatorial Pacific nutrient tongue. Regardless of past changes in the equatorial nutrient tongue, the foraminifera-bound δ¹⁵N values of ≤8‰ before 9 My could not have been reached without a smaller ODZs at that time.

Our deep-ocean [PO₄³⁻] record was derived from the South Pacific, and while this reflects the largest body of water on the planet, it is important to consider whether this record may be caused by a redistribution of phosphate between the Atlantic and Pacific without requiring a change in the global ocean’s nutrient content. It has long been thought that the formation of the Isthmus of Panama would have changed ocean circulation (23–25) and affected the distribution of nutrients between the Atlantic and Pacific Oceans (26). Before the closure of the Central American Seaway (CAS), a significant amount of North Atlantic deep water could flow through the CAS to the Pacific, which has been shown in models to decrease the nutrient content of the deep Pacific (25, 26). Although the timing of the formation of the Isthmus of Panama is debated, most paleoceanographic data point to a timing more recent than 5 My (27–29). Furthermore, the reconstructed [PO₄³⁻] before 8 My was only 0.5 μM, less than 50% of minimum [PO₄³⁻] in the modern deep ocean. A redistribution of nutrients between the Atlantic and Pacific Oceans would not be adequate to cause such low [PO₄³⁻]. Consequently, the increase in deep Pacific [PO₄³⁻] during the late Miocene is best explained as an increase in the mean [PO₄³⁻] of the global ocean. This interpretation is supported by the δ¹³C gradient reconstructed between the Atlantic and Pacific Oceans, which showed minimal change during this time (30).

We used an idealized two-box model to quantify the extent of suboxia and its relationship with the rising nutrient content of the ocean since the late Miocene (SI Appendix, Fig. S8 and SI Text). In this model, the combined ODZs (<1% of total ocean volume) are separated from the rest of the ocean (mean ocean). Fixed nitrogen (i.e., “biologically available nitrogen”) is lost from the ODZs through water-column denitrification and from the mean ocean through sedimentary denitrification. Fixed nitrogen is added to the ocean through N₂ fixation at a rate equivalent to that of total denitrification. Assuming a constant sedimentary denitrification rate of 130 Tg N/y (31) and a constant Redfield ratio (mean ocean [NO₃⁻]/[PO₄³⁻] = 16), we derive the relationships between ODZs nitrate δ¹⁵N, water column denitrification rate, and mean ocean nutrient content (contours in Fig. 3). The reconstructed [PO₄³⁻] and δ¹⁵N records are used to calculate the rates of water-column denitrification (suboxia) from the late Miocene to the present (symbols and arrows in Fig. 3).

From 9 to 8 My, the mean ocean nutrient content did not change, while ODZs δ¹⁵N increased by 5 to 6‰. However, the model indicates that this 5 to 6‰ increase reflected only a modest increase in water-column denitrification rate (∼10 Tg N/y, or ∼17% of modern value). This is mainly because of the low nutrient content and near absence of water-column denitrification before 9 My. In such an ocean, the ODZs δ¹⁵N is highly sensitive to small changes in water-column denitrification rates (Fig. 3). Since 8 My, both water column denitrification and mean ocean nutrient content gradually increased to their modern values, with the largest change between 8 and 6 My. The evolution of the water-column denitrification rate and mean ocean nutrient content were linearly well correlated (r² > 0.98) since 8 My (Fig. 3), with a slope of ∼31 Tg N/y/μM [PO₄³⁻]. This sensitivity of water-column denitrification rates to mean ocean nutrient change falls within the range of estimates from more complex biogeochemical models (7, 32). The strong linear correlation between water-column denitrification rate and mean ocean nutrient content since 8 My suggests that most of the ODZs expansion was driven by the nutrient increase. Indeed, the nutrient content of the ocean plays a dominant role in the biological pump, the strength of which is inversely related to the O₂ concentration of the ocean interior through export production and microbial respiration (33, 34). Consistent with this expectation, an increase in productivity metrics have been observed since the late Miocene, with a “biogenic bloom” between 8 and 6 My (35–37). Although the decreasing SST and increasing O₂ production may have supplied more O₂ to the ocean since the late Miocene (SI Appendix, SI Text), the rising oceanic nutrient content and resulting higher productivity appear to have overwhelmed these effects, driving ocean deoxygenation over the past 8 My. Whereas studies of global warming have tended to cast ocean O₂ content as dependent on climate and ocean circulation (1, 2), our study indicates that the modern ODZs are underpinned by a historically high oceanic concentration of phosphate.

The rise in the water-column denitrification rate since the late Miocene has important implications for the rate of N₂ fixation in the global ocean. As the residence time of fixed nitrogen is only 3 to 5 ky in the ocean (38), the rate of N₂ fixation is set equal to the rate of total denitrification in the model such that fixed N losses and inputs are always balanced on million-year timescales. As a result, the N₂ fixation rate has increased since the late Miocene in the model, in response to the increasing watercolumn denitrification rate. In today’s ocean, the highest N₂ fixation rates are observed in the western Pacific, with little N₂ fixation near the Pacific ODZs (39). Our reconstructed ODZs history would imply that the N₂ fixation rate in the western Pacific before 8 My was lower than today.

According to our data, the late Miocene experienced the largest jump in ocean [PO₄³⁻]. The late Miocene interval marks a time of substantial change in both the climate and biosphere. On land, plants using the C4 photosynthesis pathway—a metabolic strategy for concentrating CO₂ in plant tissues—expanded, changing landscapes and ecosystems in subtropical and tropical regions (40). Africa and East Asia became more arid, and large deserts such as the Sahara emerged (41). These vegetation and landscape shifts also coincided with the divergence between lineages that would give rise to humans and chimpanzees (42). In the oceans, SST records from both hemispheres show a dramatic decline from 7 to 4 My, with a concomitant increase in the meridional temperature gradient (37). These environmental changes have long been linked to evidence for an increase in ocean productivity (36) and a strengthening of the biological carbon pump of the oceans (30).

One striking feature of the [PO₄³⁻] record is that the late Miocene jump coincided with the global benthic δ¹³C decline and the expansion of C4 plants on land, implying a common driver of these changes (Fig. 2). Although recent studies suggested that the timing of the expansion of C4 plants may have been different at different places, most of these changes occurred between 8 and 6 My (43). One hypothesis proposed for the late Miocene benthic δ¹³C decline is that the expansion of C4 plants changed the land surfaces to less dense grasslands, triggering greater soil erosion and higher soil organic matter input to the oceans (36). The subsequent remineralization of low-δ¹³C soil organic matter in the ocean would have caused a decline in the δ¹³C of the mean ocean carbon reservoir (37). Because soil organic matter contains a significant amount of P, the remineralization of soil organic carbon would have also increased the P input into the ocean, causing a transient imbalance between the source and sinks of oceanic P reservoir. The net result would be a rise in the mean ocean’s nutrient content. Our reconstructed deep-ocean [PO₄³⁻] rose by ∼0.8 μM from 8 to 6 My, consistent with the soil erosion hypothesis.

The deep-ocean [PO₄³⁻] continued to rise after 6 My, however, while the global benthic and soil δ¹³C records stabilized at 6 My (Fig. 2). Also, the soil organic matter reservoir is a temporary carbon reservoir on geological timescales, so faster soil erosion alone cannot sustain the accelerated delivery of nutrients to the ocean. It is thus unlikely that the nutrient increase after 6 My was caused by the same soil erosion mechanism due to C4 expansion. We propose that the increase in the deep-Pacific [PO₄³⁻] after 6 My was caused by an increased weathering rate. Recent work suggests that the global weatherability increased since the mid- to late Miocene due to the emergence of the Southeast Asian islands as a result of arc-continent collision (44). If correct, then such an increase in weathering could have contributed to the observed increase in the deep ocean [PO₄³⁻], particularly the gradual deep-ocean [PO₄³⁻] increase after 6 My.

Beyond the expansion of marine ODZs, the increase in the nutrient content of the oceans has important implications for the atmospheric pCO₂. The extra nutrients would have increased primary productivity in much of the world’s oceans, sequestering additional CO₂ from the atmosphere in the deep ocean. This ocean nutrient content mechanism was originally proposed to explain the pCO₂ variation seen on glacial-interglacial timescales (33). Although no convincing evidence has emerged for such a change in glacial-interglacial timescales, our results suggest that this mechanism was at work earlier in the history of the Earth. The ∼0.8-μM [PO₄³⁻] increase in the late Miocene could have caused an abrupt ∼60 ppm decrease in atmospheric pCO₂, depending on the responses of pelagic CaCO₃ production and high latitude productivity (45). Since 6 My, increased weathering of P and of alkalinity (46) may have worked together to further reduce atmospheric pCO₂. Although Miocene-age pCO₂ reconstructions have yielded mixed results (SI Appendix, Fig. S7), some studies have shown a pCO₂ decline during this interval (47, 48). If confirmed, such a decrease in atmospheric CO₂ would help to explain the observed global cooling and other environmental and evolutionary changes since the late Miocene, such as the increase in aridity in Asia and Africa and the evolution and expansion of C4 photosynthesis.