The concentrations of seawater phosphate, dissolved silicate, nitrate ( N O 3 − ), and nitrite ( N O 2 − ) were measured using the spectrophotometry method (Hansen and Koroleff, 2007) on a flow injection analyzer (QuAAtro, Seal Analytical, Germany). The detection limits were 0.04, 0.09, 0.02, and 0.005 μmol/L, respectively. Ammonium ( N H 4 + ) was measured immediately upon the retrieval of water samples following a procedure described by Holmes et al. (1999). However, N H 4 + in about 30% of samples was below the detection limit of the method (<0.006 μmol/L). Given our inability to determine the δ¹⁵N of the low-concentration N H 4 + in our samples (averaging 0.20 ± 0.33 μmol/L), we did not subtract the concentration or the δ¹⁵N of N H 4 + from DON and δ¹⁵N of DON. Thus, the DON concentration here was defined as the difference between TDN and N O 3 − + N O 2 − . The concentration of N O 3 − + N O 2 − was quantified using the chemiluminescence method (Braman and Hendrix, 1989) using a N O x − . analyzer (200 EU, Teledyne, USA); the analytical precision is better than ± 0.1 μmol/L. The TDN was first converted into N O 3 − using the persulfate oxidation method (Hansen and Koroleff, 2007). Briefly, 1 ml of persulfate oxidizing reagent (POR, 1.5 g NaOH + 3 g H₃BO₃ + 5 g K₂S₂O₈ in 100 ml low-DON Milli-Q water) was added to 10 ml of seawater sample in a precleaned and precombusted (500°C for 5 h) 15-ml quartz tube with a Teflon-lined cap. The screw caps were closed tightly, followed by autoclaving for 30 min. After digestion, the sample was quantified using the chemiluminescence method. The DON concentration was calculated by subtracting the concentrations of N O 3 − + N O 2 − and the N-blank of POR.

The δ 15 N − N O 3 − , δ 18 O − N O 3 −   and δ¹⁵N-TDN were determined using the denitrifier method (Sigman et al., 2001; Weigand et al., 2016) on an isotope ratio mass spectrometer (PreCon-IRMS, IsoPrime, UK) coupled with a nitrous oxide extraction system. Three reference materials, IAEA-NO-3 ( δ 15 N − N O 3 − = 4.7 %   and   δ 18 O − N O 3 − = 25.8 % ) , USGS-34 ( δ 15 N − N O 3 − = 1.8 %   and     δ 18 O − N O 3 − = 27.9 % ) , and USGS-35 ( δ 15 N − N O 3 − = 2.7 %   and   δ 18 O − N O 3 − = 57.7 % ) , were used for calibration. For data quality assurance, an internal standard of hadal water from the Mariana Trench (8,000 m in depth, δ 15 N − N O 3 − = 5.0 ± 0.2 %   , and δ 18 O − N O 3 − = 2.0 ± 0.2 % was analyzed with each batch of samples. Replicate measurements were performed. The analytical precision was ± 0.2‰ or better for δ¹⁵N and ± 0.5‰ or better for δ¹⁸O. Since both N O 3 − and N O 2 − were reduced to N₂O by the “denitrifier method,” we analyzed the samples for δ 15 N − N O 3 − and δ 18 O − N O 3 − by using sulfamic acid to remove N O 2 − prior to N O 3 − conversion to N₂O (Granger and Sigman, 2009). To evaluate the oxidation efficiency and quantify the isotope deviation caused by the POR-associated N blank, the reference material IAEA-600 (δ¹⁵N = 1.0‰) and the POR were autoclaved and measured with each batch of δ¹⁵N-TDN samples. The δ¹⁵N–DON was calculated from the isotopic mass balance equation (Eqs. 1 and 2).

where the TDN concentration was calculated by subtracting the N-blank concentration of POR.

The SPM was determined as the weight difference between the dried filters and their counterparts before filtration. The C and N content and the corresponding δ¹³C (relative to V-PDB) and δ¹⁵N (relative to air N₂) were measured at the University of California Davis Stable Isotope Facility. The mean standard deviation for reference material replicates was ± 0.06‰ for δ¹³C and ± 0.08‰ for δ¹⁵N.

Optimum multiparameter (OMP) analysis, described by Tomczak and Large (1989), is a water mass inverse model based on linear mixing that is used to find the relative abundances of different source water types (SWT). It finds a solution that best reproduces the observed data and minimizes the residuals in a non-negative least square sense. The OMP was constructed based on two main physically realistic constraints: (i) the study area can be fully represented by the combination of predefined SWT and (ii) the contribution of SWT must be non-negative and add up to 100%. The OMP can be summarized in a matrix form as follows:

where G is the matrix with the property values of SWT, X is the vector of proportional contribution of each SWT to be solved, d is the vector of the observed data, and R is the residual vector.

For the analysis, two conservative (temperature and salinity) and four semi-conservative (oxygen, phosphate, nitrate, and silicate) properties were used to define SWT (Poole and Tomczak, 1999; Budillon et al., 2003; de Carvalho Ferreira and Kerr, 2017). However, as the water masses evolve, the concentrations of these parameters are modified by biogeochemical processes, such as assimilation, remineralization, and respiration processes in the CEECS. The impact of biogeochemical processes can no longer be disregarded and have to be included in the analysis. The extended OMP analysis was applied to minimize these effects by adding the Redfield ratio item supplied by OMP files (a detailed description is available in Supplementary Information Text S2 ). The SWT parameter values (a detailed description is available in Supplementary Information Text S3 and Supplementary Information Text S4 ) used are summarized in Table 1 .

Sensitivity tests were performed to test the impact of uncertainties on the results. First, we simply changed the weights. The weights were from literature (de Carvalho Ferreira and Kerr, 2017; Zhou et al., 2018), experience, OMP files, and calculation (a detailed description is available in Supplementary Information Text S2 ). The best results were obtained by using the calculated weights ( Table 1 ), which not only effectively represent the varieties in the formation region but also hold the minor mass balance residuals. All mass balance residuals were lower than 8%, only four residuals were higher than 5%, and 74% of the residuals were lower than 1% ( Supplementary Figure S2 ). According to the mass balance residuals, the results of the OMP analysis were reliable (Tomczak and Large, 1989). Furthermore, we analyzed the same dataset by adding and subtracting the standard deviation based on the observed variability ( Table 1 ) from the properties of the SWT (matrix G) to assess the robustness of the solution. Most of the standard deviations of the solutions were <10%. Thus, the uncertainties in the OMP runs were considered acceptable (Gasparin et al., 2014; Zhou et al., 2018).

Based on the relative abundances of different SWT, the expected DON concentration and δ¹⁵N–DON affected only by the physical mixing can be calculated as follows:

We further calculated the deviation of observed values from calculated values. The deviation is defined in terms of Error-DON.

The steady-state model is frequently used to interpret the N isotope data from the ocean (Sigman et al., 2009).

For the reactant N pool:

For the product N pool:

where f is the fraction of the reactant remaining.

The OMP analysis was used to determine the relative contribution of the core water masses in the CEECS ( Supplementary Figure S3 ). CDW flowed in the east, southeast, and especially northeast directions in the surface layer, which is consistent with previous studies (Liu et al., 2021). In the surface water of the northern study area (transect A1–A4), there is a weak bidirectional northwestward extrusion of CDW over Subei Shoal. From there, the contribution of CDW gradually decreases with increasing distance from the coastal area. YSMW and KSW are present as high contributions in the northern nearshore area and offshore area, respectively. In the southern study area, KSW contributes more than 50% of the total water in the entire surface layer. The remaining contribution comes mainly from TWC and CDW in the nearshore area. The bottom water in the study region is a mixture of TWC, KBBCNT, and YSMW. KBBCNT contributes the most. YSMW is confined to the northern part of the study area.

The difference between measured values and calculated values (Eq. 4) based on OMP results was defined as Error-DON. When the measured values are equal to the calculated values, it means that the physical mixing of the water masses is the dominant factor controlling the spatial distribution of DON. A positive Error-DON corresponds to the net production of DON, and a negative Error-DON corresponds to the net consumption. To identify the DON dynamics in different water masses, Error-DON versus measured DON was plotted ( Figure 6A ). We found that both net production and net consumption of DON occurred in every water mass. In the study area, the Error-DON averaged 10.7 ± 8.6% for DON net production and −9.0 ± 7.0% for net consumption. Interestingly, in all the sampling sites (n = 138), there were equal numbers of DON net production (n = 69) and consumption sites (n = 69). This result could explain why there is a misleading impression that DON exhibited a conservative behavior in general over the CEECS (from the river mouth to the offshore area up to 800 km; Kwon et al., 2018). In actuality, we captured the non-conservative behavior (active biochemical cycling) of DON via the OMP analysis.

The positive Error-DON values were distributed mainly in the southern coastal areas and northern part of the study area (see the diamonds in Figure 6B ). In these DON net production sites, the relationship between δ¹⁵N–DON and DON concentration was categorized into two patterns. In the N O 3 − rich euphotic layer (9.1 ± 9.2 μmol/L, P-zone 1), the DON correlated negatively with δ¹⁵N–DON (R ² = 0.57, p ≤ 0.001, n = 29, Figure 7A ), while in the coastal middle and bottom waters where the SPM was high (64.0 ± 90.1 mg/L, P-zone 2), the DON correlated positively with δ¹⁵N–DON (R ² = 0.29, p ≤ 0.001, n = 37, Figure 7B ). The two patterns indicated that there were two different mechanisms of DON production in the two zones. In P-zone 1, the δ¹⁵N of the produced DON decreased the δ¹⁵N–DON in the water column. Conversely, the δ¹⁵N of the produced DON elevated the δ¹⁵N–DON in the water column in P-zone 2. A detailed explanation is provided in Section 4.2.

The negative Error-DON values were concentrated mainly in the southern part of the study area and part of the northern bottom layer (see dots in Figure 6B ). Based on the distribution of net consumption sites and the apparent oxygen utilization (AOU), we classified the net consumption sites into two zones: (i) C-zone 1, where the AOU values were negative (AOU = −0.12 ± 0.68 mg/L, n = 17), and (i) C-zone 2, where the AOU values were positive (AOU = 3.17 ± 0.69 mg/L, n = 36). A negative correlation was observed between the AOU and DON concentration (R ² = 0.28, p = 0.001, n = 36, Supplementary Figure S5A ). With the assistance of multiple N-isotopic signatures, we determined the most plausible turnover processes for DON in each zone.

DON can be released through the autochthonous biological process of extracellular exudate production by phytoplankton (Bronk, 2002). The surface photosynthetic sources for DON in the Eastern Tropical South Pacific euphotic zone (Knapp et al., 2018) and the South China Sea surface waters (Zhang et al., 2020) were suggested through a positive correlation between the DON stock and chlorophyll-a at these sites. Unfortunately, we lack chlorophyll-a data. However, the N O 3 − dual-isotope variation pattern in P-zone 1 is capable of characterizing the phytoplankton activity. In P-zone 1, the observed N O 3 − concentrations (9.1 ± 9.2 μmol/L) were lower than the expected values (15.7 ± 10.7 μmol/L), while both δ 15 N − N O 3 − ( 12.1 ± 3.8 % ) and δ 18 O − N O 3 − ( 8.9 ± 4.1 % ) were heavier than the expected isotope ratios δ 15 N − N O 3 − : 6.0 ± 0.04 % and δ 18 O − N O 3 − : 2.8 ± 0.7 % In other words, N O 3 − was partially consumed accompanied by the enrichment of δ 15 N − N O 3 − and δ 18 O − N O 3 − in P-zone 1. This suggests the dominance of N O 3 − assimilation by phytoplankton (Zhong et al., 2020). The co-variation of N O 3 − dual isotopes fits a 1:1 line (R ² = 0.96, p ≤ 0.001, n = 29, Supplementary Figure S6A ), as also expected for N O 3 − assimilation by phytoplankton (Granger et al., 2004). A positive correlation between N O 3 − consumption and DON concentration (R ² = 0.58, p ≤ 0.001, n = 29, Supplementary Figure S6B ) suggested the net DON production at these sites as well. In N O 3 − rich surface waters like P-zone 1, phytoplankton preferentially consume ¹⁴N relative to ¹⁵N (Sigman and Fripiat, 2019). The isotope effect for N O 3 − assimilation produces photosynthetic biomass (which, in turn, produces DON) with a relatively low δ¹⁵N (Knapp et al., 2018). Thus, the phytoplankton-sourced DON may hold a relatively low δ¹⁵N value in P-zone 1.

The isotope effect is a key parameter that associates N O 3 − assimilation to the ¹⁵N/¹⁴N of N O 3 − and DON. Given the fact that N O 3 − (the reactant) can be continuously replenished by nitrification (Zhong et al., 2020), a steady-state model (Eq. 7, applied to the product N pool) is more appropriate to estimate the isotope effect (¹⁵ϵ_phy) during DON production in P-zone 1. The slope of the linear regression between the δ¹⁵N–DON and the fraction of N O 3 − remaining in the water column after N O 3 − assimilation yielded an estimate for ¹⁵ϵ_phy of 2.6 ± 0.6‰ ( Figure 6C ). The δ¹⁵N_initial (y-intercept) was 3.8 ± 0.4‰, which was consistent with the average value (3.8 ± 1.0‰) of the expected δ¹⁵N values (theoretical value under physical mixing; Eq. 5) for DON in P-zone 1. This result validated the application of the steady-state model, and the addition of phytoplankton-sourced δ¹⁵N–DON decreased the δ¹⁵N–DON in the water column.

Approximately 87% of the annual particulate load from Changjiang River is discharged in summer (Zhu et al., 2011). The P-zone 2 is a high-SPM (64.0 ± 90.1 mg/L) region. High SPM reduces light penetration in the water column and suppresses the photosynthesis process and primary production (He et al., 2017). From the perspective of δ 15 N − N O 3 − , the δ 15 N − N O 3 − (average of 6.6 ± 0.9‰) exhibited a homogeneous behavior in P-zone 2, suggesting that it was unlikely to be an occurrence of intensive phytoplankton assimilation. Thus, the phytoplankton-sourced autochthonous DON was not the dominating source of DON in P-zone 2, even though nutrients were rich in this area. PN and DON are the two main organic N species in the water column. In CEECS, PN was significantly controlled by SPM transport (R ² = 0.98, p ≤ 0.001, n = 40, Supplementary Figure S7 ). SPM (PN) transport and estuarine physical mixing are two independent processes (Gao et al., 2020). More than half of the fine‐grained particles are temporarily deposited near the estuary in P-zone 2 (Chen et al., 2017). Thus, there are possibilities that PN and DON pools interacted during transportation in P-zone 2. An inter-transformation between the two most important organic N pools may have occurred (Gebhardt et al., 2005) and altered the distribution and variation of DON observed in this area.

PN produces DON mainly through cell lysis, exudation, and particle solubilization (Bronk and Steinberg, 2008). These processes are more likely to occur with the rupture of the C–C bond but are much less likely to occur with the rupture of the C–N bond. Thus, the process of producing DON from PN would not cause a significant isotope fractionation, and the produced δ¹⁵N–DON should be similar to the δ¹⁵N of PN from which DON is derived (Knapp et al., 2011). The δ¹⁵N–PN ranged from 0.7 to 7.1‰, with an average of 3.2 ± 1.8 ‰. This was slightly higher than that of δ¹⁵N–DON (0.1 to 5.5‰, average 2.9 ± 1.6 ‰) in P-zone 2. The mean values of the δ¹⁵N–PN and δ¹⁵N-DON in P-zone 2 showed a statistically significant difference (ANOVA, p < 0.01). This difference explains why the DON in P-zone 2 correlates positively with δ¹⁵N–DON (R ² = 0.29, p ≤ 0.001, n = 40, Supplementary Figure S4B ). The addition of PN-sourced DON elevated the δ¹⁵N–DON in the water column.

DON is known as a potential source of N-nutrient for phytoplankton (Knapp et al., 2018). Hu et al. (2012) demonstrated the ability of Prorocentrum donghaiense, which frequently occurs in CEECS, to assimilate DON. Moreover, Zhang et al. (2015) suggested that DON is bioavailable during the diatom to dinoflagellate bloom succession. Thus, in C-zone 1 where N O 3 − is completely consumed, DON may play an important role in fueling phytoplankton growth. The bio-active N released from DON in C-zone 1 likely occurs with isotope fractionation when hydrolysis reactions break common C–N bonds (Knapp et al., 2011). Since DON can be continuously supplied by water mixing in the horizontal direction, a steady-state model (Eq. 6 applied to the reactant N pool) would fit the continuous DON supply and removal. The isotope effect (¹⁵ϵ_ass-DON) during DON consumption in C-zone 1 was estimated to be 3.7 ± 1.2‰ ( Figure 7A ), which is in accordance with the isotope effect (3–10‰) of the conversion of amide to amine and amine to ammonia (O′Leary and Kluetz, 1972; Macko et al., 1986; Silfer et al., 1992). The result was slightly lower than the isotope effect of 4.9 ± 0.4‰ during DON consumption in the upper layers of South China Sea (Zhang et al., 2020) and 5.5 ± 1.2‰ in the upper 50 m of Eastern Tropical South Pacific (Knapp et al., 2018). This could result from the residence time of DON that determines the extent of DON consumption. Previous studies suggested that the residence time ranged from days to months in C-zone 1 (Gu et al., 2012; Tan et al., 2018; Wang et al., 2018). In other words, DON has a shorter residence time in this area than in the upper layers of South China Sea (∼3 years; Liu and Gan, 2017) or the upper ocean of Eastern Tropical South Pacific (months to years; Knapp et al., 2011). Hopkinson et al. (2002) showed that the average half-lives of very labile and labile DON are 12 and 113 days, respectively. We infer that the labile DON is not completely consumed due to the short residence time in C-zone 1. Therefore, the isotope fractionation of labile DON consumption in C-zone 1 is weaker than that of refractory DON consumption.

In C-zone 2, the enrichment ratio of ¹⁸O and ¹⁵N from the potential N O 3 − sources deviated from the slope of 1.0, spreading over the range from 1.0 to 2.0 ( Supplementary Figure S8A , Granger et al., 2004). This suggests that continuous nitrification occurred in this area (Umezawa et al., 2014; Wang et al., 2016). The negative correlation between the AOU and DON concentrations (R ² = 0.28, p = 0.001, n = 36, Supplementary Figure S5A ) and the significantly positive correlation between AOU and N O 3 − concentrations (R ² = 0.64, p < 0.001, n = 36, Supplementary Figure S5B ) suggest the aerobic microbial degradation of DON to be associated with N O 3 − regeneration and a coupling between the consumption of DON pools and the production of N O 3 − in C-zone 2. The DON remineralization and nitrification elevate the NO ₃ ⁻ concentration (Liu et al., 2020).

We use the steady-state model (Eq. 6 applied to the reactant N pool) to estimate the isotope effect during DON consumption. The slope of the linear regression yielded an estimate for ¹⁵ϵ_ntr-DON of 5.6 ± 2.6‰ ( Figure 7B ). The ¹⁵ϵ_ntr-DON was consistent with the isotope effect of 5.5 ± 1.2‰ for the DON degradation estimated from data for the upper 50 m of Eastern Tropical South Pacific, where the residence time was many months to years (Knapp et al., 2018). In C-zone 2, KBBCNT makes the highest contribution (>50%, Supplementary Figure S3 ). The shelf water exchange with Kuroshio Current is 1–2.3 years [summarized in Chen (1996)], which indicates that DON in C-zone 2 undergoes a longer period of degradation than in C-zone 1. As a result, DON degrades more thoroughly, and the remaining DON is more refractory. The δ¹⁵N_initial (y-intercept) was 4.2 ± 0.3‰, which is slightly lower than the average value (4.7 ± 0.2‰) of the expected δ¹⁵N values (theoretical value derived from physical mixing; Eq. 5) for DON in C-zone 2. The 0.5‰ difference may result from DON originating from “fresh” PN in C-zone 2. The isotope δ¹³C is a useful proxy to discriminate the sources of POC. In C-zone 2, δ¹³C averaged −21.9 ± 1.1‰, suggesting that marine phytoplankton was the predominant origin of POC (−22 to −18‰; Middelburg and Nieuwenhuize, 1998). The linear relationship between POC and PN (R ² = 0.91, p ≤ 0.001, n = 36, Supplementary Figure S8B ) suggests that PN was strongly associated with POC. The C/N molar ratio between POC and PN was 5.5 ± 0.3, which is lower than the Redfield ratio of 6.6 (Redfield, 1960), suggesting that the PN was “fresh” (Liu et al., 2019). The solubilization of PN is an important source of DON (Benner and Amon, 2015). The δ¹⁵N–PN ranged from 0.6 to 6.1‰, with an average of 3.3 ± 1.6 ‰, which was lighter than that of δ¹⁵N–DON (3.0 to 6.7‰, averaged 4.7 ± 0.9 ‰) in the C-zone 2. The mean values of δ¹⁵N–PN and δ¹⁵N-DON in the C-zone 2 showed a statistically significant difference (ANOVA, p < 0.01). As mentioned in Section 4.2.2, the δ¹⁵N of DON produced was similar to the δ¹⁵N–PN. Thus, the addition of DON from PN would decrease the expected δ¹⁵N–DON in C-zone 2. Therefore, we suggested that the expected δ¹⁵N–DON had already been substantially modified by DON originating from “fresh” PN. In this section, we estimated the isotope effect during DON consumption and clarified the confounding coupling of N O 3 − , DON, and PN in the N cycling in C-zone 2.

The remaining sites not included in these four zones are concentrated in coastal surface waters (n = 17, Figure 6B , gray dots), where the Error-DON was relatively low (−3.3 ± 3.5%). Substantial challenges still exist for studying DON cycling in this zone. On the one hand, it is still uncertain whether this deviation is caused by biological processes or uncertainties for calculation. On the other hand, both the N O 3 − and PN concentrations were high in this zone, and the variation caused by DON consumption was relatively small. Understanding of the particularly complex and overlapping interactions between N O 3 − , PN, and DON remains elusive. Further studies are required to investigate this issue.