Climate-driven Oxygen Minimum Zone (OMZ) expansions in the geologic record provide an opportunity to characterize the spatial and temporal scales of OMZ change. Here we investigate OMZ expansion through the global-scale warming event of the most recent deglaciation (18-11 ka), an event with clear relevance to understanding modern anthropogenic climate change. Deglacial marine sediment records were compiled to quantify the vertical extent, intensity, surface area and volume impingements of hypoxic waters upon continental margins. By integrating sediment records (183-2,309 meters below sea level; mbsl) containing one or more geochemical, sedimentary or microfossil oxygenation proxies integrated with analyses of eustatic sea level rise, we reconstruct the timing, depth and intensity of seafloor hypoxia. The maximum vertical OMZ extent during the deglaciation was variable by region: Subarctic Pacific (~600-2,900 mbsl), California Current (~330-1,500 mbsl), Mexico Margin (~330-830 mbsl), and the Humboldt Current and Equatorial Pacific (~110-3,100 mbsl). The timing of OMZ expansion is regionally coherent but not globally synchronous. Subarctic Pacific and California Current continental margins exhibit tight correlation to the oscillations of Northern Hemisphere deglacial events (Termination IA, Bølling-Allerød, Younger Dryas and Termination IB). Southern regions (Mexico Margin and the Equatorial Pacific and Humboldt Current) exhibit hypoxia expansion prior to Termination IA (~14.7 ka), and no regional oxygenation oscillations. Our analyses provide new evidence for the geographically and vertically extensive expansion of OMZs, and the extreme compression of upper-ocean oxygenated ecosystems during the geologically recent deglaciation.

Funding: Support was provided by the National Science Foundation (OCE 0825322 and OCE 1255194 to TMH), the University of California Multicampus Research Programs and Initiatives (to TMH), University of California Davis REACH IGERT (to SEM), Mia Tegner Historical Ecology Grant (to SEM), the EPA STAR Fellowship (to SEM) and Switzer Environmental Fellowship (to SEM). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Copyright: © 2015 Moffitt et al. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited

Geochemistry . Isotopic and trace element geochemical proxies of oxygenation are complex, as the cycling, preservation and attribution of these elements can be controlled by processes such as of surface productivity and flux, water column processes, sediment-water interface flux and diffusion, and sub-surface sedimentary processes. Denitrification in continental margin regions occurs under hypoxic conditions in both the water column and sediment ([O 2 ]<0.23 ml L -1 ) [ 19 , 20 ], and changes in the nitrogen isotopic signal are a product of regional changes in the biological available N pool [ 35 ]. Although intense denitrification is often indicative of severe hypoxia, more subtle changes in δ 15 N can be attributed to a range of related localized factors, such as source water and nutrient availability [ 19 , 36 – 39 ]. Redox sensitive trace elements such as cadmium (Cd), uranium (U), chromium (Cr) and rhenium (Re) occur in bulk sediment and carbonate fossils, and can be used to characterize the redox state of the seafloor as well as functioning as proxies of paleoenvironmental redox changes [ 40 – 43 ]. Paleoproductivity proxies are useful supportive proxies in OMZ reconstructions, due to the mechanistic coupling between surface export and subsurface respiration [ 44 ]. Additionally, the δ 13 C record of planktonic carbonate reflects the surface ocean isotopic pool, and can be used to reconstruct surface water productivity and carbon export [ 45 – 47 ].

Microfossils . Oxygen is a physiochemical parameter that acts as a limiting factor, thereby determining the distribution of where organisms can live based on their physiological requirements. Within the optimal range of a limiting factor, a species will reach maximum competitiveness and maximum abundance [ 30 ], resulting in tight faunal affiliations with specific oxygenation ranges. Benthic foraminifera have species-specific oxygenation thresholds and therefore marker taxa function as oxygenation proxies [ 31 ]. Metrics of foraminiferan community structure, such as diversity and density, also record changes in seafloor oxygenation [ 22 , 23 , 32 ]. Due to their opportunistic responsiveness to environmental change, Foraminifera are ideally suited for high-resolution oxygenation reconstructions [ 32 – 34 ].

Seafloor hypoxia proxies for paleoceanographic reconstructions, partitioned by the thresholds and capacity each proxy has to record fine-scale changes in seafloor hypoxia, as well as organic flux to the seafloor [ 24 ].

Sedimentary structures . Severe hypoxic conditions preclude benthic macrofauna, preventing bioturbation and thereby allowing for the preservation of laminated sediments ( Table 2 ) [ 25 – 27 ]. For example, annual sediment laminations (i.e. unmixed, fine-grained sediment displaying distinct, continuous layering) are formed in modern California margin basinal features at [O 2 ]<0.1 ml L -1 [ 25 , 28 ]. Laminations are formed under the absence of bioturbating invertebrates [ 29 ] and sufficiently high organic carbon export from the surface [ 16 ], and are one of the clearest indicators of severe hypoxia in benthic environments.

These hypoxia categories are detailed in Table 1 , and follow Hofmann et al., [ 24 ]. Hypoxia proxies include [Re], [Mn], [U], [Cd], [Mo], δ 15 N, foraminiferan communities, and sedimentary laminations. Units for each proxy reflect the cited literature, which constrains the proxy to a specific oxygenation category.

We employ a multi-proxy approach to oxygenation reconstructions here, and include core data across sedimentary, faunal and geochemical proxies ( Table 1 ). Fig. 3 depicts, in schemata form, how regional multi-proxy oxygenation data are interpreted, and we discuss each proxy below.

Salinity (34 psu) and hydrostatic pressure (P = 10 bar) are assumed constant. Temperature columns indicate the temperature used for partial pressure and saturation calculations at the associated concentrations. Data table adapted from Hofmann et al., [ 24 ]. The intermediate category is described as “coastal” in Hofmann et al., [ 24 ]. We refrain from using this term here, to prevent confusion between hypoxia categories and offshore habitat locations

For this work, we follow the hypoxia thresholds and categories defined in [ 24 ], which synthesizes the existing hypoxia vernacular, to draw thresholds that are biologically meaningful ( Table 1 ). Mild hypoxia begins at [O 2 ]<2.45 ml L -1 and is the threshold where sensitive species exhibit avoidance reactions. Intermediate hypoxia, often referred to as “coastal hypoxia”, occurs at [O 2 ]<1.4 ml L -1 and is the threshold wherein ecosystems are dominated by organisms with adaptive features. Severe hypoxia ([O 2 ]<0.5 ml L -1 ) is a threshold at which mass mortality is induced for most organisms, past which only highly specialized species can survive [ 24 ].

The geospatial distributions of severely hypoxic [O 2 ] minimums (of [O 2 ] = 0.5 ml L -1 and [O 2 ] = 0.2 ml L -1 ) are depicted on both panels as white lines. For the upper panel, regional blocks are defined by black lines to highlight where paleoxygenation reconstructions were completed. Data from World Ocean Atlas [ 192 ].

OMZs are tightly coupled to upwelling systems and Eastern Boundary Currents, such as the California Current, the Humboldt Current and the Benguela Current, as well as the Oman and Pakistan Margin in the Indian Ocean ( Fig. 2 ). In these regions, respiration within the pycnocline depletes dissolved oxygen and simultaneously enriches seawater in the carbon and nitrogen byproducts of respiration [ 18 ]. Marine denitrification occurs within OMZs (e.g., [ 19 , 20 ]) therefore the physical extent and intensity of OMZs is inherently coupled to the oceanic nitrogen cycle. OMZs form at shelf and upper slope depths, and are considered to be unique biological, geochemical and evolutionary environments, analogous to cold seep or deep-sea vent environments [ 21 ]. As continental margin ecosystems transition from well oxygenated surface waters to the hypoxic core of the OMZ ([O 2 ] = 0.5–0.1 ml L -1 ), faunal diversity, trophic structures and physiological strategies change (e.g., [ 22 , 23 ]). OMZ oxygenation gradients produce successional biological zonation and are fundamental habitat barriers for benthic and pelagic organisms [ 21 ].

Here we synthesize published continental margin sediment core records to investigate Oxygen Minimum Zone (OMZ) changes through the last deglaciation. We build on previous syntheses of oxygenation proxy records (e.g., [ 7 , 16 , 17 ]), and provide a focus on regional-scale sensitivity. By integrating sediment records, sea level change, and high-resolution bathymetry, we provide geospatially analyzed paleoceanographic data that are interpretive baselines for modern oceanography and global environmental change.

Hypoxia substantially degrades ecosystems through mass mortality events, the alteration of food-web structures and the loss of habitat (e.g., [ 11 ]). Changes in [O 2 ] in the ocean interior have broad consequence for global biodiversity, marine economic resources and ocean management [ 12 ]. To grasp the scale of future hypoxia disturbance, the paleoceanographic (pre-instrumental) record of natural variability provides a critical analytical and interpretive window. Ocean sediments are climatic and environmental archives, which preserve geochemical, microfaunal and sedimentary evidence that record globally relevant Earth system events, similar to other climate records such as the Greenland Ice Sheet Project [ 13 ]. Recent investigations reveal that paleoceanographic investigations hold valuable insight into modern environmental conservation and management [ 14 , 15 ].

Glacial Termination IA (14.7 ka) is an event of rapid warming in the Northern Hemisphere, which initiates the warm interval of the Bølling-Allerød (B/A) from 14.7–12.9 ka. The Younger Dryas (YD), a reversal towards cool conditions from 12.9–11.7 ka, follows the B/A. The YD ends with glacial Termination IB (11.7 ka), a subsequent rapid warming event. Deglacial warming in the Southern Hemisphere begins at 18 ka.

Resolving records of global change through the most recent deglaciation event (18–11 ka) is one of the primary challenges to developing cohesive and robust theories regarding rapid climate change [ 1 ]. The last deglaciation was a profound event in the global climate system, wherein atmospheric [CO 2 ] increased by 80–100 ppmv [ 2 , 3 ], global average temperature rose 3–4°C [ 4 ], and sea levels rose ~110 m ( Fig. 1 ) [ 5 , 6 ]. This change was also accompanied by the pervasive loss of dissolved oxygen in the upper ocean [ 7 ], with unknown impacts on large marine ecosystems. The coupling between climate, carbon emissions and subsurface dissolved oxygen is a salient and critical element of anthropogenic climate change. The global inventory of ocean oxygen is predicted to decline between 1 and 7% by the year 2100, through stratification, ventilation reduction and decreased O 2 solubility (e.g., [ 8 ]), and the hypoxic water volume in the global ocean is predicted to increase by 50% [ 9 ]. Climate model results beyond the 100-year window reveal extensive oceanic deoxygenation, on thousand-year timescales, under “business-as-usual” carbon emission scenarios, and show that oceanic deoxygenation is a fundamental and long-lasting property of anthropogenic carbon perturbation (e.g. [ 10 ]).

Continental margin volume is corrected for eustatic sea level rise [ 5 , 6 ] and calculated using the SRTM30_PLUS global bathymetry dataset [ 50 ]. Hypoxic and oxic margin volume (km 3 ) and percent (%) are calculated for the water column above a region-specific isobath: SP (0–3,300 mbsl), CC (0–2,400 mbsl), BM (1–1,200 mbsl), and HC (1–3,300 mbsl). Upper ocean oxic volume (km 3 ), from the surface ocean to the upper subsurface hypoxic boundary, is included.

Continental margin surface area is corrected for eustatic sea level rise [ 5 , 6 ] and calculated using the SRTM30_PLUS global bathymetry dataset [ 50 ]. Total vertical extent of hypoxia is stated, including intermediate and severe. Hypoxic and oxic margin surface area (km 2 ) and percent (%) are calculated for the seafloor above a region-specific isobaths: SP (0–3,300 mbsl), CC (0–2,400 mbsl), MM (1–1,200 mbsl), and HC (0–3,300 mbsl). Upper ocean oxic surface area (km 2 ), from the surface ocean to the upper subsurface hypoxic boundary, is included.

We include four regions for which geospatial analyses were conducted: the Subarctic Pacific (SP; 0–3,200 mbsl; southern latitude limit at 38° 30′ N in the Western Pacific and 49° 30′ N in the Eastern Pacific), the California Current (CC; 0–2,400 mbsl; 31° 40′-49° 30′ N), Mexico Margin (MM; 0–1,200 mbsl; 20°-30° N) and the Humboldt Current and Equatorial Pacific (HC; 0–3,300 mbsl; 10° 30′ N-32° S) ( Fig. 4 ). The Benguela Current (BC) and the Oman and Pakistan Margin (OPM) are discussed, however sediment records from these regions did not meet the criteria established to warrant geospatial analysis. We limit our statements of specific volume and surface area changes to Tables 4 and 5 . These analyses represent the state of paleoceanographic records of deglacial OMZ expansion, including caveats associated with proxy records and limitations of the spatial resolution of core sites; our work may also be used to highlight existing knowledge gaps.

Data visualizations applied to bathymetry, or submarine topography, provide unique perspectives into the distribution of seafloor ecosystems. Here we apply bathymetric masks associated with reconstructed paleoxygenation, wherein mask depths were chosen based on the extent of regional hypoxia across two or more cores at key temporal deglacial intervals. If only a single hypoxia record was available for a deglacial event, a hypoxic mask ±50 m from the core depth was applied. Oxic and hypoxic benthic surface area (km 2 ) and water volume (km 3 ) were quantified for temporal intervals through the deglaciation. To produce geospatial maps of paleoxygenation, we synthesized regional scale patterns in oxygenation, took into account eustatic sea level change, and analyzed the deglacial changes in both hypoxic seafloor surface area and hypoxic water volume from the SRTM30_PLUS global bathymetry dataset [ 50 ] using ArcGIS geospatial software [ 51 ]. Analyses were limited to the continental margin within a 400 nautical mile buffer offshore of the continental coastline. One exception to this is the inclusion of the seafloor within 400 nautical miles around the Galapagos Islands, in order to capture the associated cluster of deep core sites nearby. Within each region, the analysis was limited to the seafloor and water column above a specified isobath, selected below the deepest regional core depth.

Sedimentary archives were then assigned paleodepths based on the modern water depth for the core site, the age model of that core, and deglacial eustatic sea level change. Paleodepths were calculated using estimated global eustatic sea level fluctuations constructed from multiple deglacial sea level records ( Fig. 1 ) [ 5 , 6 ]. We acknowledge eustatic sea level is a simplification of local sea level trends through the deglaciation, which are also impacted by glacio-hydro-isostatic effects [ 48 , 49 ].

Paleoxygenation was assessed based on proxy evidence available for each for core. If proxy evidence indicated hypoxia, we partitioned that signal into intermediate or severe hypoxia groups, based on the biologically-relevant classification scheme identified by Hofmann et al., [ 24 ]. To ensure hypoxia designations are as conservative as possible, we refrained from any hypoxic designation without evidence. In practice, this meant we only designated hypoxia where the proxy evidence was explicit. In the absence of that evidence, we did not designate any hypoxia (i.e., in the absence of laminations, we did not designate a potential range from oxic-intermediately hypoxic). Importantly, we found no instance where multiple proxies within one archive produced conflicting hypoxia reconstructions.

We rely on published chronologies for the deglacial and post-glacial core material, and do not reinterpret any chronologies. Regional archives are assembled on a unifying age axis, which is determined only by previously published chronologies. Sedimentation rates for sites <1000 mbsl span a minimum of 9 cm kyr -1 to a maximum of 500 cm kyr -1 , with the majority sedimentation rate ranging from 20–100 cm kyr -1 . Sedimentation rates for sites >1000 mbsl span a minimum of 4 cm kyr -1 to a maximum of 16 cm kyr -1 . Mass Accumulation Rates (MARs) are published for deep sites (>2200 mbsl) and they range from 2–9 g cm -3 kyr -1 .

Associated data provided here include modern water depth of core extraction (m), latitude and longitude, sedimentation (cm kyr -1 ) or mass accumulation rate (g cm -2 kyr -1 ), the presence of 14 C dating or δ 18 O stratigraphy, and the presence of published hypoxia proxy records (including laminations, microfossils, or geochemistry).

Results

Benguela Current The BC system, also referred to as the Angola-Benguela Current, is the equatorward flowing Eastern Boundary Current of the South Atlantic subtropical Gyre (Fig. 2), associated with high productivity, organic-rich sediments, and seasonal upwelling [152]. Upwelling events, and the export of surface productivity, are linked directly to modern subsurface OMZ structure [153, 154]. The HC OMZ is significantly shallower and less geographically extensive than other systems reviewed here (Fig. 2). The core of the BC OMZ is between 300–400 mbsl, and hypoxic waters extend to as shallow as 50 m [155]. Sediment records from the BC reveal mixed deglacial productivity signals, wherein a few single sites indicate decreased productivity during the LGM [156], while more recent and contradictory work indicates a decrease in surface productivity in the Holocene as compared to the LGM [157–160]. Contributing to this mixed signal, it appears that the productivity center may have moved offshore since the LGM [160]. From the evidence available, it does not appear that the OMZ associated with the BC follows a glacial/interglacial cycle like that which dominates the Pacific, but is both more heterogeneous and directly linked to regional cycles of productivity and upwelling [160].