ENSO drives near-surface oxygen and vertical habitat variability in the tropical Pacific

El Niño-Southern oscillation (ENSO) is the leading cause of sea surface temperature variability in the tropical Pacific with known impacts on tuna geographic range, but its effects on oxygen and available oxygenated habitat space are less clear. Variations in oxygenated vertical habitat space in the upper-ocean can alter interactions between predator and prey, as well as drive changes in the vulnerability of economically important tuna and other pelagic fish to surface fishing gear. Using in situ measurements, we show that ENSO is the primary driver of upper-ocean oxygen partial pressure (pO2) variability on year-to-year time scales in the tropical Pacific. Mechanistically, these pO2 variations are primarily caused by vertical shifts in thermocline depth, which alternately elevate and depress cold, hypoxic waters from the ocean interior depending on the ENSO phase and location. Transport-driven, isopycnal pO2 variations within the thermocline also play an important but secondary role. In the western tropical Pacific, waters within the exclusive economic zones of Palau, Micronesia, Nauru, and the Marshall Islands undergo the greatest variations in oxygenated tuna vertical habitat extent: approximately 19.5 m, 23.9 m, 19.5 m, and 19.3 m, respectively, between El Niño and La Niña phases. Oxygen thus plays an important role in altering available tuna vertical habitat space between different phases of ENSO.


Introduction
. Because of the critical role that tuna fisheries play in ensuring economic development and food security in the region (e.g. Barnett 2011), it is important to better understand how tuna respond to ENSO-driven changes in their habitats. Past work has focused on temperature-related variations in observed tuna habitat extents associated with ENSO (Lehodey et al 1997), but habitat favorability can also be affected by changes in oxygen availability (e.g. Barkley et al 1978, Bushnell and Brill 1991, Brill 1994. With predicted declines in global oceanic oxygen content under greenhouse warming (Frölicher et al 2009, Keeling et al 2009, Bopp et al 2013, Cocco et al 2013, Cabré et al 2015, Schmidtko et al 2017, oxygen-driven constraints on available habitat will likely become even more important.
Changes in oxygenated vertical habitat space may have particularly important and far-reaching effects on tuna fisheries and ecosystem dynamics. Hypoxiainduced vertical habitat compression can lead to crowding of tuna closer to the surface, where they are more easily captured and potentially overfished (Green 1967, Barkley et al 1978, Prince and Goodyear 2006. Differential variations in vertical habitat extent among organisms with different oxygen requirements can decouple or enhance predator-prey interactions, altering previously established food webs (Prince et al 2010, Bertrand et al 2011, Koslow et al 2011, Mislan et al 2017. Horizontal migrations away from traditional fishing grounds can also result from changes in vertical habitat suitability (Block et al 2011, Pinsky et al 2013. A better understanding of ENSO-driven variations in vertical habitat extent can therefore help improve seasonal predictions of species spatial distributions, allowing both fishers and fisheries managers to make better informed short-term decisions (Hobday et al 2016). Quantifying the impacts of ENSO on vertical habitat extent can also inform strategies to address longer-term future changes, which include projected increases in the frequency of extreme ENSO phases (Cai et al 2014(Cai et al , 2015. More frequent and/or intense extremes could lead to pronounced year-to-year changes in available tuna habitat that may swamp projections of relatively modest long-term mean trends (Mislan et al 2017). Thus, a deeper understanding of oxygen-defined vertical habitat variability could also increase the accuracy of current and future tuna stock assessments, which are frequently based on catch per unit effort (CPUE) estimates normalized by climatological habitat favorability maps that do not take interannual oxygen variability into account (Brill and Lutcavage 2001, Bigelow et al 2002, Goodyear 2003, Maunder et al 2006, Prince and Goodyear 2006, Maunder and Piner 2015. More recent ecological tuna models do account for the effects of oxygen on tuna movements and populations, using either energy budgets of individual organisms (Dueri et al 2014) or habitat favorability indices computed with theoretical functions bounded by thresholds measured in the lab or field (e.g. Lehodey et al 2008, Lehodey et al 2010, Lehodey et al 2013, Lehodey et al 2015. However, a lack of knowledge regarding oxygen variability and movements within real-world tuna habitats makes it difficult to precisely judge the accuracy of oxygen-influenced spatial dynamics within these models. Coupled ocean-atmosphere models with biogeochemistry can provide some insight into oxygen variability within the tropical Pacific (e.g. Frölicher et al 2009, Cabré et al 2015, Resplandy et al 2015, but without a comprehensive mechanistic understanding of observed variability, these models can only be judged by their skill in reproducing (approximate) climatological signals. These coupled models are also used to project changes in ocean oxygen content with future warming, where the statistical significance of future trends is typically determined by how large the trends are compared to natural variability within the models. If the magnitude of natural, short-term variability such as ENSO in these models is not in line with observations, then their ability to provide robust projections of future conditions may be diminished. Oxygen variability and trends within the tropical Pacific have proven to be especially difficult to model, given the complicated ENSO dynamics, ocean circulation patterns, and particle export fluxes in this region (e.g. Cabré et al 2015). Quantification of observed present-day oxygen variability within the tropical Pacific is therefore needed to better evaluate and improve both ecological and coupled ocean-atmosphere models.
Here we quantify observed ENSO-driven upperocean (0-300 m) oxygen variability and analyze its mechanistic drivers throughout the tropical Pacific. The impact of these ENSO-driven oxygen variations on hypoxia-defined tuna vertical habitat extent is then examined and related to potential effects on fisheries, food security, and economic development mainly in the western tropical Pacific, where tuna catches are greatest.

Methods
2.1. Oxygen partial pressure (pO 2 ) and tuna hypoxic depth (THD) calculations To accurately represent the oxygen conditions experienced by marine animals such as tuna, we report oxygen content in terms of partial pressures (pO 2 ) because it is pO 2 , rather than dissolved oxygen concentration (O 2 ), that drives oxygen transfer and provision to animal tissue (Hofmann et al 2011, Seibel 2011). To compute pO 2 from input variables of temperature, salinity, and O 2 , we first convert O 2 to percent oxygen saturation (Garcia and Gordon 1992). We then divide percent oxygen saturation by the fractional atmospheric concentration of oxygen (21%) and as a final step, correct for hydrostatic pressure at depth (Enns et al 1965). An additional advantage of using pO 2 is that because it is a function of percent oxygen saturation, it automatically accounts for changes in surface oxygen solubility owing to changes in sea surface temperatures (SST); in other words, pO 2 is not affected by variations in solubility, allowing us to rule out this potential confounding factor driving variations in oxygen content. To describe changes in the availability of oxygenated vertical habitat in the upperocean to tuna and other pelagic predators with similar physiologies, we use THD, defined as the shallowest depth at which pO 2 first falls below 15 kPa. This 15 kPa pO 2 threshold is approximately equivalent to the 3.5 ml l −1 dissolved concentration threshold (assuming a salinity of 32.5 psu, a temperature of 23.
Raw profiles of temperature, salinity, and O 2 were binned onto 5-by-5 horizontal degree, monthly mean maps with standard WOD vertical levels (5 m depth resolution down to 100 m, 25 m depth resolution between 100 and 500 m, and 100 m depth resolution beneath 500 m). These monthly mean maps were then used to calculate corresponding pO 2 maps at grid points where all three input variables were available. The mean seasonal cycle was then computed at each grid point for all variables. No lower threshold was placed on the number of years required to compute a climatological monthly mean (e.g. if a given grid cell had only one monthly mean data point in January over the entire time series, then that data point became January's climatological mean at that grid cell). Anomalies were then calculated by subtracting the mean seasonal cycle at each corresponding grid point. The resulting 4D (x, y, z, time) anomaly maps were used for all subsequent map-based analyses (figures 1-4; supplementary figures 3-6, 8,9), where missing data is denoted with empty grid cells.
For all other analyses (table 1; figure 5; supplementary figure 7), raw profiles of the variables were grouped by month and region of interest, including the Western Equatorial Pacific box (WEP), the Eastern Equatorial Pacific box (EEP), and exclusive economic zones (EEZs) (defined in figures 2(a) and (c)). Raw profiles of pO 2 were computed whenever corresponding raw profiles of all three input variables were available. The mean seasonal cycle was then computed within each region of interest, and anomaly profiles were subsequently calculated by subtracting the mean seasonal cycle within each corresponding region.
All aforementioned mapped and regional anomalies were also computed by subtracting out World Ocean Atlas 2013 climatologies (Locarnini et al 2013, Zweng et al 2013, Garcia et al 2013, which did not yield any notable differences in our results. To maintain data source consistency, we thus show only anomalies computed by subtracting out the WODcomputed climatologies discussed above. See supplementary information-data coverage for a discussion of the effects of spatially and temporally variable amounts of data on our results.     figure 2(a)). From left to right, the PNA countries are roughly arranged from east to west. Percentage labels indicate relative deviations of El Niño (red) and La Niña (blue) composite means from the overall mean. Asterisks indicate when El Niño and La Niña composite THD are significantly different from one another within a given country's EEZ, using a two-sided Wilcoxon rank-sum test and a multiple-hypothesis-test false discovery rate of 0.05. The total number of THD measurements available from individual profiles in each EEZ is listed in black below the country names, while the number of measurements available during El Niño and La Niña are listed in red and blue, respectively. Table 1. Thermocline depth (TCD) and tuna hypoxic depth (THD) variations and relationships. Mean and ENSO composite TCD and THD values within the Western Equatorial Pacific (WEP) and Eastern Equatorial Pacific (EEP) regions, as defined in figure 2(a). All-profile THD versus TCD regression coefficients are computed by relating THD and TCD anomalies from all simultaneous profiles (using iteratively reweighted least-squares with a bisquare weighting function and corresponding pseudo R 2 values-see Willett and Singer 1988), and thus represent the magnitude of THD changes driven primarily by TCD changes alone. ENSO-driven THD versus TCD regression coefficients are computed as the difference between El Niño and La Niña THD values divided by the difference between El Niño and La Niña TCD values, and thus represent the magnitude of ENSO-related THD changes driven by both ENSO-associated TCD and isopycnal pO 2 changes. One potential driver of near-surface oxygen variability is changes in the vertical structure of differentiated water masses. Subsurface temperature data is much more plentiful than that of O 2 (see supplementary information-data coverage), and thermo-and oxycline depths covary strongly throughout the tropical Pacific; we therefore use variations in TCD alone to represent the vertical movements of water masses with distinct temperature and oxygen signatures. TCD was calculated from depth-resolved monthly mean temperature maps (figures 1, 3-4; supplementary figures 4, 6(d), 9(b), (d), (f), (h)) and profiles (table 1) using the variable representative isotherm method, as recommended for tropical waters by Fiedler (2010), who compared several objective methods of computing TCD with empirical TCD estimates made by eye.
Other objective methods of computing TCD, such as the depth of the 20°C isotherm or the depth of the maximum vertical gradient in temperature, led to similar results here. In addition to its computational efficiency, the advantage of the variable representative isotherm method is that it accounts for the wide range of temperatures in which the thermocline occurs, since there is no single isotherm that well represents empirical TCDs (e.g. empirical thermocline temperatures range from 16°C to 26°C within the tropical Pacific alone) (Fiedler 2010). According to the variable representative isotherm method, the thermocline spans the water column from the base of the mixed layer to the depth at which temperature has dropped halfway toward the temperature at 400 m; TCD is then defined as the midpoint of this layer. Mathematically, the isotherm at which the TCD occurs is as follows: where T(MLD)=SST − 0.8°C.

ENSO definition and index
All variables (either binned onto 5°-by-5°monthly mean maps or spatially averaged over specified regions) were examined for ENSO-related variability using the Oceanic Niño Index (ONI). ONI is calculated as the 3 month running mean of SST anomalies in the Niño 3.4 region (5°N-5°S, 120°-170°W) (Huang et al 2016). El Niño and La Niña periods occur when ONI is greater than 0.5°C and less than −0.5°C, respectively, for at least five consecutive 3 month running mean periods (supplementary figure 2). Time series of the ONI were downloaded from http://origin.cpc.ncep.noaa.gov/ products/analysis_monitoring/ensostuff/ONI_v5. php on 4 June, 2018.

Statistical significance of spatial data
To avoid incorrect or overstated interpretations, we control our false discovery rate (FDR) whenever analyzing the results of multiple-hypothesis tests. The FDR is defined as the statistically expected fraction of local (i.e. at a single grid point in the case of map analyses or within a single region in the case of comparisons across EEZs or other specified areas) null hypothesis test rejections for which the respective null hypotheses are actually true. Controlling FDR in the context of multiple-hypothesis tests thus requires smaller p-values compared to single hypothesis testing in order to reject local null hypotheses. The procedure for determining significance under a controlled FDR is as follows (Wilks 2016): (1) sort the collection of p-values from N (i.e. the total number of grid points in the case of map analyses and the total number of regions in the case of comparisons across EEZs or other specified areas) local hypothesis tests pi , with i=1, K, N, in ascending order.
(2) Reject local null hypotheses if their respective p-values are no larger than a threshold level p FDR , which is equal to the largest pi that is no larger than the fraction of α FDR specified by i/N. p FDR is thus calculated as follows: where α FDR is the chosen FDR (expressed as a fraction).
This method assumes that the multiple local tests are statistically independent, but is also valid even when the results of the multiple tests are strongly correlated (such as is the case with geospatial data). Indeed, when spatial correlations are high, the achieved FDR will be even smaller (stricter) than the chosen FDR (Ventura et al 2004, Wilks 2006, so that the chosen α FDR should be approximately double the desired level when analyzing highly spatially correlated data grids (Wilks 2016).

Code availability
The MATLAB code required to make all of the figures and tables generated here can be found at https://doi. org/10.5281/zenodo.2648131.

Data availability
Data in the form of * .mat files required to the make all of the figures and tables generated here can be found at http://doi.org/10.5281/zenodo.2648124.

Results and discussion
3.1. ENSO-driven upper-ocean pO 2 and oxygenated vertical habitat variability In the equatorial Pacific (between 7.5°N and 7.5°S), mean near-surface pO 2 differs substantially from east to west, with higher pO 2 values at any given depth in the west compared to the east ( figure 1(a); supplementary figures 3(a), (b)). The mean oxycline, where pO 2 changes most quickly with depth, is located at around 100 m depth near the western boundary, descends to 150 m near the dateline, and shoals to around 50 m near the eastern boundary ( figure 1(a)). Mean THD (figure 1 solid gray contour) lies close to the oxycline across the entire basin. The within-oxycline vertical pO 2 gradient is strongest in the east, coincident with the top of the shallow oxygen minimum zone there.
Monthly variability in equatorial Pacific pO 2 shows distinct east-west differences as well ( figure 1(b); supplementary figures 3(c), (d)). Near-surface pO 2 variability is greater in the east due to the presence of sharp vertical gradients in pO 2 , with maximum variability occurring within the relatively shallow oxycline between 140°W and 80°W ( figure 1(b); supplementary figures 3(c), (d)). pO 2 variability in the west, between the western boundary and 160°W, is about half as large as that in the east and is again greatest surrounding the local oxycline ( figure 1(b); supplementary figures 3(c), (d)).
Interannual variability in near-surface pO 2 is driven primarily by ENSO throughout the tropical Paci-  (Yang et al 2017). In the west, the effects of ENSO on near-surface oxygen conditions have not been wellcharacterized before this study. In both the east and west, ENSO-driven pO 2 variability is greatest within the vicinity of the oxycline (around 100-150 m depth in the west and 50 m depth in the east) (figures 1(c), (d)), as is the case for total interannual pO 2 variability ( figure 1(b)). The ENSO-associated pO 2 variability dominates the total variability, in agreement with model results from Resplandy et al (2015), who showed that ENSO is the major driver of interannual variations in tropical Pacific air-sea oxygen fluxes within three Coupled Model Intercomparison Project 5 (CMIP5) Earth System Models (CESM1-BGC, GFDL-ESM2G, and GFDL-ESM2M).
Interannual variability in tropical Pacific THD is also driven primarily by ENSO. Relative ENSO-driven THD variations are greatest from the western boundary (around 125°E) to 160°W and from 130°W to the eastern boundary (around 80°W) (figures 2(b), (c)). Within the western region (125°E to 160°W), THD is as much as 22 m shallower during El Niño versus mean conditions, and 34 m shallower during El Niño versus La Niña conditions (figures 1(c), (d) contours). Averaged over the WEP defined in figure 2(a), THD variations compress vertical habitat space by 10.3% (14.5 m) during El Niño (figure 2(b); table 1) and expand it by 9.5% (13.3 m) during La Niña (figure 2(c); table 1), compared to mean conditions. Vertical habitat space is thus compressed by 27.8 m during El Niño relative to La Niña within the WEP (p<0.05, unequal variances t-test). In the east, ENSO acts in the opposite direction, shoaling THD during La Niña and deepening it during El Niño. Averaged over the EEP defined in figure 2(a) and compared to mean conditions, hypoxia-induced vertical habitat compression of 2.9% (1.1 m) occurs during La Niña (figure 2(c); table 1), while expansion of 12.9% (4.9 m) occurs during El Niño (figure 2(b); table 1). Vertical habitat space is thus compressed by about 6 m during La Niña relative to El Niño within the EEP (p<0.05, unequal variances t-test).

ENSO-associated mechanistic drivers of oxygen variability
Based on the observations analyzed here, ENSOrelated variations in oxygenated tuna vertical habitat space appear to be driven primarily by the up-anddown motions of distinct water masses separated by well-defined thermo-and oxyclines. ENSO has different effects on the vertical positioning of these water masses and therefore TCDs depending on the region. In the WEP (as defined in figure 2(a)), El Niño shoals TCD by 12.2% (18.1 m) (supplementary figure 4(b); table 1), while La Niña deepens TCD by 9.0% (13.3 m) compared to mean conditions (supplementary figure 4(c); table 1). TCD is thus approximately 31.5 m shallower during El Niño relative to La Niña within the WEP (p<0.05, unequal variances t-test). Within the EEP (as defined in figure 2(a)), El Niño deepens TCD by 20.4% (10.5 m) (supplementary figure 4(b); table 1), while La Niña shoals TCD by 15.6% (8.1 m) compared to mean conditions (supplementary figure 4(c); table 1). TCD is thus about 18.5 m shallower during La Niña relative to El Niño within the EEP (p<0.05, unequal variances t-test).
The strongest temporal correlations between pO 2 and TCD (figure 3(a)) occur at depths surrounding the mean thermocline, where both pO 2 and TCD variability are largest (figures 1(b)-(d); 3(a)). The prevalence of these strong positive correlations between pO 2 and TCD, in addition to the strong positive correlations between THD and TCD throughout the tropical Pacific (figure 3(b); table 1), suggest the following: as the thermocline deepens, hypoxic waters beneath the thermocline are pushed downward, expanding oxygenated vertical habitat space and increasing pO 2 in the upper water column; conversely, as the thermocline shoals, hypoxic waters beneath the thermocline are pulled upward, compressing oxygenated vertical habitat space and decreasing pO 2 in the upper water column.
Though first-order variations in THD anomalies are well-explained by the up-and-down motions of the thermocline and accompanying water masses, concurrent changes in pO 2 along isopycnals within the thermocline-caused by changes in primary production and physical transport-can either reinforce or dampen this primary effect (isopycnals are roughly equivalent to isothermals in the tropical Pacific because density below the mixed layer is primarily determined by temperature there). During El Niño, decreases in equatorial primary and export production in the east should decrease microbial oxygen consumption and thus increase pO 2 along isopycnals; in the west, increases in local production and export should lead to a decrease in pO 2 along isopycnals. During La Niña, these conditions should reverse (Chavez et al 1999, Behrenfeld, et al 2001, Lehodey 2001. Contrary to these expectations, however, isopycnal pO 2 within the thermocline decreases in the east and increases in the west during El Niño ( figure 4(b); supplementary figure 5(b)), with the opposite occurring during La Niña (figure 4(c); supplementary figure 5(c)). The isopycnal pO 2 variations observed here are therefore likely caused by ENSO-associated changes in lateral ventilation and transport within the thermocline rather than local changes in primary production. ENSO-associated lateral transport changes are, in turn, likely driven by variations in the strength of subsurface zonal currents, the most important being the Equatorial Undercurrent (EUC). The EUC carries oxygen-rich water from the subtropics and western Pacific to the central and eastern Pacific along the thermocline (Stramma et al 2010), and undergoes large variations in strength associated with ENSO (Johnson et al 2000, Izumo 2005. During the 1997/98 El Niño, for instance, EUC mass transport (measured along the equator at 140°W) dropped to nearly zero from a mean of around 34 Sv (Izumo 2005). Decreases in EUC strength of this magnitude greatly reduce within-thermocline transport of oxygen from the western to the central and eastern equatorial Pacific during El Niño. Increases in EUC strength during La Niña, on the other hand, bring anomalously large amounts of oxygen from the west to the east, raising and lowering within-thermocline isopycnal pO 2 values in the east and west, respectively. During both ENSO phases, these transport-driven isopycnal oxygen changes dampen and oppose the effects of TCD variations, but because they have much a weaker impact on pO 2 and THD, the vertical movement of water masses remains the dominant factor influencing ENSO-associated pO 2 and THD anomalies.
The effects of transport-driven changes on ENSOrelated THD and pO 2 variations are particularly small in the western Pacific (figures 4(b), (c)). As a result, ENSO-induced TCD variations alone can accurately predict concomitant THD variations, leading to comparable all-profile and ENSO-driven THD versus TCD regression coefficients within the WEP (table 1). All-profile THD versus TCD regression coefficients are computed by relating THD and TCD anomalies from all simultaneous profiles, and thus represent the magnitude of THD changes driven primarily by changes in TCD. ENSO-driven THD versus TCD regression coefficients are computed as the difference between El Niño and La Niña THD values divided by the difference between El Niño and La Niña TCD values, and thus represent the magnitude of ENSOrelated THD changes driven by both ENSO-associated TCD and isopycnal pO 2 changes. In contrast to the situation in the west, ENSO-driven variations in eastern Pacific, within-thermocline isopycnal pO 2 noticeably dampen the THD variations induced by vertical shifts in water masses and TCD (figures 4(b), (c)). Because of this, variations in TCD alone (that is, not accounting for within-thermocline isopycnal pO 2 changes) slightly overpredict differences in THD between ENSO phases within the EEP, as evidenced by a larger all-profile THD versus TCD regression coefficient compared to the ENSO-driven one here (table 1). Despite the noticeable dampening effect of physical transport changes on pO 2 and THD variations in the EEP, the up-and-down movements of water masses remain the primary factor driving ENSO-related pO 2 and THD variations, as can be seen from the net direction and sign of these variations. In sum, based on observations alone, the vertical motions of the thermocline appear to be the main driver of ENSOinduced variations in near-surface pO 2 and THD throughout the tropical Pacific, while ventilation-driven variations in isopycnal pO 2 within the thermocline play an important (especially in the east) but secondary role. This is again in agreement with model results from Resplandy et al (2015), who found that the main drivers of interannual variations in tropical Pacific air-sea oxygen flux were changes in physical circulation and vertical transport rather than primary production.

Asymmetries between El Niño and La Niña
In both the WEP and EEP regions (as defined in figure 2(a)), El Niño induces greater deviations of THD from mean conditions compared to La Niña (table 1), which is consistent with observed asymmetrical effects of ENSO on other well-studied climate phenomena (Hoerling et al 1997, Kang and Kug 2002, An et al 2005. Averaged over the WEP, El Niño shoals THD by 14.5 m, while La Niña only deepens THD by 13.3 m (table 1). In the EEP, El Niño deepens THD by 4.9 m, while La Niña only shoals THD and thus compresses tuna vertical habitat space by 1.1 m (table 1). These asymmetric variations in THD are caused by corresponding asymmetric variations in TCD between opposing ENSO phases, such that ENSO-driven TCD changes are also larger during El Niño compared to La Niña within both the WEP and EEP regions (table 1).

Oxygenated vertical habitat variability within western tropical Pacific EEZs
We now narrow our focus to the western tropical Pacific both because it is the most tuna-rich region in the world and ENSO-driven oxygenated vertical habitat variations are most pronounced here. THD variations of the magnitude observed here can significantly affect catchability and CPUE using common industrial surface fishing gear such as longlines and purse seines, potentially resulting in overfishing in years when vertical habitat space is compressed (Green 1967, Barkley et al 1978, Evans et al 1981, Nakano et al 1997, Prince and Goodyear 2006, Bigelow and Maunder 2007, Prince et al 2010. Though tuna occur throughout the Western and Central Pacific Ocean, approximately 60% of the total tuna catch in this region is taken from within the EEZs of only eight Pacific island nations that jointly enacted the Nauru Agreement to better manage their fisheries (Havice 2013). Of the eight Parties to the Nauru Agreement (PNA), the EEZs belonging to Micronesia and the Marshall Islands undergo the greatest reductions in mean vertical habitat space during El Niño (14.0 m or 9.9% and 13.4 m or 7.9%, respectively), while the EEZs belonging to Palau, Papua New Guinea, and Nauru exhibit the greatest expansions in mean vertical habitat space during La Niña (14.3 m or 12.5%, 11.2 m or 8.1%, and 14.5 m or 11.1%, respectively) (figure 5). The EEZs of Palau, Micronesia, Nauru, and the Marshall Islands undergo the largest total mean ENSO-driven variations in vertical habitat space, with 19.5 m, 23.9 m, 19.5 and 19.3 m differences in THD between El Niño and La Niña, respectively (figure 5).

Conclusions
Past work has demonstrated that temperature is an important predictor of tropical Pacific tuna habitat suitability and migration towards the east during El Niño years (Lehodey et al 1997), but the role of oxygen was previously underappreciated due to a lack of knowledge regarding its large-scale ENSO-driven variability. Here we have shown that oxygen-induced vertical habitat expansion in the east and compression in the west acts in the same direction as temperature to potentially push skipjack and yellowfin tuna eastward in times of El Niño (and especially so during an extreme El Niño), as they seek out available food sources and more favorable oxygen conditions. Temperature as well as pO 2 therefore likely play important roles in determining tuna spatial distributions with ENSO. The large ENSO-driven variations in oxygen content within the tropical Pacific suggests that it may be important to account for interannual oxygen variability in predictions of skipjack and yellowfin tuna spatial distributions and stock assessments calculated from CPUE. Depending on the level of accuracy needed, the use of climatological mean oxygen values or long-term mean oxygen trends, may not be accurate or representative enough to capture ongoing tuna spatial dynamics, thresholds, or preferences of interest.
Here we have specifically defined vertical habitat availability in the tropical Pacific using THD computed with a minimum oxygen threshold applicable to skipjack and yellowfin tuna (15 kPa). THD can, however, also be used to define vertical habitats in other areas of the ocean and for other tuna species when computed with different species-specific thresholds. Bigeye tuna, for example, can spend considerably longer in low-oxygen waters than skipjack and yellowfin tuna (e.g. Lowe et al 2000, Evans et al 2008, Schaefer and Fuller 2010, while albacore tuna have oxygen requirements more similar to those of skipjack and yellowfin tuna (e.g. Laurs and Lynn 1977, Childers et al 2011, Williams et al 2015. THD variability for these other species may look very different from those for skipjack and yellowfin tuna, with variations in ocean oxygen content impacting their vertical habitats, catch rates, and stock assessments in different ways within different ocean regions. Within the tropical Pacific, however, it is likely that ENSO impacts the oxygenated vertical habitat suitability of many tuna species, given its large effects on oxygen throughout the region.
With climate change, the frequency of extreme ENSO events is predicted to increase (Cai et al 2014(Cai et al , 2015, potentially leading to greater variability in oxygenated vertical habitat space and therefore less stable year-to-year supplies of tuna and other large pelagic fish. Future warming may also shift the tropical Pacific towards a more El Niño-like mean state (Vecchi et al 2006, Zhang and Song 2006, Vecchi and Soden 2007, Collins et al 2010, Huang and Ying 2015, Ying et al 2016. If these predictions are accurate, vertical habitat space will shrink in the west, pushing skipjack and yellowfin tuna eastward and upward and rendering them more vulnerable to overfishing with industrial surface fishing gear.
In the western tropical Pacific, whether during El Niño years or more permanently under future climate change scenarios, habitat suitability changes and eastward migration of skipjack and yellowfin tuna out of western EEZs threaten food security from local tuna catch, economic development from lucrative tuna canning and port operations, and government income from licensing fees paid by foreign fishing interests to access EEZ waters (Hamnett and Anderson 2000, Moss 2007, Guillotreau et al 2012. Continued and increased in situ monitoring of oxygen in the tropical Pacific is therefore crucial for detecting current and future changes in oxygen conditions that will affect the spatial distributions and populations of tuna and other pelagic fish. Because of the oftentimes dramatic societal impacts brought about by ENSO, changes in ENSO frequency and intensity, on top of changes in mean conditions, are of great immediate concern for tropical Pacific societies. Future work should therefore analyze model-projected, climate warming-driven changes in ENSO-related oxygen variability and extremes.
Washington Research Foundation Fund for Innovation in Data-Intensive Discovery and the Moore/ Sloan Data Science Environments Project. PMEL contribution no. 4953.