Paleoceanographic Insights on Recent Oxygen Minimum Zone Expansion: Lessons for Modern Oceanography

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.


Introduction
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]).
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 (preinstrumental) record of natural variability provides a critical analytical and interpretive Figure 1. Deglacial changes in Antarctic temperature (Vostok ice core record, purple line) [146,147], Greenland temperature (GISP2 ice core record, blue line) [106,107], sea level (black line) [5,6] and atmospheric pCO 2 (red line) [3]. 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.
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].
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.

The role and importance of OMZs
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].
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].
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

Oxygenation proxies in paleoceanographic records
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.
Sedimentary structures. Severe hypoxic conditions preclude benthic macrofauna, preventing bioturbation and thereby allowing for the preservation of laminated sediments (Table 2) [25][26][27]. For example, annual sediment laminations (i.e. unmixed, fine-grained sediment displaying distinct, continuous layering) are formed in modern California margin basinal features For the upper panel, regional blocks are defined by black lines to highlight where paleoxygenation reconstructions were completed. Data from World Ocean Atlas [192].  [26,27].

Foraminifera
Above 2 ml L -1 , foraminiferan composition is likely not altered by changes in bottom water oxygenation [199].
Foraminiferan communities in intermediate-severely hypoxic sediments are dense and associated with opportunistic taxa and specific morphologies [200][201][202]. Marker species, affiliated with a narrow oxygenation range, can be used to reconstruct seafloor oxygenation on a very fine-scale [31,203].
Organic flux to the seafloor alters the composition and density of foraminiferan communities [204,205], however interpretation of these community traits is not straightforward in a low oxygen setting (e.g., [206]).
Foraminifera are well adapted to the extreme chemical heterogeneity of oxygenation, methane enrichment, organic flux and sulfur-reducing environments found on continental margins [207].  [19]. The isotopic signal of denitrification reflects regional changes in biologically available N pool [35].
δ 15 N is an indirect record of organic flux, as photosynthesis isotopically fractures the nitrogen pool.
The δ 15 N of particulate nitrogen varies with the degree of surface nutrient utilization, commonly termed productivity [37,38], and water column denitrification [128]. For continental margins, sediment denitrification contributes to the isotopic signal [19,20]. Denitrifying bacteria are facultative anaerobes, meaning they are able to respire either oxygen or nitrogen oxides.
Mo is a biologicallyessential, and is conservative in behavior. The pathway for authigenic Mo accumulation in sediment is unresolved, however it is clear that sulfide concentrations must be >2 μM for active Mo formation [43] Concurrent enrichment of Mo and Re is an indicator of an anoxic depositional environment [41].
Significant Cd enrichment (>2 ppm) is associated with laminated depositional records, a suboxic to anoxic seafloor and the presence of trace presence of H 2 S [42].
Composition of Cd in sinking particles is similar to plankton material; therefore changes in productivity may affect levels of enrichment [41].
Cd is involved with nutrient cycling and scavenging. It is released in sediment during the diagenesis of organic matter, forming insoluble sulfides when trace amounts of H 2 S are present [209]. Cd enrichment occurs in surface sediments [41].
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. A diagram for how multiproxy records are interpreted into paleoxygenation categories is shown in Fig. 3.
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][33][34].
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][37][38][39]. Redox sensitive trace elements such as cadmium (Cd), uranium (U), chromium (Cr) and rhenium (Re) occur in bulk sediment and carbonate No accumulation increase U is enriched in "suboxic" and anoxic sediments as a result of reduction and precipitation (>5μg g -1 ) [210].
Unclear, as the composition of U in plankton is not known.
U is conservative in behavior, and the accumulation of U may be kinetic in nature. U accumulation may be promoted by kinetic effect, in areas of low sedimentation rates [17,211,212]. U enrichment occurs in subsurface sediments [41].
Unclear, as the composition of Re in plankton is not known.
Re is conservative in behavior, diffuses across the sediment-water interface, and the direct mechanism for accumulation is unknown. Free H 2 S is not critical for the accumulation of Re [210]. Concurrent enrichment of Mo and Re is an indicator of an anoxic depositional environment [41].
fossils, and can be used to characterize the redox state of the seafloor as well as functioning as proxies of paleoenvironmental redox changes [40][41][42][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][46][47].

Sedimentary archives
We include cores in this review that met criteria for further geospatial analysis (Table 3). These criteria dictated sediment records that: • exist in a geographic region that includes a modern OMZ, • range in age from approximately 20-0 ka, • have well-constrained chronologies (radiocarbon and/or isotope stratigraphy), • include one or more primary proxies for hypoxia, as described above.
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).
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 .
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.
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].

Bathymetric and Geospatial Analyses
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. We include four regions for which geospatial analyses were conducted: the Subarctic Pacific (SP; 0-3,200 mbsl; southern latitude limit at 38°30 0 N in the Western Pacific and 49°30 0 N in the Eastern Pacific), the California Current (CC; 0-2,400 mbsl; 31°40 0 -49°30 0 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 0 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.

Subarctic Pacific
The SP is dominated at the surface by the eastward flowing North Pacific Current, which branches into the Alaskan Gyre and the southward-flowing CC (e.g., [52]). North Pacific Intermediate Water (NPIW) forms in the Sea of Okhotsk and mixes laterally into the Pacific subtropical gyre [53]. Deep water in the North Pacific is poorly ventilated, nutrient-rich and oxygen-depleted [18], and plays a critical role in maintaining extensive intermediate water OMZs [54,55]. The SP has a modern, seasonal OMZ in the Gulf of Alaska (Fig. 2), attenuated in the summer, with the core of the OMZ from 670-1060 m in the water column [56]. Hypoxic waters spread westward, south of the Aleutian ridge and the Russian Kamchatka Peninsula, during winter months and form a seasonal OMZ that spans the Subarctic Pacific [56]. The SP is a region with inherent complexities, due to the influence of the Cordilleran Ice Sheet [57], the variability of sea ice formation [58], and circulation sensitivity to the episodic closures of the Bering Straight [59]. Previous investigations have characterized the Last Glacial Maximum (LGM) (Fig. 1) as a time of low surface ocean productivity, high ice-rafted debris sediment flux, and cold surface temperatures [60][61][62][63][64][65]. Archives track the deglacial oscillations of the warm Bølling-Allerød (B/A) and cool Younger Dryas (YD), consistent with the expected atmospheric teleconnections between the North Pacific and the North Atlantic [66]. The B/A and Holocene are warmer and more productive across the region, from the Gulf of Alaska [67] to the Western margin of Japan [68]. In step with the changes in surface productivity, seafloor hypoxia developed during the warm, productive intervals of the deglaciation across intermediate (600 mbsl, site RAMA 44PC) and deep (2,900 mbsl, site EW0408-85JC) water depths [58,[67][68][69][70][71][72]. The productivity and oxygenation oscillations of margin and shelf environments to Termination IA, the B/A and the YD reveal the substantial changes this contiguous high latitude environment recently underwent.
Paleoxygenation reconstructions for the Subarctic Pacific. Four sediment cores meet the criteria for deglacial reconstructions: EW0408-85JC, EW0408-11JC, RAMA 44PC, and CH84-14 ( Fig. 4; Table 3). Cores EW0408-85JC and EW0408-11JC were both collected from the Gulf of Alaska, RAMA 44PC was collected east of Kamchatka Peninsula, and CH84-14 was collected east of Hokkaidō Island, Japan. Shifts in seafloor oxygenation throughout the SP exhibit a cohesive, though limited, hydrographic picture (Figs. 5 and 6). We limit our SP reconstruction to mid-way through the deglaciation (14-4 ka), wherein the regional margin is extensively hypoxic, followed by a contraction to shallow, upper intermediate water hypoxia in the Holocene. The Bering Sea and Sea of Okhotsk are not included in the analysis, due to their unique, regional-scale oceanography.  [73]. Although lacking a high enough sedimentation rate necessary to be included in this analysis, ODP Site 882 indicates that the deep SP (>3000 m) gained oxygen through the deglaciation. At Termination IA, the three deepest sites record the development of severe hypoxia, associated with laminations [67] and high authigenic [U] and [Mo] [74] in EW0408-85JC, and high sedimentary δ 15 N values in CH84-14 and RAMA 44PC [68]. These proxy records indicate that at 14 ka severe hypoxia ranged from 596-2,894 mbsl (Fig. 5). Regional hypoxia was followed by subsequent oxygenation in intermediate and deep waters through the YD, Termination 1B and the early Holocene. In the mid-Holocene, the shallowest site recorded the development of severe hypoxia in upper intermediate waters ([Mo]) [74].
Paleoceanographic reconstructions from the SP indicate extensive and severe hypoxia developed during the warming of the deglaciation, and suggest severe hypoxia as shallow as~600 mbsl and across~2,298 m of contiguously hypoxic water column (Fig. 6). The absence of hypoxia during the YD cooling of the Northern Hemisphere indicates that the subsurface SP is sensitive to rapid climatic oscillations. Coherent oxygenation oscillations, from the Gulf of Alaska to Eastern Japan, indicate that geographically widespread changes occur in response to global climate change.

California Current
The CC is the eastern limb of the North Pacific Subtropical Gyre (e.g., [75][76][77]), and is driven onshore in the Northern Hemisphere spring [54]. The CC is characterized by seasonal  [80][81][82]. ODP Site 893A (576.5 mbsl) stands out amongst the region's extensive suite of cores [27,[83][84][85][86][87][88] as a globally recognized site with high-resolution synchroneity to Northern and Southern Hemisphere climate events [40,Figure 5. Subarctic Pacific (SP) deglacial core data synthesized into hypoxia categories. Changing deglacial core depths reflect global eustatic sea level change. The encircled number adjacent to each core label corresponds to the number of available oxygenation proxies, which are enumerate in Table 3. Vertical grey bars correlate to temporal intervals in OMZ geospatial reconstructions for this region.
doi:10.1371/journal.pone.0115246.g005 88,89]. The expansion and contraction of the regional CC OMZ occurred with remarkable fidelity to the glacial terminations, warming intervals, and cooling oscillations of the Northern Hemisphere, with clear regional-scale deoxygenation associated with Termination IA, the B/A, Termination IB and the start of the Holocene (e.g., [27,29,40,42,90,91]). Sediment cores Figure 6. Subarctic Pacific bathymetric seafloor masks and surface area (km 2 ) histograms of deglacial hypoxia impingement for (a) 14 ka, (b) 10 ka, and (c) 4 ka. The Bering Sea and Sea of Okhotsk are not included in the analysis, due to their unique, regional-scale oceanography. Seafloor is selected between 0-3,200 mbsl, with a southern latitude limit at 38°30 0 N in the Western Pacific and 49°30 0 N in the Eastern Pacific with a northern limit along the Aleutian Arc. Analyses were limited to the continental margin within a 400 nautical mile buffer offshore of the continental coastline and the Aleutian Arc. The changing gray shoreline through the panels depicts the paleo-shoreline. At 14 ka, severe hypoxia ranged from 596-2,894 mbsl. At 10 ka, the water column was oxygenated, and at 4 ka severe hypoxia was found between 132-232 mbsl. from Vancouver Island record OMZ expansion from 13.5 to 12.6 ka, slightly delayed from more southern CC sites [92].
Intermediate waters (>1,000 mbsl) in the CC were oxygenated during the LGM, up until Termination IA at 14.7 ka [27,29,40,88,93,94]. The three deepest sites exhibit evidence of hypoxia in the LGM, including ODP Site 1017E and ODP Site 1019E, and F-8-90-G21 (Fig. 7). The ODP cores have elevated [Mo] values in the LGM [42], and hypoxia-associated foraminiferan marker species are preserved in the deepest core [95][96][97]. . California Current (CC) deglacial core data synthesized into hypoxia categories. Changing deglacial core depths reflect global eustatic sea level change. The encircled number adjacent to each core label corresponds to the number of available oxygenation proxies, which are enumerate in Table 3. Vertical grey bars correlate to temporal intervals in OMZ geospatial reconstructions for this region.
Intermediate waters returned to oxygenated conditions midway through the deglaciation (12.9-11.7 ka), synchronous with the ephemeral YD cooling observed in other Northern Hemisphere climate records (Fig. 8) [106,107]. This margin-wide oxygenation is evidenced by the unanimous absence of laminations across all depths, foraminiferal assemblages that reflect increased oxygen concentrations, and low concentrations of redox metals in all regional cores. Termination IB (11.7 ka) initiated a return to regional hypoxia, detected in eight of the cores. Sites within the Santa Barbara Basin preserved laminations at Termination IB [29] and were dominated by hypoxic, OMZ-associated foraminiferal communities [93,102,103]. Termination IB was a secondary expansion of the CC OMZ, and is associated with a narrow band of severely hypoxic water at 436-525 mbsl, bracketed by intermediate hypoxia from 417-525 mbsl and 625-954 mbsl (Fig. 8).
Paleoceanographic reconstructions from the CC reveal the extraordinary shallowness (severe hypoxia <300 mbsl) and extensity (1,233 m of contiguously hypoxic water column) of the regional OMZ during the recent deglaciation. Oxygenated upper ocean ecosystems are dramatically compressed at both Termination events to <300 m from the ocean surface. The YD cooling mid-way through the deglaciation is a remarkably ephemeral event of regional oxygenation. The analysis presented here reveals how geographically pervasive and temporally responsive subsurface oxygen concentrations in the CC system are to global climate change.

Mexico Margin
The MM OMZ is remarkable in thickness and intensity (Fig. 2), and is a product of high surface production, a sharp pycnocline that inhibits local ventilation of subsurface waters, and the "sluggish and convoluted deep circulation" of regional subpycnocline waters [108]. The eastern Pacific sea surface temperature warm pool is centered off of southern Mexico and Guatemala, a product of the large seasonal net heat flux and weak wind mixing of the region [108,109]. The CC and NPIW travel equatorward along Baja California and at 25°N turn westward with gyral circulation [110][111][112], maintaining an oceanographic transition zone between 20°-30°N. In the north, the OMZ sits between 500-1,000 mbsl [56] and in the south shoals upwards (100-700 mbsl), due to the influence of southern sourced intermediate water (including AAIW and "Equatorial Waters") [113]. The oxygen concentration of these waters is low (~0.2 ml L -1 ), and the upper hypoxic boundary is extremely shallow. The term "Equatorial Waters" refers to a composite of Subtropical Underwater found below the equator [114,115], north and south subtropical surface water [108], 13°C waters [116,117] and subtropical mode water [118].
Sediment records in this region indicate major reductions in bottom water oxygenation through the last deglaciation, however the deglacial intensification of the OMZ does not appear to be synchronous across multiple locations [16]. Primary production and carbon export is Oxygen Minimum Zone Expansion thought to have been low during the LGM [119], potentially due to reduced upwelling [17,120]. Changes in the depth of the equatorial thermocline and nutricline have been hypothesized to drive deglacial oscillations in oxygenation [121]. Unlike SP and CC records, oxygenation oscillations are not a comprehensive feature of the margin record, and only appear in records from the Sea of Cortez, which is isolated from the open margin and considered a more continental record of surface and atmospheric conditions [16,26]. The equatorial Pacific is a major global denitrification zone. δ 15 N records from MM sites exhibit glacial-interglacial variability [122] and synchroneity to sites along the Peru-Chile margin [123][124][125].
Paleoxygenation Reconstructions for the Mexico Margin. Seven paleoceanographic records were selected for deglacial reconstructions, including two cores within the Sea of Cortez (GGC-55/JPC-56 and DSDP Site 480), three cores along the southwestern margin of Baja California Sur (MV99-PC14, MV99-PC08, MD02-2508 and GC31/PC08) and two cores along the Mexican Margin (NH8P, NH15P) ( Fig. 4; Table 4). The MM exhibits reduced regional coherency, with asynchronous hypoxia and oxygenation oscillations (Fig. 9). However, this region Figure 9. Mexico Margin (MM) deglacial core data synthesized into hypoxia categories. Changing deglacial core depths reflect global eustatic sea level change. The encircled number adjacent to each core label corresponds to the number of available oxygenation proxies, which are enumerate in Table 3. Vertical grey bars correlate to temporal intervals in OMZ geospatial reconstructions for this region.
doi:10.1371/journal.pone.0115246.g009 exhibits a clear trend of higher oxygenation in the LGM, transitioning to severe hypoxia in intermediate waters at the start of the Holocene.
The two deepest, southernmost sites from the Mexican margin suggest the seafloor was hypoxic during the LGM, as indicated by the presence of laminations (Site NH8P) [126,127] and high [Mo] concentrations (Site NH22P) [119]. At 18 ka, hypoxia was detected at only one site in intermediate waters (898 mbsl). Hypoxic waters spatially expanded along the Mexican margin and in the Sea of Cortez at 14 ka, and ranged from 334-932 mbsl. Site NH8P records the expansion of severe hypoxia along the margin, with laminations preserved between 15-7.5 ka [126,127]. Site GGC-55/JPC-56 preserved laminations from 14.8-12.5 ka, concomitant with high δ 15 N values (>14‰), and high opal accumulations rates [128]. Lamination preservation occurred from 13-11.2 ka in DSDP Site 480, subsequently followed by bioturbation from 11.2-10.5 ka [26].
Paleoceanographic reconstructions for MM provide evidence for the regional intensification of subsurface hypoxia through the recent deglaciation. The MM exhibited reduced sensitivity (Fig. 10) to the first Northern Hemisphere glacial termination event (14.7 ka), as compared to the NP (Fig. 6) and CC sites (Fig. 8). A regional-scale expansion of hypoxia was recorded at 10-11 ka, wherein~600 m of the water column became severely hypoxic (Figs. 9 and 10).

Eastern Equatorial Pacific and Humboldt Current
The HC, also known as the Peru Current, is the eastern limb of the South Pacific subtropical gyre, characterized by upwelling and extreme biological productivity [130]. A thick (~500 m), intense ([O 2 ]<0.2 ml L -1 ) and shallow (upper boundary~50 mbsl) OMZ characterizes the HC (Fig. 2) [131]. Cold HC surface waters move north along the South American continental margin, and then become a part of the equatorial cold tongue. The southward-flowing Peru-Chile Undercurrent is found 75-500 mbsl, in association with the OMZ [132]. This equatoriallyderived water is transported as far as 48°S [133], and is derived from Equatorial 13°C Waters [116,117], Subtropical Underwater [115], and Eastern South Pacific Intermediate Water [134]. The spatial distribution of OMZ thickness is correlated with upwelling conditions [135]. An intense upwelling system is located off of Peru [136,137], which brings nutrient-rich AAIW to the surface and stimulates high production in the equatorial cold tongue [138][139][140]. A more seasonal upwelling cell exists southward off of Chile [141]. A functional break in deep water properties exists at 15-25°S, where tropical and subtropical deep waters are dynamically separated, as deep waters are directly connected in the Western Pacific [142]. Abyssal waters (>4,000 m) are a mixture of water from the Weddell Sea and the North Atlantic, and are relatively high in oxygen concentrations [143]. The deglaciation is thought to have co-occurred with extensive ventilation of the deep-sea, renewing oxygen concentrations in the deep ocean interior (e.g., [2,63,73,144]). Recent syntheses document global increases in deep ocean [O 2 ] through the deglaciation, reflecting the transfer of respired carbon from the deep ocean to the atmospheric and surface ocean carbon pool [7]. Numerous high quality records of deep-sea oxygenation are present in the Eastern Equatorial Pacific, and are included in the HC analysis. Intermediate water records from the HC are limited; however, it is clear that intermediate water (~300-1,400 mbsl) deoxygenated at these depths prior to the Northern Hemisphere glacial Terminations [124,145].  Table 4). Sites exhibit regional synchroneity from deep to upper intermediate water depths, lack coherent timing with Northern Hemisphere climate records (Fig. 11), and exhibit a deglacial temporal signature similar to the Southern Hemisphere (e.g., [146,147]). Broad deglacial patterns include increasing deep ocean oxygenation, severe and shallow hypoxia shoaling at~17 ka in upper intermediate water, and a~6 kya temporal overlap between deep and intermediate hypoxia (Fig. 12).
For the eight equatorial deep waters sites, from 2,203-3,091 mbsl, the LGM is associated with high concentrations of [U] [144,148,149]. These records reconstruct the presence of intermediate hypoxia at 18 ka from 2,082-3,088 mbsl. The deepest, most southerly record (GeoB7139-2) stands out as anomalous to the deep equatorial cores, wherein δ 15 N values Figure 11. Equatorial Pacific and Humboldt Current (HC) deglacial core data synthesized into hypoxia categories. Changing deglacial core depths reflect global eustatic sea level change. The encircled number adjacent to each core label corresponds to the number of available oxygenation proxies, which are enumerate in Table 3. Vertical grey bars correlate to temporal intervals in OMZ geospatial reconstructions for this region.  Seafloor is selected between 0-3,300 mbsl and latitudinally constrained between 10°30 0 N-32°S. Analyses were limited to the continental margin within a 400 nautical mile buffer offshore of the continental coastline and the Galapagos Islands. The changing gray shoreline through the panels depicts the paleo-shoreline. At  Oxygen Minimum Zone Expansion rapidly increased (~5-6‰) at~17 ka, suggesting a shift towards denitrification and an associated decrease in oxygen concentration [123,150]. Deep-sea hypoxia continues to be a regional feature through the deglaciation and into the mid-Holocene, although the timing of post-glacial deep ocean oxygenation is debated. Concentrations of [U] remain high into the mid-Holocene in the equatorial cores [144], beyond 14.7 ka, which differs from the timing of oxygenation in the North Pacific [73]. This temporal difference may be an artifact of increases in productivity and organic matter flux, or it may be a true signal of the temporal differences in oxygenation across the Pacific. For our purposes here, we follow existing interpretations of high [U] in the post-glacial Equatorial Pacific as a signal of intermediate hypoxia (Fig. 11).
Deoxygenation in intermediate waters occurred at~17 ka, as exhibited by laminations, ã 5‰ δ 15 N increase along the Chilean Margin (W7706-41K, W7706-40K, W7706-37K) [145,151], and high δ 15 N values adjacent to the Panamanian Margin (ODP Site 1242; Fig. 11) [124]. Sedimentation is discontinuous through the deglaciation in W7706-37K and may be due to sediment disturbance rather than a signal of oxygenation reversal [145]. Together, these cores indicate severe hypoxia at 13 ka ranged from 108-331 mbsl and intermediate hypoxia ranged The HC provides a unique record of synoptic changes occurring in both deep and intermediate water through the events of the deglaciation. The HC exhibits a coherent regional signal of oxygenation, while lacking synchroneity to the Northern Hemisphere. Striking features of the record include the regional deoxygenation of upper intermediate waters at 17 ka, the extreme vertical expansion of hypoxic water at 13 ka, and the~6 kyr overlap of deep and intermediate hypoxia creating~3,000 m of contiguous hypoxic water column.

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][158][159][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].

Oman and Pakistan Margin
The OPM, within the Indian Ocean, has the most globally extensive OMZ in terms of water column vertical extent (>1000 m) [161] (Fig. 2). Subsurface dissolved oxygen content of the region is controlled by nutrient additions, changes in intermediate and deep water ventilation, and convective mixing of deep waters [162]. Oxygen depleted waters sourced from small inlet seas contribute to the maintenance of severe hypoxia in the region [163]. Monsoonal seasonality of the Arabian Sea drives OMZ oscillations, wherein strong onshore winds force the upwelling of nutrient-rich waters [164,165]. OPM sediments are more organic rich than other lowlatitude OMZ regions, attributed to high primary productivity and reduced remineralization at depth [166].
The OPM is a sensitive recorder of climate variability and an archive of the connectivity between low latitude monsoonal dynamics and high latitude temperature variability [167,168]. Paleoceanographic investigations show that the strength of the OPM monsoon has highresolution, sub-millennial cycles that exhibit synchrony to Greenland climate variability [168]. Regional surface water productivity has co-varied with OMZ intensity, such that OMZ intensity is especially weak during cold stadial events when summer monsoon effects are reduced [167,[169][170][171][172], while OMZ intensity is enhanced during warm interstadials (e.g., the B/A) [171]. Deglacial cores demonstrate rapid (10 2 year timescales) shifts from intermediate to severe hypoxia due to climate forcing [173]. Weakening of the OMZ occurred during Heinrich 1 (~18-14 ka) and the YD [174,175]. Monsoon intensity peaked with a corresponding increase in surface productivity [172,176] and ocean hypoxia from 9.5 to 5.5 ka [177,178]. Geospatial reconstructions of deglacial hypoxia were not conducted for the OPM, as the limited highresolution deglacial records were not sufficient to warrant further analyses.

Mechanisms and implications of OMZ variation
Large-scale physical and biogeochemical processes in the ocean drive OMZ formation, and have high-latitude and low-latitude derived features. As such, the deglacial changes in dissolved oxygen may be explained by proximal and distal mechanisms. High latitude sources for deglacial OMZ variability have been hypothesized to include changes in intermediate water production/ventilation [17,27,40,68,71,84,93,150,179] and changes in deep ocean circulation [70,123,170]. Deglacial micronutrient [Fe] enrichment has been hypothesized to intensify hypoxia, particularly in the Subarctic Pacific [58,67,151]. Additionally, oxygen consumption at the site of NPIW subduction has been hypothesized as a high latitude mechanism [68]. Low latitude mechanisms for deglacial OMZ variability include the relative importance of equatorial counter currents [87,122], the intensity and location of upwelling systems [127,180] changes in the oceanic preformed nutrient inventory [73,124,126,128,181], surface ocean productivity [119,121], atmospheric structural changes [17,85,182], and the strength of monsoonal systems [167,168].
In some cases, independence between mechanisms that force hypoxia has been demonstrated, for example where local productivity and the high-latitude ventilation of intermediate waters are shown to decouple during the deglaciation [90,183]. However, many of the physical and biological processes that drive the development of subsurface hypoxia are ultimately linked. Indeed, circulation models from the latter half of the 20 th century have revealed that changes in the rates of surface productivity can ultimately be driven by physical perturbations in circulation, through interactions between thermocline depth, nutrient flux, and particulate export into hypoxic waters [184].
Major changes in the distribution of [O 2 ] in Eastern Boundary Currents occurred during the recent deglaciation. The substantial oxygen changes arose concomitantly with increasing atmospheric carbon concentration, surface warming and rising sea levels [2][3][4][5][6]. The analyses presented here reveal that rapid oscillations in oxygen distribution are an inherent feature of large ocean regions, as subsurface [O 2 ] exhibits extreme sensitivity to hemisphere-specific warming and cooling. Mid-way through the deglaciation, hypoxic waters shoaled vertically, shifting the upper OMZ boundary towards the ocean surface, with shoaled depths ranging from 596 mbsl in the SP (Fig. 6), 332 mbsl in the CC (Fig. 8), 334 in the MM (Fig. 10), to 108 mbsl in the HC (Fig. 12). Subsurface hypoxia expanded and intensified, creating OMZ environments which ranged in thickness from 2,298 m in the SP, 1,233 m in the CC, 598 m in the MM, to 3,022 m in the HC (Table 4). These data provide evidence of the capacity for OMZs to exhibit extreme shallowness and water column extensity to states that have no analogue in the modern ocean.
The deglaciation offers an extremely informative case study of the sensitivity and coupling between OMZs and global climate oscillations. Our analyses show that during global-scale warming events, vast expanses of the upper global ocean deoxygenate, resulting in the vertical compression of oxygenated ocean ecosystems. These analyses of past changes in dissolved oxygen are fundamentally relevant to the loss of [O 2 ] observed in the modern ocean [185][186][187], and provide bounds and reasonable expectations for expanding OMZs in the future.
Changes in the distribution of oxygen translate directly to the structure of marine ecosystems (e.g., [11,12]) and OMZ expansions have implications for modern oceanography, biodiversity conservation, ocean management and sustainable fisheries. Modern oceanographers can anticipate that Eastern Boundary Current OMZs have the capacity to expand to hydrographic structures that have never been instrumentally observed. These analyses underscore the continued need for high-quality instrumental hypoxic time series data and predictive models of modern oceanographic change. The challenge to resource managers and conservationists is how to plan for the great deal of uncertainty introduced when ecosystems undergo changes that only have precedent in the geologic record.

Modern subsurface oxygen variability
In the modern ocean, the vertical expansion and intensification of OMZ regions has been detected in the California Current [185,188], equatorial waters [186], the Subarctic [189,190] and Subtropical Pacific [191], the North Atlantic [192], the Indian Ocean [193] and the Southern Ocean [194]. The loss of subsurface dissolved oxygen is an acute perturbation to coastal ecosystems and both benthic and pelagic communities [187]. These data indicate that deoxygenation in the 20 th and 21 st Centuries is a feature of every global ocean basin. However, unique physical and biological processes within oceanographic provinces provide a more nuanced and complicated view of climate-forced ocean deoxygenation.
OMZs exhibit high-frequency variability on diurnal and semi-diurnal timescales [195], on the intra-annual timescales of upwelling and relaxation cycles [196], and in tight association in the California Current with La Niña events [197]. These OMZ oscillations illustrate the potential for undescribed scales of variability in the coastal ocean. Additionally, as longer instrumental time series are developed, novel questions regarding intrinsic high-frequency oceanographic variability arise. For example, conflicting interpretations of California Current oxygen trends exist, such that the data can be analyzed to reveal a long-term deoxygenation trend [185] or a 20-25 year undescribed oscillation [188]. This interpretive conflict highlights the importance of understanding high frequency variability of oxygen concentrations, which may be overprinted by anthropogenic climate forcing. Complicating the picture even more so, recent analyses of historical sediments from the Eastern Equatorial Pacific indicate that the regional OMZ has contracted, and that this contraction is correlated with the slackening of trade winds in the tropical Pacific [198]. These regional trends in dissolved oxygen illustrate the potential for OMZ oceanographic provinces to respond differently to anthropogenic climate forcing, however 20 th century oxygenation in the Equatorial Pacific is thought to be anomalous in the context of a broader global picture. On centennial time scales, processes associated with a warm surface ocean, including gas solubility reduction and physical stratification, are predicted to substantially reduce dissolved oxygen concentrations in the ocean interior [8,184].

Conclusions
We integrate existing deglacial geochemical, sedimentary, and microfossil oxygenation proxies to reconstruct the timing, depth and intensity of seafloor hypoxia in Eastern Boundary Currents, principally in the Eastern Pacific. These analyses illustrate the high degree of coupling between the global climate system and OMZ environments and provide the most comprehensive window, to date, into the spatial capacity of OMZ ecosystems to expand and contract due to climate change. The recent deglaciation was accompanied by the dramatic shoaling of the upper hypoxic boundary toward the ocean surface, the compression of upper ocean oxygenated habitat, and the expansion of the subsurface hypoxic water column. Subarctic Pacific and California Current continental margins exhibit tight correlation to the oscillations of Northern Hemisphere deglacial events, whereas the Mexico Margin and the Equatorial Pacific and Humboldt Current exhibit hypoxia expansion prior to Termination IA (14.7 ka), and no regional oxygenation oscillations. Oxygenation changes occurred in synchrony across ocean basins, revealing the extensive sensitivity of upper ocean systems to changes in global climate.