Climate and Human-Driven Variability of Summer Hypoxia on a Large River-Dominated Shelf as Revealed by a Hypoxia Index

Coastal hypoxia has become common especially in large river dominated coastal ecosystems. To better quantify the severity of hypoxia and the contribution of hypoxia drivers, we applied principal component analysis (PCA) on observable properties from eight summer hypoxia events in the East China Sea and defined the first principal component as the hypoxia index (HI). Multiple linear regression showed that the HI significantly correlated with three direct hypoxia drivers including water column stratification, subsurface water residence time, and respiration rates, which accounted for 5.7, 55.3, and 34.5%, respectively, of the total variance of PCA derived HI. We further reconstructed the HI over the past 60 years using available long-term data of stratification, model-derived residence times and respiration rates. The results show that summer hypoxia has become more severe since the 1960s. ENSO and global warming may have exacerbated hypoxia by affecting the river discharge, resulting in freshening in the plume-impacted shelf area, while anthropogenic activities may have exacerbated hypoxia by elevating fluvial nutrient concentrations, resulting in higher respiration rates. In addition, warming of the bottom water from the Kuroshio Current accounts for an additional increasing rate for HI, which made hypoxia more severe by means of decreasing oxygen solubility. Overall, our results indicate that stratification, water residence and oxygen solubility resulting from climate change can explain about 80% while higher respiration resulting from higher nutrient inputs can explain about 20% of the variation in the severity of hypoxia during the past half century.


INTRODUCTION
Coastal hypoxia has become common especially in large river dominated coastal ecosystems, where a substantial amount of oxygen is consumed in subsurface waters via oxidation of organic matter coming from either external sources or internal primary production (Breitburg et al., 2018). Recent studies show hypoxia has become more wide-spread or severe in the last two decades; the number of reported global "dead zones" is increasing (Breitburg et al., 2018), the area of hypoxia is expanding, and the minimum values of dissolved oxygen (DO) concentration are decreasing, even approaching anoxia. However, using a single property, such as hypoxic area, may be biased and insufficient to elucidate the severity of hypoxia. Other properties, such as thickness and minimum and average oxygen concentrations, are also crucial for exploring the severity of the hypoxia problem (Matli et al., 2018;Scavia et al., 2018).
Climate change and human activities influence coastal hypoxia via immediate or direct drivers and peripheral or indirect drivers. By immediate or direct drivers, we mean those factors that consume O 2 or limit the supply of O 2 to the bottom water (Fennel and Testa, 2018). The former includes water column and benthic organic carbon respiration (Resplandy, 2018), while the latter includes vertical advection and diffusion Markus Meier et al., 2018) (as characterized by water column stratification), lateral advection (Fennel and Testa, 2018) [as characterized by subsurface water residence time (τ)], and oxygen solubility (Lui et al., 2014). By peripheral or indirect factors, we mean those that lead to the direct factors. For example, eutrophication or high primary production in surface waters can lead to high respiration and thus hypoxia in bottom waters, and temperature and/or salinity gradients can lead to water column stratification. Most of the research on the coastal hypoxia mechanisms has focused on indirect factors, such as eutrophication driven by human activity (Rabalais et al., 2014;Wang et al., 2016) and ensuing enhancement of primary production and organic matter input Li et al., 2018), or temperature and non-temperature effects driven by climate change (Altieri and Gedan, 2015;Li et al., 2016a). We argue that the relative contribution of each of these factors is still not clear. Because these indirect factors influence hypoxia via direct pathways, we must disentangle the contributions of the direct drivers in order to quantify the contributions of climate change and human activities and to delineate the pathways impacting the long-term variation of coastal hypoxia.
The relative contributions of direct drivers controlling longterm hypoxia also vary due to location and climatological regimes of local marine ecosystems. In this study, we have developed a hypoxia index (HI) as defined by the first component of principal component analysis (PCA), using data of observable properties from eight cruises in the Changjiang plume-impacted shelf area in the East China Sea (ECS). We then use a regression analysis of the derived HI and three direct hypoxia drivers (stratification, τ, and biological respiration) to quantify their relative contributions to hypoxia severity. As these direct factors were not available for long-term trend analysis, we used long-term observed data and model results, to reconstruct the HI and to explore the contribution of climate change (associated with stratification, τ, and oxygen solubility) and human activities (associated with biological respiration) to variation in hypoxia over the past 60 years. This study also has implications for understanding hypoxia in estuaries and bays affected by substantial nutrient runoff and climate change.

Study Area and Sampling Collection
Changjiang River (also Yangtze River) is the largest river in the Euro-Asian continent, freshening the area continuously (Park et al., 2011) with its annual discharge of 924.8 × 10 9 m 3 yr −1 . Changjiang plumeimpacted shelf area is one of the most productive continental shelf sea in the world due to the mixing of nutrient-laden Changjiang River plume and Kuroshio Current. Summer hypoxia occurs with strong stratification and prolonged residence time for subsurface water (Zhang et al., 2019), and also fueled by marine organic matter exporting from the surface (Wang et al., 2016). Hypoxia has been recorded since 1950s, and has become more severe recently with DO minimum decreasing and hypoxic area expanding (Zhu et al., 2011;Wei et al., 2017).
We  Figure 1). For this study, we focused on an area between 30.5 • N and 32 • N, 122.5 • E and 125 • E, and assembled oceanographic data over the past 60 years. This is the only transect in the region where long-term data are available.
At each station, temperature and salinity were measured with a Sea Bird model 911 conductivity-temperature depth (CTD) recorder. Water samples for measuring concentrations of nutrients, chlorophyll a (Chl a), and DO were collected with Niskin bottles attached to the CTD rosette.

Observable Properties for Extended Definition of Hypoxia Severity
Since the status of hypoxia is not merely reflected by the area or oxygen concentration, it is necessary to find an index to express the hypoxia severity more comprehensively. In this study, we derive the index based on a PCA of six observable properties (Supplementary Table 1): sectional area, mean thickness, the longitude span of the hypoxia area (Figure 1), DO minimum, mean bottom DO, and mean DO concentration when saturation (DO%) was below 100%. Since the conventional hypoxia threshold of 2 mg L −1 is below the empirical sublethal and lethal DO thresholds for many organisms (Vaquer-Sunyer and Duarte, 2008), we defined hypoxia area as the vertical sectional area surrounded by the 94 µmol L −1 (3 mg L −1 ) contour along the 31 • N transect, the same threshold for hypoxia as previous reports defined Levin et al., 2009;Zhou et al., 2017). The mean thickness of the hypoxia area (D hy , m) is the averaged vertical distance between the 94 µmol L −1 contour and the bottom of the stations. The longitude span of the hypoxia FIGURE 1 | A conceptual model for the dimension of sectional hypoxia area and water and salt budget. V Q , volume of Changjiang River runoff; V g , volume of groundwater runoff; V z , vertical exchange volume between layer 1 and 2; V ent , volume of advection-induced entrainment; V s , surface advective outflow volume from layer 1; V d , deep advective inflow volume to the layer 2; S 1 , salinity of layer 1; S 2 , salinity of layer 2; S 1oc and S 2oc , are oceanic endmember salinities of layer 1 and 2, respectively. area is the maximum distance of the 94 µmol L −1 contour along 31 • N transect.

Calculating the Direct Drivers Using the LOICZ Model
In order to quantify the contributions of the direct drivers in regulating hypoxia, we calculated τ and respiration rates during each cruise, using a two-layer river dominated LOICZ (Land-Ocean Interactions in the Coastal Zone) model (Figure 1), which was based on water and salt balances within each layer at steady state (Gordon et al., 1996). Layers 1 and 2 were divided by the mixed layer depth (MLD). We assumed the Changjiang River and the Kuroshio Current are the two major water masses mixing in the studied system. As groundwater data is very limited and the distribution is patchy, we assumed fresh groundwater only contribute to the layer 1 and a constant flux (V g ) of 0.02 m 3 m −2 d −1 (Guo et al., 2020) with average salinity (S g ) of 5 and phosphate of 3.3 µmol L −1 . Thus, water and salt balances in the entire system are: The salt contributed by Changjiang River water and groundwater were far less than deep advective inflow (V d ) [salinity of Changjiang River water (S Q ) at the Xuliujing station (Supplementary Figure 1) was usually <0.15 (unpublished data)] and the total discharge of Changjiang River water and groundwater was far less than V d , then the last two items (V Q × S 2Q + V g × S 2g ) of Eq. 2 could be ignored.
A combination of Eqs1 and 2 gives: For the bottom layer, the salt balance is: The variables are defined in Figure 1.
Stratification ( σ) was defined as the average difference in the density anomaly between bottom and surface (0-2 m depth) at stations along the transect. We used the simple linear equation of seawater state (Intergovernmental Oceanographic Commission [IOC], 2010) to estimate the density difference between bottom and surface seawater, where T s and T b are surface and bottom temperature ( • C), S s and S b are surface and bottom absolute salinity, and p s and p b are surface and bottom pressure (dbar), respectively. For reducing the observational errors of S s and S b , we used historical data from the Japan Meteorological Agency, Oceanographic Data Center (JODC), the Korea Oceanographic Data Center (KODC), and the World Ocean Atlas (WOA), which were averaged in the specified box (122.5 • E-125 • E, 30.5 • N-32 • N, Supplementary  Figure 1), besides of our own data. MLD is defined as the depth at which the density differs by 0.125 from the surface value (Huang and Russell, 1994). Average water residence time is typically expressed as the ratio of system volume to the average rate of outflow through the system (Gordon et al., 1996). For layer 2, there are two types of water mixing; one is the vertical exchange with layer 1 including diffusion and advection (including entrainment) (V z ), and the other is layer 2 water exchange with the offshore subsurface water (V d ). The latter often accounts for most of the water mixing. If so, then τ is the total volume of subsurface water divided by the water inputs (outputs) or V 2 /V d , can be simply presented by combining Eq. (3) for canceling out V d : The long-term τ can be calculated based on Eq. 6, V 2 was obtained by multiplying the depth below MLD and assuming a fixed study area of 1.8 × 10 10 m 3 along the transect (122 • E-124 • E, 30.5 • N-31.5 • N). The last direct driver is the respiration rates in the water (R h , mmol C m −2 d −1 ), which we derive from the difference between NPP (net primary production) and NEP (net ecosystem production), as NPP = GPP−R a and NEP = GPP−(R a + R h ). GPP is gross primary production, and R a is autotrophic respiration rates. NPP data were acquired from the Marine Satellite Data Online Analysis Platform (SatCO 2 ) 1 , including monthly SeaWiFS VGPM NPP from 1997 to 2007 and MODIS Aqua VGPM NPP from 2002 to 2016. These NPP products were corrected by euphotic layer depth including the effect of colored dissolved organic matter absorption (He et al., 2017). NEP was calculated using the LOICZ model, as the net biological uptake or release of carbon calculated from phosphate stoichiometrically (C/P = 106) (Gordon et al., 1996) within each layer.
Thus, we have, R h of 37.0 mmol C m −2 d −1 from in situ incubation was applied for 2014 due to too much cloud for deriving satellite NPP in the study area. The mean values and uncertainties (standard deviation) of satellite NPP and LOICZ model derived NEP are 68.8 ± 7.9 and 28.2 ± 68.9 mmol C m −2 d −1 , respectively, based on the eight cruises with an error propagation method (Taylor, 1997). The mean value and uncertainty of R h is estimated to be 44.7 ± 70.2 mmol C m −2 d −1 . The volumetric respiration rates (mmol C m −3 d −1 ) was assumed constant over the whole water column in the hypoxia area. This approach is verified and deemed acceptable based on our results during 2011 cruise at five stations of 24 h-incubation for respiration along the transect in Changjiang plume-impacted shelf area (averaged 30.1 ± 24.8, 35.5 ± 19.2, and 26.4 ± 21.5 for 2 m, 10 m, and bottom waters, respectively). The integrated respiration rates in the hypoxia area (R hy , mmol C m −2 d −1 ) depends on the thickness of the hypoxia area and the water depth (D, m) over longitudes (L) spanning the sampling stations from west (Lw) to east (Le): The contribution of sediment to R hy was assumed to be 50% (Zhou et al., 2017) and an average 8.43 mmol O 2 m −2 d −1 of sediment oxygen consumption (Song et al., 2016) compare to our 16.26 mmol O 2 m −2 d −1 of R hy in 2011. Redox-mediated phosphorus desorption from inorganic particles were assumed to be very limited due to adsorption-desorption equilibrium at mineral surfaces (Meng and Yao, 2014) and very high phosphorus burial efficiency  in the study area.

Long-Term Data
Changjiang River discharge (Q) data were downloaded from the website of the Ministry of Water Resources of the People's Republic of China. 2 The discharges were measured at the hydrological gauge station, Datong, which has no tidal influence and is 624 km upstream from the river mouth. Dissolved inorganic nitrogen (DIN) and dissolved inorganic phosphate (DIP) fluxes in August were calculated as the product of Changjiang River discharge in August and nutrient concentrations at Datong gage station or Xuliujing station (Supplementary Figure 1) in historical records (Gu, 1980;Quan et al., 2005;Wang, 2006;Li et al., 2007Li et al., , 2016bDai et al., 2011;Gao et al., 2012) and in situ observations during summers of 2009, 2010, 2011, 2013, 2014, and 2016 at Xuliujing station.
Sea surface temperature (SST) in August along the 31 • N transect was provided by the NOAA-ESRL Physical Sciences Division, Boulder Colorado from their website at https://www. esrl.noaa.gov/psd/, named as COBE SST2.
In situ observations of long-term temperature and salinity in August within the study box were from JODC, 3 KODC, 4 and WOA version 2013. 5 Wind speed data were acquired from the International Comprehensive Ocean-Atmosphere Data Set (ICOADS) with 2degree spatial resolution at the website of the Royal Netherlands Meteorological Institute (KNMI 6 ).

Data Analysis
To quantify the hypoxia severity, we applied PCA to the six observable properties. PCA is a useful technique for reducing the dimensionality of large datasets, increasing interpretability but at the same time minimizing information loss. Numerous variables including those correlated with each other can be synthesized to create new uncorrelated variables (as principal components), for successively maximizing the variance of datasets (Jolliffe and Cadima, 2016). Briefly, a data matrix of six observable properties for eight cruises (Supplementary Table 1) was standardized by substracting the mean and dividing their standard deviation to give a new data matrix, which then was multiplied by the eigenvector to give the variances matrix of six properties for 8 years. We then summed the variances for each year to give the first component. The differences in order of magnitudes of the six properties will not affect the first principal component. The first principal component accounted for 74% of the total variance of the observable properties, and was significantly positively correlated with sectional area (Supplementary Table 1), mean thickness and the longitude span of the hypoxia area, but significantly negatively correlated with DO minimum, mean bottom DO, and mean DO concentration when DO < 100%, representing the status of the "severity" of the hypoxia. Therefore, we use the first principal component as a HI to indicate the severity of hypoxia along the transect. The larger the index, the more severe the hypoxia is.
To estimate the contribution of the direct controlling factors to the HI, multiple regression was applied to generate coefficients of the individual variables, which were then standardized by multiplying by the ratio of the standard deviation of each independent variable to the standard deviation of the dependent variable in order to compare them. We used SPSS R v.13 software for PCA and multiple regression, as well as the outlier identification based on a boxplot procedure (Sim et al., 2005).

DO and Nutrients
The core of the under-saturated water was usually located between 122.5 • E and 123.25 • E (Supplementary Figure 2). The area where DO < 42% (∼3 mg L −1 in DO concentration) varied in shape over the 8 years, and the largest area was in summer of 2016. The lowest DO saturation of 12% was found in 2017 in the bottom water of station B8 (unpublished data).
Dissolved inorganic nitrogen along the 31 • N transect was high in the west and low in the east (Supplementary Figure 3) Figure 2B).
The HI showed that hypoxia has been becoming more severe ( Figure 2C) in terms of decreasing DO concentration and increasing area, thickness and longitude span from 2006 to 2014, but alleviated a bit after 2014.

Direct Drivers of Hypoxia
Seawater density anomaly of the 8 years showed strong stratification along the 31 • N transect (Supplementary Figure 5). The density anomaly difference between bottom and surface ( σ) indicated that strong stratification usually occurred between 122.5 • E and 124 • E. The highest was 13.06 kg m −3 , at 122.65 • E in summer 2009, while the lowest value of 3.68 kg m −3 was from the summer 2006.
τ was the lowest in the 2006, and maximal in 2016 (Supplementary Table 2). Surface NEP in 2013 was the highest and subsurface NEP was the lowest among the eight cruises, indicating a high export production and respiration rates. Temperature associated with oxygen solubility changed little so we did not present it for the eight cruises.

Direct Driver Contribution to Hypoxia Severity
Stratification of the water column plays an important role in hypoxia development . Over the 8 years, the bottom DO was significantly negatively correlated with the density difference ( σ) between bottom and surface (r = −0.66, p < 0.01, n = 47) (Figure 3A), indicating the stronger stratification the lower bottom DO. HI was significantly correlated with averaged σ between 122.3 and 124 • E (r = 0.82, p < 0.01, n = 8) (Figure 3B), demonstrating an important role of stratification in regulating hypoxia.
Hypoxia index also significantly correlated with τ (r = 0.83, p < 0.01, n = 8) (Figure 3C), indicating the longer the water stayed under the MLD, the more severe the hypoxia became, due to insufficient oxygen supply from deep water or the upper mixed layer (Fennel and Testa, 2018). This is consistent with previous studies on river-dominated ocean margins (Rabouille et al., 2008).
In R hy and HI dataset, an outlier (R hy for 2014) was identified. After removing the outlier, the significant correlation between R hy and HI ( Figure 3D) (r = 0.82, p < 0.03, n = 7) showed that community respiration of organic matter by microorganisms such as bacteria (Chen et al., 2006) was a direct crucial factor that controls the extent of hypoxia. Biogenic particles from the euphotic zone dominated the organic matter as a substrate respired in the aphotic layer in the Changjiang plume-impacted shelf area (Wang et al., 2016). In addition, temperature may play a role in determining the respiration rates (Hopkinson and Smith, 2005), but there was no significant correlation between R hy and bottom temperature, perhaps because temperature changed little in the bottom waters over the eight summers (21.29 ± 0.91 • C). If we use the relationship between pelagic respiration and temperature derived by Hopkinson and Smith (2005), respiration would be 391.9 ± 12.9 mmol C m −2 d −1 which varied less than 3% over this narrow temperature range. In fact, respiration may not follow temperature in some estuaries (Jensen et al., 1990;Satta et al., 1996), especially when there is substantial input of external organic matter (Iriarte et al., 1997). While temperature plays an important role in regulating seasonal respiration rates variation, its influence cannot be distinguished during warm summer periods.
According to a multiple regression analysis, σ, τ, and R hy together explained 95.6% of the total variance of HI (r = 0.98, p = 0.02, n = 7). Area and DO min are the most common observable properties used to characterize hypoxic zones. The correlation coefficients between the drivers and area or DO min were lower than for the drivers and HI. The correlations for area and DO min are 0.79 (p = 0.23, n = 8) and 0.75 (p = 0.30, n = 8), respectively, suggesting HI was more robust in describing the severity of hypoxia. The regression equation linking HI to the direct drivers is: The standardized coefficients were 0.06, 0.57, and 0.35, for σ, τ, and R hy , respectively. Based on these coefficients, the three factors contributed 5.7, 55.3, and 34.5%, respectively, to the total variance of HI. If we consider σ and τ as two related physical parameters and combine them into one category, then, we could argue that physical stability of the bottom water contributes more than does respiration in forming hypoxia (61.0 versus 34.5%).

Reconstruction of Hypoxia Index for the Past Six Decades
In this section, we reconstruct the HI during the past 60 years, using Eq. 9 and data about direct drivers from field observations and/or correlations between direct and indirect parameters.
To calculate σ over the long-term, we used averaged salinity and temperature of surface and bottom layers in the study area based on historical data from JODC, KODC, and WOA. σ showed an increase at a rate of 0.025 kg m −3 yr −1 during 1961-1996. After a slight decrease during 1997-2002 it rose after 2003 at a higher rate of 0.047 kg m −3 yr −1 . σ was consistent with the variation of Changjiang River discharge during 1957-2017 (r = 0.32, p < 0.05, n = 55) (Figure 4A), which freshening the sea surface continuously. Previous studies also suggested a long-term freshening in the northern ECS (Siswanto et al., 2008;Park et al., 2011). Although the August discharge temporarily decreased by ∼10% at June 2003 when the Three Gorges Dam began to impound water (Chai et al., 2009), it has been increasing before and after 2003 ( Figure 4A). The more stratified conditions allow more severe hypoxia in the study area (Wei et al., 2007;Zhu et al., 2016).
The long-term τ has been increasing at a rate of 0.11 d yr −1 during 1961-1996, and increased at a much higher rate of 2.52 d yr −1 during 2003-2017 ( Figure 4B). This is also associated with the discharge but more closely associated with salinity of the upper mixed layer (Figure 4B) according to Eq. 6. The fresher the sea surface is, the less water exchange, which allows for longer time for decomposition of organic matter and more consumption of oxygen.
Because respiration data in the hypoxia area are lacking for this long time period, we used NPP based on the significant correlation between NPP and R hy during the 8-year record. As previous studies have found, R h significantly correlated with NPP (Hopkinson and Smith, 2005), so we can use NPP to estimate R hy assuming a fixed averaged hypoxic water thickness over longitude to water depth ratio (0.52) derived from the eight cruises, and assuming 50% contribution of sediment to hypoxia: High NPP was due to higher river DIN flux in recent decades . While summer discharge varied relatively little during the past 60 years (p > 0.05; Figure 4A), nutrient concentrations have increased most of the period (Figures 4C,D) due to anthropogenic activity such as agricultural use of fertilizers and industrial and household pollution (Wang, 2006;Li et al., 2007Li et al., , 2016bDai et al., 2011). After 2010, both DIN and DIP concentrations and fluxes decreased. Together DIN and DIP fluxes accounted for 56.9% of the variation of NPP, whereas individually DIN and DIP contributed 33.3 and 23.6%, respectively. We can build a relationship between NPP and DIN and DIP fluxes (p < 0.05, n = 12): NPP = 4.59 × 10 −8 ± 1.85 × 10 −8 × F DIN + 1.27 × 10 −6 ± 7.25 × 10 −7 ×F DIP + (36.83 ± 7.96) .
Then the long-term R hy ( Figure 4E) can be calculated by incorporating Eq. 11 into Eq. 10. Finally the long-term HI variation can be reconstructed with long-term σ, τ, and R hy using Eq. 9. The reconstructed HI varied between −4.42 and 2.68 (Figure 4F), and the standard deviation was from ±1.95 to ±3.58. The highly significant correlation between the reconstructed HI values and PCA-derived HI values (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) verifies the HI reconstructing model (r = 0.89, p < 0.01, n = 8). The reconstructed HI was highest in 2016 while PCA-derived HI was highest in 2014 ( Figure 4F). This difference could be from the uncertainty of reconstructed HI or the limited sampling locations and periods, suggesting we should be careful when comparing hypoxia severity between specific years, and more historical field data are needed for the verification of HI. However, in summer of 2016 DO was nearly zero north of the 31 • N transect during the same cruise, indicating hypoxia was most severe in 2016, consistent with our reconstructed HI. Therefore, we use reconstructed HI for trend analysis.
FIGURE 4 | Long-term variation of (A) Changjiang River discharge in August and averaged density difference between bottom and surface ( σ); (B) averaged temperature (T) and salinity (S) in layer 1 using historical data and water residence time in layer 2 derived from LOICZ model. (C) DIN (DIN = NO 3 − + NO 2 − + NH 4 + ) concentration and flux and (D) DIP concentration and flux in the fresh water endmember in Changjiang River; (E) averaged NPP derived from satellite along 31 • N transect using the SatCO 2 platform, modeled NPP based on Eq. 11, and modeled respiration in the hypoxia area (R hy ); (F) plane area of bottom hypoxia region (<3 and <2 mg L −1 ) and DO concentration minimum in historical records (Gu, 1980;Tian et al., 1993;Li et al., 2002;Quan et al., 2005;Wei et al., 2007;Pei et al., 2009;Li et al., 2018), compared with reconstructed (recon) and PCA-derived hypoxia index (HI); (G) MEI (multivariate ENSO index) and averaged August SST (COBE-SST2 from NOAA) at 31 • N transect in Changjiang plume-impacted shelf area.
The hypoxia index showed an overall increase over time accelerating in the last two decades ( Figure 4F). Specifically, HI increased slowly during 1961-1996 (with the rate of 0.042 yr −1 ), then it increased during 1997-2002 and during 2003-2017 at higher rates (0.222 and 0.298 yr −1 , respectively). Long-term variation of HI is consistent with the general trends of σ, τ, and R hy (Figures 2E,F, 4B), suggesting all three direct drivers also regulate the long-term HI variation.

Climate Change and Human Activities Effects on Inter-Annual Variation of Hypoxia
ENSO and global warming contributed to long-term variation of hypoxia in terms of fluvial discharge freshening Changjiang plume-impacted shelf area, resulting in stronger stratification and longer residence time in the study area ( Figure 5). Indeed, the average salinity in layer 1, were significantly negatively correlated with σ (r = −0.90, p < 0.01, n = 53) and residence time (r = −0.697, p < 0.01, n = 53). The influence of ENSO on the Asian summer monsoon propagates from winter until early summer (Yang et al., 2018). Warm and cold ENSO phases were represented by MEI, consistent with the surface temperature along the 31 • N transect ( Figure 4G). The warm phase (El Niño events) will elevate freshwater discharge of Changjiang River, through affecting precipitation and evaporation in the Changjiang River drainage area during the summer monsoon (Gong and Ho, 2002). This is supported by the significant positive correlation between MEI (November-December) and August discharge of the next year for Changjiang River during 1955-2017 (Supplementary Figure 6), similar to the results of Bai et al. (2014). Meltwater from glaciers in the Qinghai-Tibetan Plateau in western China, which has increased due to global warming, is another probable freshwater source contributing to the runoff increase (Ding et al., 2006;Siswanto et al., 2008).
Surface warming would also play a role in regulating stratification in Changjiang plume-impacted shelf area (Park et al., 2011); the correlation between SST in layer 1 and σ was significant (r = 0.39, p < 0.01, n = 53), suggesting that global warming and ENSO also contributed to the change in stratification in terms of temperature over the past 60 years (Figure 5). Water masses such as Taiwan Warm Current (Pei et al., 2009), the Kuroshio Current and Ekman transport could also play a role in regulating the long-term variation in stratification (Siswanto et al., 2008), but this needs more data and further exploration.
Moreover, the fourth direct driver-oxygen solubility, was closely associated with subsurface water temperature, and rapid warming in the ECS (Belkin, 2009) will decrease oxygen solubility. Although salinity also affects oxygen solubility, we ignored its effect due to very little variation during the past 60 years (∼0.007 yr −1 ) and ∼1/3 less significant impact on oxygen solubility in seawater than temperature when they change to the same extent. Our analysis showed bottom water has been warmed up by 4.44 • C from 1961 to 2017 [(0.078 • C yr −1 , much faster than SST increasing at a rate of 0.026 • C yr −1 for the whole ECS (Belkin, 2009) and 0.011 • C yr −1 of Kuroshio Current (Wu et al., 2012)]. Since the bottom water mainly comes from Kuroshio Current and will be cooled down and oxygenated each winter, we infer the warming of Kuroshio Current not the bottom water per se contributes to the long-term oxygen solubility decrease. Warming rate of 0.011 • C yr −1 corresponds to oxygen solubility decrease rate of 0.317 µmol L −1 yr −1 at salinity of 34, then increasing rate of HI will be higher by an additional 0.022 yr −1 (based on the correlation between HI and average bottom DO); namely, warming will make hypoxia ∼34% (∼10%) more severe by means of decreasing oxygen solubility during 1961-1996 (1997-2017) in the study area.
If we use multiple linear regression again between the longterm reconstructed HI and the three direct drivers, we can differentiate the contribution of σ, τ, and R hy to the total variance of HI, as 6.1, 63.6, and 30.2%, respectively. Considering the oxygen solubility as the fourth direct driver, the contribution of four direct drivers would be 4. 0, 41.8, 19.8, and 34.4% (5.5, 56.9, 27.0, and 10.6%) during 1961-1996 (1997-2017). Although many studies focus on stratification Zhang et al., 2018) or organic matter-induced respiration Zhu et al., 2011;Wang et al., 2016) based on short time observations, τ played a more important role than the other three drivers in controlling HI variation in Changjiang plume-impacted shelf area during the past decades. This can be attributed to both advection rate and low DO context of Taiwan Warm Current (Qian et al., 2017) or Kuroshio Intermediate Water (Lui et al., 2014). Other studies also addressed the water residence time as a dominant driver of hypoxia, e.g., in Bellingham Bay (Wang and Yang, 2015), Gulf of St. Lawrence (Fennel and Testa, 2018), and comparison study between ECS and Gulf of Mexico (Rabouille et al., 2008), but few quantified its importance.
Wind plays a less important role in regulating interannual stratification, although eastward Ekman transport, forced by the prevailed southerly or southwesterly wind, may contribute to the eastward extension of the Changjiang River plume (Siswanto et al., 2008). If we use the wind mixing parameter h/W 3 (h is the water depth, W is the wind speed) with a value of 0.056 m −2 s 3 (Xuan et al., 2012) and a highest (mean) wind speed of 5.33 (3.41 ± 0.86) m s −1 from 1961 to 2017 along the transect during summer, the MLD is estimated to be 8.47 (2.26 ± 0.035) m. This calculation indicates that wind may influence only the very surface, which is consistent with our 8 years of observations. On the other hand, a typhoon will break the hypoxia, but it would not last long compared with the time scale of weeks to months for hypoxia. Moreover, local hypoxia would persist and even be exacerbated after a typhoon as biological production is enhanced by bring nutrients from bottom to surface water due to the storm mixing (Ni et al., 2016).
Enhanced respiration was mainly caused by anthropogenic activities releasing nutrients to Changjiang plume-impacted shelf area, resulting in elevated NPP (Figures 4F, 5). This is consistent with results from the northern Gulf of Mexico (Rabalais et al., 2014), Chesapeake Bay (Li et al., 2016a), and other eutrophic coastal waters (Breitburg et al., 2018). Increasing temperatures may also lead to higher NPP (Altieri and Gedan, 2015;Li et al., 2016b), but there was no significant correlation between MEI and satellite NPP during 1997-2017, nor was there a significant correlation between SST and NPP. Unlike the eight cruises, the long-term variability in respiration was associated with warming to some extent, as demonstrated by a weak correlation between bottom R hy and temperature (r = 0.376, p < 0.01, n = 53).

CONCLUSION
Overall, in coastal waters or shelf seas, hypoxia severity in the subsurface water is directly regulated by σ, τ, R hy , and oxygen solubility. The year-to-year variation in reconstructed HI indicates that summer hypoxia has become more severe since 1960s, because fluvial nutrients concentrations have increased due to increases in anthropogenic activities, and stratification has been stronger and water residence time longer than before. ENSO and global warming may have indirect affects through regulating the fluvial discharge and altering the stratification and residence time by freshening and warming the surface sea, as well as reducing oxygen solubility in the subsurface water. In addition, higher nutrients inputs, as one of the most significant human activities, will enhance primary production and eventually the respiration rates in the subsurface waters; these nutrient fluxes will also be enhanced by elevated discharge. In conclusion, our analyses suggest that climate change accounts for about 80% of the forces regulating hypoxia in the Changjiang river-impacted shelf while human activity accounts for the remaining 20%.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
KW, JC, and W-JC conceived the study. KW, DH, BW, and WF conducted the field studies. W-JC and DK polished the manuscript. KW processed the data. All authors prepared the manuscript.