Nutrient co‐limitation in the subtropical Northwest Pacific

Nutrients limiting phytoplankton growth in the ocean are a critical control on ocean productivity and can underpin predicted responses to climate change. The extensive western subtropical North Pacific is assumed to be under strong nitrogen limitation, but this is not well supported by experimental evidence. Here, we report the results of 14 factorial nitrogen–phosphorus–iron addition experiments through the Philippine Sea, which demonstrate a gradient from nitrogen limitation in the north to nitrogen–iron co‐limitation in the south. While nitrogen limited sites responded weakly to nutrient supply, co‐limited sites bloomed with up to ~60‐fold increases in chlorophyll a biomass that was dominated by initially undetectable diatoms. The transition in limiting nutrients and phytoplankton growth capacity was driven by a gradient in deep water nutrient supply, which was undetectable in surface concentration fields. We hypothesize that this large‐scale phytoplankton response gradient is both climate sensitive and potentially important for regulating the distribution of predatory fish.

Surface waters of the western subtropical North Pacific have chlorophyll a (Chl a) and nutrient concentrations that are among the lowest globally Longhurst 2010). The mean state of this ecosystem is characterized by a dominance of small picophytoplankton (<2 μm), with tightly balanced rates of growth and mortality (i.e., via grazing, viral lysis). Efficient remineralization of cell-bound nutrients explains how picophytoplankton can maintain elevated growth rates in an environment where the standing stocks of nutrients are low (Goldman et al. 1979). Less clear, however, is how these types of system can sustain export of sinking particulate organic matter or host large predatory fish (Lehodey et al. 1997;Longhurst 2010;Tréguer et al. 2018;Landry et al. 2019). Both processes (i) require external supply of new nutrients, and (ii) are enhanced by a shift in phytoplankton community structure from picophytoplankton to larger diatoms (Irigoien et al. 2002;Tréguer et al. 2018). The latter increases phytoplankton cell size by tens of thousands (by volume), which both enhances particle sinking and shortcuts several trophic transfer stages, increasing energy flow to higher predators (Irigoien et al. 2002).
In the subtropical Northwest Pacific, fixed nitrogen (N) supply is generally assumed to be the main limiting factor for phytoplankton (Moore et al. 2004;Hashihama et al. 2009;Longhurst 2010;Moore et al. 2013;Li et al. 2015). In this case, increased phytoplankton growth rates and a shift from picophytoplankton to diatoms would depend primarily on the enhanced supply of N. However, a recent data synthesis has highlighted that, alongside N, concentrations of phosphate (P) can also be strongly depleted in this region (Martiny et al. 2019). Additionally, further south at the western extent of the Equatorial Pacific upwelling zone, measurements of phytoplankton photophysiology and the results of a nutrient addition experiment both independently pointed toward a system approaching co-limitation by both N and iron (Fe) (Behrenfeld et al. 2006;Li et al. 2015). While these different lines of evidence argue that Fe and/or P could be playing important roles as (co-)limiting nutrients alongside N in the subtropical Northwest Pacific, there have been no studies systematically testing for this . Conducting replicated experiments with the required factorial nutrient amendment design is challenging, but affords the advantages of (i) defining the identity of (co-)limiting nutrient(s) and (ii) quantifying phytoplankton responses to individual and combinations of nutrient supply in a way that can be mapped at ocean scales (Moore et al. 2008;Browning et al. 2017). Here, we present the results of 14 such experiments conducted along an~8000 km cruise track throughout the subtropical Northwest Pacific.

Methods
Experiments were conducted onboard the RV Tan Kah Kee from 25th April 2019 to 13th June 2019 (KK1903). Data are available in Browning (2021). Experiments directly followed previously published protocols (Browning et al. 2017). Following collection of whole, unfiltered seawater (1 L acid-washed polycarbonate bottles; Nalgene) with a towed sampling device (~2 m depth; Zhang et al. 2019), seawater was spiked in triplicate with N (1 μM NO 3 À + 1 μM NH 4 + ), P (0.2 μM) and Fe (2 nM) in factorial combinations. At several of the sites, additional treatments also conducted were NO 3 À (2 μM) alone, and NO 3 À + Fe (2 μM and 2 nM, respectively). Samples were incubated for~48 h in on-deck incubators (Lee Filters "Blue Lagoon" screening and flushed with surface seawater). Samples for Chl a were filtered (250-500 mL), extracted in 90% acetone, and analyzed fluorometrically on a 10AU Turner Designs fluorometer (Welschmeyer 1994). Flow cytometry samples were preserved at 1% paraformaldehyde concentration, frozen at À80 C, and later analyzed using a CytoFLEX (Beckman coulter; Marie et al. 2000). Phytoplankton groups were identified and enumerated with FlowJo software. Fast Repetition Rate fluorometry measurements were undertaken using a FastOcean with integrated Act2 laboratory system (Chelsea Technologies group) yielding minimum and maximum fluorescence (F o and F m , respectively) that were blank corrected using fluorescence of 0.2 μm filtrates. Samples for diagnostic pigment analyses (initial t = 0 samples 4.4-5.3 L; pooled seawater from replicates for some treatments~2 L) were filtered onto glass fiber filters, frozen at À80 C, then later analyzed using a Shimadzu LC-20A HPLC system (Heukelem and Thomas 2001). Macronutrient and dissolved Fe samples were collected filtered (0.8/0.2 μm pore size; AcroPak, Pall). Analyses followed previously described protocols: phosphate was analyzed on ship (Ma et al. 2008); nitrate + nitrite were preserved (À20 C) and determined following Zhang (2000); and ammonium was analyzed on-ship (Zhu et al. 2013). Dissolved Fe samples were acidified and stored for >6 months and analyzed following Rapp et al. (2017) except that a Preplab instrument (PS Analytical) was used for preconcentration and standard addition was used for concentration determination. Analysis of GEOTRACES intercalibration standard GSP91 yielded Fe concentrations similar to determinations by other laboratories (0.13 nM in the same analytical run as the samples, n = 2; 0.16 AE 0.045 nM over multiple analytical runs, n = 8; Wuttig et al. 2019) (https://www.geotraces. org/standards-and-reference-materials/). Vertical profiles of dissolved Fe were collected using a trace-metal-clean rosette with Niskin-X bottles and analyzed by chemiluminescence-based flow injection (King et al. 1995) after preconcentration by a PA-1 column.
Rates of N 2 fixation were determined using the 15 N 2 gas dissolution method (Mohr et al. 2010), conducted in duplicate 4.5 L Nalgene polycarbonate bottles. The 15 N 2 pre-dissolved seawater (2 L) was prepared with 15 N 2 gas (98.9 atom %, Cambridge Isotope Laboratories) following Lu et al. (2018). Samples were spiked with 100 mL 15 N 2 enriched filtered seawater prepared at the same site and incubated on-deck for 24 h. The final 15 N 2 enriched seawater concentration in the incubation bottles (~3 atom %) was calculated based on an assumed 100 atom % concentration in the pre-dissolved seawater. Measurement of enriched seawater prepared in an identical manner by the same investigators on a previous cruise via GasBench-IRMS confirmed this prepared seawater was always >95 atom %; n = 6; Wen et al. 2017). The exact enrichment of individual incubation bottles was not determined for this research cruise. Incubated samples were filtered onto pre-combusted (450 C, 4 h) GF/F filters. Surface seawater particulate organic matter was also collected to determine the natural, background 15 N-PON abundance. Samples were analyzed using an elemental analyzer coupled to a mass spectrometer (EA-IRMS, vario PYRO cube-Isoprime 100). The N 2 fixation rates were then calculated according to Montoya et al. (1996), taking 4‰ as the minimum acceptable change in the δ15N of particulate nitrogen.
To estimate nitrate transfer into the euphotic zone via turbulent diffusion, the nitrate concentration gradient at the 1% light depth was calculated and multiplied by a turbulent diffusion coefficient of 1 Â 10 À5 m 2 s À1 (Lee et al. 2007;Kaneko et al. 2021). This was undertaken for both depth profiles of nitrate determined on the research cruise and for the World Ocean Atlas annual average climatology.

Results and discussion
All investigated sites were strongly and equally oligotrophic (Chl a 0.01-0.04 mg m À3 ) and depleted in fixed nitrogen (N; nitrate + nitrite ≤10 nM; ammonium <15 nM; Table 1; Fig. 1). In line with the extreme N depletion, experimental stimulation of any increase in Chl a required the supply of at least N ( Fig. 1; one-way ANOVA, α = 0.05, followed by a Tukey posthoc test). The default N treatment was nitrate plus ammonium, but additional further treatments at several sites with nitrate alone showed similar results ( Fig. 1c; Supporting Information S1). At Sites 1-7, located in the northern part of the study region, significant increases in Chl a were observed following supply of N alone. Within this zone, further increases in Chl a were observed following combined supply of N + P at Sites 1, 2, and 4; here, we define this type of response as serially N-P limited.
Further south at Site 7, supplying N alone also increased Chl a, however, the serial limiting nutrient switched identity from P to Fe. Furthest south at Sites 8-13, Chl a increases following N-only supply were small and not significant with the applied statistical test. Chl a responses at these sites instead only became significantly enhanced following the combined supply of N + Fe, indicating these nutrients were co-limiting. Finally, Sites 5-6 and 14, located geographically in between the serially N-P limited sites to the north and the N-Fe colimited sites to the south, exhibited primary N limitation followed by serial responses to N + P + Fe supply (i.e., primary N limitation followed by serial P + Fe co-limitation).
The region of N and Fe co-limitation identified by the Chl a responses in the southern portion of the study region was further distinguished in a number of cases by distinct photophysiological changes following N supply (Figs. 1c and S1; Behrenfeld and Milligan 2013). At N (and serial N-P) limited sites to the north, the apparent photosystem II photochemical efficiency parameter, F v /F m , was high and generally varied insignificantly following supply of any nutrient combination (mean = 0.45, sd = 0.06, n = 222). Values of F v /F m at N-Fe serially/co-limited sites 7-11 were also elevated (mean = 0.47, sd = 0.05, n = 93) and remained so following Fe addition in any treatment combination, but in contrast to the northern sites, often declined following N supply (either as +N or +N + P; mean = 0.37, sd = 0.03, n = 33). These F v /F m reductions were a result of tipping the N-Fe co-limited communities into stronger Fe limitation following the artificial supply of N (Behrenfeld et al. 2006;Behrenfeld and Milligan 2013;Browning et al. 2017).
In addition to the qualitative classification of the sites into different nutrient limitation regimes, the magnitudes of Chl a responses to supply of limiting nutrients between experimental sites also showed strong geographic variability (Table 1). This variability was coherent with the identity of the limiting nutrients. Chl a based net growth rates were 0.8-0.9 d À1 (mean 0.9 d À1 ) at N (and serial P) limited sites following supply of all three nutrients in combination (N + P + Fe; Sites 2-4). In contrast, this full nutrient treatment resulted in net growth rates of 1.4-2.1 d À1 (mean 1.8 d À1 ) at N-Fe co-limited sites (Sites 8-12). Corresponding net growth rates at serially N-Fe limited sites fell in between these ranges (1.2-1.8 d À1 , mean 1.5 d À1 ). Therefore, in addition to arranging into geographically coherent categories in terms of the identity of limiting nutrients, overall magnitudes and rates of Chl a biomass responses generally matched such categorization, with two-fold higher mean net growth rates following nutrient supply at the N-Fe co-limited sites in the south than at the N (or serially N-P) limited sites to the north.
This marked geographic variability in Chl a growth responses following the supply of limiting nutrient(s) was strongly tied to changes in phytoplankton community structure (Fig. 2). At N (and serially N-P) limited sites in the northern part of the study region, Chl a responses following supply of limiting nutrients were mostly attributable to increases in numbers and intracellular Chl a concentration of phytoplankton communities that dominated the initial seawater, that is, cyanobacteria and photosynthetic picoeukaryotes (Fig. 2a,b). A combination of low maximum specific growth rates of these phytoplankton and their rapid consumption by small, fast responding grazers were probably key in restricting rates and magnitudes of observed biomass enhancements (Landry et al. 2000;Moore et al. 2008;Browning et al. 2017; net growth rates derived from cell concentrations following N + P treatment at N/serially N-P limited Sites 2-3 ranged 0.57-1.04 d À1 (Prochlorococcus); 0.59-0.8 d À1 (Synechococcus); 0.23- Table 1. Initial conditions, identified limiting nutrients, and responses to nutrient supply in experiments.   0.61 d À1 (picoeukaryotic phytoplankton)). In contrast, the N-Fe co-limited sites to the south exhibited extreme shifts in phytoplankton community structure following N + Fe supply, with diatoms increasing from undetectable contributions to an estimated 68-83% of total Chl a ( Fig. 2b; net growth rates of the diatom-associated pigment fucoxanthin were 2.55-3.32 d À1 ; Fig. S2). Synechococcus and photosynthetic picoeukaryote concentrations and cellular Chl a content also demonstrated evidence for N-Fe serial/co-limitation at these sites ( Fig. 2a; net growth rates derived from cell concentrations at N-Fe co-limited Sites 8-13 ranged 0.12-0.52 d À1 for +N and 0.86-1.19 d À1 for +N + Fe (Synechococcus); and 0.19-0.52 d À1 for +N and 0.43-0.78 d À1 for +N + Fe (picoeukaryotic phytoplankton)); however, unlike the picophytoplankton, expected higher specific growth rates of diatoms together with their larger, slower responding grazers, are suggested to have enabled higher magnitude biomass accumulation following N + Fe supply (Landry et al. 2000). Accordingly, in addition to overall  . Upper panels are derived from cell concentrations, lower panels fluorescence per cell (F cell À1 , a proxy for intracellular Chl a content). PPE, photosynthetic picoeukaryotes. (b) Phytoplankton community responses for select treatments derived from diagnostic pigment analysis and converted to approximate contributions of each phytoplankton type to total Chl a. Apart from experimental sites 1 and 2, the estimated fractional contribution of diatoms in initial samples was always zero. Pigment concentrations were converted to approximate contributions of different phytoplankton types using CHEMTAX (Mackey et al. 1996)  biomass yields, mean net Chl a growth rates in experimental treatments across the dataset were highly predictable by final diatom contributions to the phytoplankton community (R 2 = 0.80, p < 1 Â 10 À16 , n = 48).
What regulated the marked north-south transition in biological responses to nutrient supply? These sites were indistinguishable in terms of initial N, Chl a, and community structure. Inspection of satellite-derived current velocities and drifting float trajectories showed that this transition instead mapped onto changes in near-surface currents (Fig. 3a,b). While Sites 2-6 in the north were located within currents having no sustained flow direction, Sites 8-12 to the south were situated in coherent westward flowing waters of the North Equatorial Current (NEC) that originate in the eastern tropical Pacific (Talley et al. 2011). No gradient in surface nitrate concentration exists across the transition between these ocean currents (Table 1); however, climatological nitrate datasets show that Ekman divergence as a result of westward acting wind stress in the NEC zone shoals the nitracline depth by~50 m (here defining the nitracline as the 1 μM nitrate contour; Fig. S3).
Using nitrate gradients calculated at the euphotic depth alongside an assumed representative value of turbulent diffusivity, we calculated the regional distribution of the turbulent diffusive nitrate flux into the euphotic zone (Figs. 3c and S3). In contrast to surface nitrate concentrations, which are ubiquitously low, mean nitrate transport upwards into the euphotic zone was estimated to be around five-fold higher in the NEC zone (N-Fe co-limited sites 9-13) in comparison to further north (N limited sites 2-6) (Figs. 3c; also see S4). Depth profiles of dissolved Fe concentrations demonstrated that, unlike nitrate and phosphate, waters below the mixed layer were not substantially enriched in Fe, and were furthermore indistinguishable between the NEC zone in the south and further to the north (Fig. S5). Consequently, entrainment of lower Fe:N waters around the NEC via Ekman divergence and nitracline shoaling, alongside low expected aerosol Fe deposition rates in this zone (Fig. S6), would have both together contributed to the development of the observed N-Fe co-limited conditions. Spatial patterns in the relative supply rates of N and Fe also offer a mechanistic framework for understanding the responses to experimental supply of P ( Fig. S7; Table 1). Following N supply, serial Chl a responses to P addition were observed in combination with Fe (i.e., serial P-Fe co-limitation) at Sites 5, 6, and 14 and alone at Sites 1, 2, and 4 (all sites away from the NEC zone in the south). A combination of (i) reduced upwelled N away from the NEC, together with (ii) a probable increase in Fe supply nearer Asian landmasses (Fig. S6), would together increase the overall Fe:N supply ratio to surface waters, creating a niche for diazotrophs (N 2 -fixing bacteria; Ward et al. 2013). In line with this theoretical expectation, mean N 2 fixation rates measured on the research cruise were >20-fold higher at sites away from NEC waters (Table 1; Fig. 3d; also see Hashihama et al. 2009;Kitajima et al. 2009). The enhanced N 2 fixation at these sites would introduce new N into the system without a corresponding addition of P, which over time would drive the observed P drawdown in surface waters (Table 1), leading to the serial P responses observed ( Fig. S7; Supporting Information S2; Wu et al. 2000;Martiny et al. 2019).
In addition to establishing a north-south spatial gradient in the nitracline, the westward advecting NEC current also transports phytoplankton zonally across the tropical Pacific (Villarino et al. 2018). This transport could retain signatures of the community structure from more productive waters upstream (Villarino et al. 2018), which typically contain more diatoms. Therefore, although not detected by diagnostic pigments (Fig. 2b), we hypothesize a westward transport of bloom forming diatoms and/or their resting spores in the NEC from the eastern subtropical Pacific. Further diatom input could result from enhanced vertical exchange with submixed layer waters (Fig. 3a;Furuya 1990). Whether laterally or vertically, any input term of bloom-forming diatoms into the surface layer would effectively reduce their equilibrium resource concentration, R* (Lévy et al. 2014), which based on the extremely depleted nutrient concentrations would be expected to become fully outcompeted by smaller picophytoplankton at all sites (Dutkiewicz et al. 2009). Together, we hypothesize that the enhanced nitrate, and potentially diatom, supply rates into the surface layer in the NEC might be key to retaining the low-level diatom seed population, which appears fundamental in regulating the subsequent rapid biological growth response to supply of limiting nutrients (Figs. 1c and 2b; Cermeño et al. 2016). In contrast to the NEC zone, weaker circulating currents further north would increase the time for competitive exclusion of bloom forming diatoms by picophytoplankton (Figs. 2b and 3a,b;Lévy et al. 2014), leading to the much weaker growth responses to nutrient supply observed there (Fig. 1c).
The identified large-scale gradient in rapid diatom growth capacity could be important for regulating predatory fish (Fig. 3e). In addition to greater enhancements of primary producer biomass (Fig. 1c), growth of larger diatoms shortcuts two or more trophic positions in the food chain between primary producers and predatory fish (that is, heterotrophic nanoflagellates and ciliated protozoa; Landry et al. 2019). Although ubiquitously oligotrophic throughout in its mean state (Table 1; Fig. 1a,b), the north-south geographic shift we observed in growth responses to nutrient supply matches the distribution of predatory tuna (Yellowfin and Bigeye), which are consistently elevated in the NEC region and even further to the south (Fig. 3d;Suzuki 1994;Longhurst 2010;Fonteneau and Hallier 2015). This tuna distribution has been recognized for decades (Suzuki 1994;Fonteneau and Hallier 2015), but has remained a puzzle because nitrate and Chl a in surface waters are ubiquitously low (Lehodey et al. 1997;Longhurst 2010). Our finding of the matching gradient in diatom growth capacity could potentially help to reconcile the counterintuitive observation of elevated concentrations of tuna in the NEC region (Lehodey et al. 1997;Longhurst 2010). A requirement for this link is that transient nutrient (i.e., N + Fe) supply events do occur in such systems, which remains to be investigated (although see Johnson et al. 2010). We have further shown that this gradient is intimately tied to the wind-driven NEC. Accordingly, if our hypothesis is correct, projected future changes in the distribution of wind stress, which controls NEC current flow and Ekman transports (Duan et al. 2017), could have a major impact on fish distribution even if changes are not registered in mean nutrient and Chl a fields.