Vertically migrating phytoplankton fuel high oceanic primary production

Marine net primary production (NPP) is remarkably high given the typical vertical separation of 50–150 m between the depth zones of light and nutrient sufficiency, respectively. Here we present evidence that many autotrophs bridge this gap through downward and upward migration, thereby facilitating biological nutrient pumping and high rates of oceanic NPP. Our model suggests that phytoplankton vertical migration (PVM) fuels up to 40% (>28 tg yr−1 N) of new production and directly contributes 25% of total oceanic NPP (herein estimated at 56 PgC yr−1). Confidence in these estimates is supported by good reproduction of seasonal, vertical and geographic variations in NPP. In contrast to common predictions, a sensitivity study of the PVM model indicates higher NPP under global warming when enhanced stratification reduces physical nutrient transport into the surface ocean. Our findings suggest that PVM is a key mechanism driving marine biogeochemistry and therefore requires consideration in global carbon budgets. Phytoplankton vertical migration has a role in nutrient pumping and primary productivity in the oceans. Here the authors quantify the total amount of oceanic net primary productivity facilitated by this bio-pumping, under present and future warming conditions.

A nnual rates of net primary production (NPP) in the ocean typically exceed 50 gC m −2 yr −1 (refs. 1,2 ). Prevalence of high NPP rates, however, disagrees with the widespread observation in low and middle latitudes of nutrient depletion in the upper water column where light intensity suffices for high photosynthesis rates. Saturating light conditions rarely extend deeper than 30-60 m, which in vast areas of the oligotrophic ocean is roughly 100 m above the nutrient-rich chemocline. This separation of the sunlit zone from the chemocline challenges unicellular autotrophs, which are believed to be short-lived and therefore to require light and nutrients simultaneously. Furthermore, high NPP rates near the surface sustain a continuous export flux out of the euphotic zone, which must be counterbalanced by nutrient inputs. However, diffusive fluxes from the nutrient-rich deep layers are often limited by a stratified pycnocline, which typically lies above the chemocline 3 . Lateral transport has been suggested by modelling studies as a nutrient source for the surface layer of oligotrophic gyres, but inorganic nutrients such as nitrate were simulated to enter only the margins of the gyres 4 . Finally, atmospheric nutrient deposition and N fixation by cyanobacteria provide roughly 10% to and maximum of 50% of the export flux and hence are insufficient to sustain productivity [5][6][7][8] , especially since phosphorus deposition rates are negligibly low 9 .
Therefore, conventional biogeochemical (BGC) models partially require weakly supported assumptions to reproduce high regional or global NPP rates (Supplementary Section 4), while detailed examinations of model results 10,11 and observations 5,9 suggest a substantial unknown source of nutrients to the surface layers.
One such alternative transport mechanism could be biological pumping by phytoplankton, carrying assimilated nutrients and migrating from the chemocline towards the surface. Vertical migration in the ocean has long been known for a few, fast-moving taxa such as mat-forming diatoms, dinoflagellates and the cyanobacterium Trichodesmium 12-15 but has so far not been included in global BGC modelling because of the limited abundance of those taxa. A much greater relevance of biological pumping has recently been proposed to follow from slow but long-term continuous migration over the scale of days and weeks by the large fraction of motile phytoplankton within the autotrophic community 16 . A Lagrangian model for phytoplankton vertical migration (PVM) applied at several marine stations indicated that substantial amounts of nutrients can be pumped up by migrating bulk phytoplankton.
In this Article, we test whether PVM can solve the enigma of high NPP rates in the nutrient-depleted surface ocean. We assess its relevance for global nutrient and carbon cycles and, more specifically, quantify the total amount of oceanic NPP facilitated by vertical migration and bio-pumping, both for present-day climatic conditions and in a warming future ocean.
Our model differentiates between two fractions of the phytoplankton community: vertical migrators (typically eukaryotes) and passive drifters (many prokaryotes) 16 . Immobile cells grow in the euphotic zone according to local light conditions while assuming moderate limitation by remineralized, deposited or fixed nitrogen and co-limitation by other nutrients as approximated by residual surface nitrate (Methods). Mobile cells are described as active Lagrangian particles, which are individual objects moving along a trajectory. In our model, cells commute between the chemocline, where they increase their nutrient/C ratios, and the upper sunlit layers, where photosynthetic carbon assimilation and utilization of stored nutrients reduce those ratios (Fig. 1). The biological upward N flux is quantified as the bulk concentration of intracellular nitrogen within mobile autotrophs above the chemocline times their mortality rate. The Lagrangian model is driven by time-dependent environmental conditions such as temperature and chemocline depth, as well as measured concentrations of chlorophyll (CHL) and nitrate at the surface (Methods, Supplementary Fig. 1 and ref. 16 ). In the present study, these input data are processed in two different set-ups. First, applications to the monitoring stations Bermuda Atlantic Time Series (BATS), Hawaii Ocean Time Series (HOT), Northwest Pacific subarctic (K2) and Northwest Pacific subtropical (S1) test the validity of our approach with respect to NPP profiles under variable but directly monitored environmental conditions. Second, for a global application, we compiled the present-day climatological seasonality of environmental data at a spatial resolution of 0.5°. In this global set-up, the impact of vertical migration is also evaluated in a sensitivity study using changes in chemocline depth Vertically migrating phytoplankton fuel high oceanic primary production Kai Wirtz 1 ✉ , S. Lan Smith 2 , Moritz Mathis 1 and Jan Taucher 3 Marine net primary production (NPP) is remarkably high given the typical vertical separation of 50-150 m between the depth zones of light and nutrient sufficiency, respectively. Here we present evidence that many autotrophs bridge this gap through downward and upward migration, thereby facilitating biological nutrient pumping and high rates of oceanic NPP. Our model suggests that phytoplankton vertical migration (PVM) fuels up to 40% (>28 tg yr −1 N) of new production and directly contributes 25% of total oceanic NPP (herein estimated at 56 PgC yr −1 ). Confidence in these estimates is supported by good reproduction of seasonal, vertical and geographic variations in NPP. In contrast to common predictions, a sensitivity study of the PVM model indicates higher NPP under global warming when enhanced stratification reduces physical nutrient transport into the surface ocean. Our findings suggest that PVM is a key mechanism driving marine biogeochemistry and therefore requires consideration in global carbon budgets.
as projected by a global Earth system model under a global warming scenario ( Supplementary Fig. 1).

effect on global NPP and nutrient inputs
Our account of PVM reproduces well the observed NPP profiles at the four marine stations (agreement between climatological profiles of model results and observations in Extended Data Fig. 1), especially at HOT and BATS (Fig. 2), and also depth-integrated NPP (NPP Σ ) calculated in the global set-up (Extended Data Fig. 2), as described in Supplementary Section 1. The global simulation also generates reasonable seasonal variations in NPP Σ in mid and high latitudes (Extended Data Fig. 3), with PVM contributing most to annual NPP at low latitudes (Fig. 3a). Despite a systematic underestimation in coastal seas and upwelling regions (Supplementary Section 1), the annual global rate of 56 PgC yr −1 calculated here exceeds current empirical estimates such as 49 PgC yr −1 by ref. 2 or 52 PgC yr −1 by ref. 17 .
Simulated NPP decreases substantially in a control model variant neglecting PVM (Fig. 3b), especially in the tropical oceans, except for the Arabian Sea and western Central Atlantic. Without vertical migration, our model predicts a global oceanic NPP of 41 PgC yr −1 , 25% lower than the prediction with migration (Extended Data Fig. 4), which is in the upper range of the 7-28% contributions by migrating primary producers previously calculated for the reference stations 16 .
However, for about one-third of the PVM areas, migration speeds below 0.7 m d −1 and migration distances below 12 m were obtained (Extended Data Figs. 5-7). At such low speeds and distances, active transport can hardly be distinguished from physical transport by diffusion. In these cases, the subsurface chlorophyll maximum (SCM) may indeed be dominated by non-migrating phytoplankton, especially when the CHL concentration at the SCM exceeds the one at the surface by a factor of only 1-3. This means that our estimated relative contributions by migrating primary producers should be taken as an upper threshold.
In addition to direct enhancements of NPP, substantial nutrient pumping by migrating phytoplankton fuels regenerated production by passively drifting phytoplankton. Biological pumping contributes an upward nutrient flux of about 5-10 mmol m −2 d −1 N in vast areas of the global ocean (Extended Data Fig. 8), which exceeds by more than an order of magnitude the estimates of combined atmospheric inputs (N deposition and N 2 fixation) 18 or of the particulate export flux 19 . Nevertheless, biological nutrient pumping generally extends tens of metres above the chemocline (see migration distance 2δz in Extended Data Fig. 5) and hence only partially reaches the surface zone where most NPP Σ occurs. Therefore, we z SCM δz δz

Fig. 1 | Biological nutrient pumping by vertically migrating phytoplankton cells.
Cells with low intracellular nutrient content (indicated by small ellipses) settle or swim down to the chemocline while nutrient-rich cells ascend to the sunlit surface layers. Their vertical habitat spans from z SCM - δz to z SCM + δz with the average position z SCM marking the depth of the subsurface chlorophyll maximum and vertical distance δz the migration amplitude. The up and down migration leads to a separation of cells having low versus high nutrient/C ratios at the same depth, depending only on their migration history. This separation is difficult to observe in situ and impossible to resolve in the Eulerian grid representation of conventional BGC models.

uplifting the SCM
Indirect evidence for the existence of bulk PVM originates from the observed accumulation of CHL at depths substantially shallower than the nutrient-rich chemocline depth zN (Fig. 4). Although the photo-acclimative variation of CHL/C ratios, which our model includes, may determine SCM formation to a large extent 20 , as also reproduced by photo-acclimation models 21,22 , this systematic uplift of SCM depths above the chemocline is difficult to explain without migration. A version of our model neglecting migration predicts SCM near the chemocline, with z SCM ≈ 0.9 (z N ) ( Supplementary  Fig. 3). In the reference run of the model, migration confers greater advantage as the chemocline deepens, so that displacements strongly increase for z N > 100 m. This pattern quantitatively agrees with data merged from the four open-ocean time-series stations (Fig. 4). SCM centre positions (z SCM ) 20-50 m above the chemocline are also frequently observed at other sites, such as in the Atlantic gyres 23 , or are found to be representative for z SCM > 100 m (ref. 24 ).
Simulations including migration generate a nonlinear relation between z SCM and z N , which is identical to the observed pattern (statistically fitted relationships are not distinguishable, Fig. 4). A very similar relationship, including a threshold around z N = 100 m, appears in an independent dataset covering the global ocean 2 , supporting the generality of the SCM displacement. This uplifting extends the classical view of the SCM as located near the top of the chemocline, which was based mainly on data for z SCM < 100 m (ref. 20 ). SCM formation in situations with z N > 100 m cannot be seamlessly reproduced by classical models using realistic parameterizations for light attenuation by water and particles and light affinity of autotrophs. Indeed, BGC models have not yet satisfactorily reproduced profiles at BATS and HOT despite the prominence of those extensive time-series observations. An exception was the early study of Doney et al. 10 , who used a climatology of BATS data before 1993. In these years, chemocline was at times unusually shallow (Fig. 2), so that averaged z N was around only 80-90 m. Similarly, station K2 in general has relatively low z N , for which our model predicts a small contribution of active migration. The model study by Li et al. 25 , which did not account for vertical migration, showed overall good skill in resolving subsurface CHL concentrations and SCM positioning at K2 (Supplementary Section 1). Like many models such as the one of Li et al. 25 , our model resolves flexible CHL/C ratios (for example, with a typical increase by a factor 2.5 from the surface to 150 m depth; Supplementary Fig. 5 16 ), which facilitate the development of an SCM and overall predict a decreasing biomass of non-migratory autotrophs with depth and a relatively small variation in the concentration of total phytoplankton C over depth down to the SCM, in agreement with a recent synthesis of particle backscatter data 24 .

Fundamental traits shaping global biogeography and diversity
Active locomotion or buoyancy regulation is found in more than two-thirds of phytoplankton species (Table 1). Another ubiquitous trait is high plasticity in nutrient/carbon ratios. These and further evidences support the widespread existence and relevance of PVM (Supplementary Section 2). Our simulation suggests that most PVM is realized at migration speeds (v) below 1 m d −1 , thus well within the capability range of motile/buoyant species (Extended Data Fig. 6 and Table 1). Yet rapid migration at v above 10 m d −1 (Extended Data Fig. 7) is simulated in areas with low NPP Σ (Fig.  3 and Extended Data Fig. 6), in proximity to the coast or with z SCM > 120 m (Extended Data Fig. 9). In particular, the combination of high migration speed and deep SCM can lead to migration distances of up to 80-100 m (Extended Data Fig. 5). Diel, long-range migration of phytoplankton in coastal proximity, is in agreement with observed migration of dinoflagellates in Oslo fjord 13

unexpected response to climate change
The distribution of functional traits well determines how the phytoplankton community reacts to environmental change 31,34 . The  biogeography of behavioural traits together with the relevance of phytoplankton behaviour compared with physical transport mechanisms may therefore be important for understanding the response of marine ecosystems to external drivers. This motivated us to conduct a sensitivity study of our global application related to a global change scenario. Oceanic production has been suggested to decrease under anticipated climate change. Warming of surface waters, increasing stratification and, concomitantly, reduced nutrient supply from deeper layers may cause NPP to decrease 35 . In a projection by the Max-Planck-Institute (MPI) Earth system model (ESM; low-resolution (LR) Coupled Model Intercomparison Project 5 (CMIP5) version) within the emission scenario representative concentration pathway (RCP) 8.5 36 , mixed-layer depth (MLD) overall decreases-with few exceptions (Supplementary Fig. 8), and within a century (from the time slice 1971-2000 to 2071-2100), the corresponding simulated oceanic production decreases by 8.5 PgC yr −1 (14%; Supplementary Fig. 5). This projection of enhanced stratification and reduced production lies well within the CMIP5 multimodel ensemble spread 35,37,38 . Accordingly, our model calculates lower NPP Σ by non-migrating phytoplankton as a result of a lower thickness of the productive water layer, but also increased subsurface productivity due to uplifted SCM and migratory phytoplankton ( Supplementary Figs. 1 and 3). This pattern is substantiated using a simple but also more mechanistic one-dimensional (1D) set-up where we iteratively run the Lagrangian PVM model with a Eulerian coupled physical-biological model, which resolves major nutrient sources and sinks such as uptake by all phytoplankton groups, N 2 fixation, regeneration and turbulent diffusive mixing (Supplementary Section 5). In this 1D set-up, altered temperature and wind forcing induces steeper vertical temperature gradient and a reduction of both biomass and productivity of non-mobile phytoplankton and of the chemocline depth due to decreased mixing; however, a reduced chemocline depth also leads to an increase in the concentration of mobile phytoplankton, especially at intermediate depths (Supplementary Fig. 6). In the global 3D application, substantial chemocline shallowing predicted for the Antarctic Circumpolar Current, North Pacific Current and North Atlantic areas south of the Greenland Current (Supplementary Figs. 7 and 8) induces only minor to moderate changes in NPP Σ because relatively high CHL concentrations at the surface (Supplementary Fig.  9) disadvantage vertical migration as a viable strategy. By contrast, the changed distribution of chemocline depths translates to much enhanced growth of vertically migrating species in regions that combine low surface CHL (in our approach, equivalent to a small standing stock of drifters) with an intermediate SCM depth of 40-60 m. For example, in the North Equatorial Current, the western South Equatorial Current, the eastern Indian Ocean, the East Greenland Current and the region between the Agulhas and the Antarctic Circumpolar Current, even moderate decreases in chemocline depth (Supplementary Fig. 8) widen the habitat of vertically migrating species towards the sunlit surface layers (Extended Data Fig. 10) and, consequently, enhance their production rates. Overall, increased NPP Σ by migrating phytoplankton more than compensates the decreased NPP Σ by non-migrating phytoplankton (Fig. 5), leading to an increase in global oceanic NPP from 56.1 (in the reference run) to 58.4 PgC yr −1 . Some of the hotspots of NPP increase are proximate to regions where conventional models, including the MPI-ESM-LR, predict the greatest NPP decline 35 , such as the Equatorial Countercurrent ( Supplementary Fig. 5). The opposing responses of ESMs and our model to the same physical forcing challenges currently widespread approaches to predicting climate change effects.

a common strategy for sequential multi-resource acquisition
Phytoplankton is in most observational studies and in almost all coupled BGC models considered a passive tracer, lacking behavioural response. Consequences of neglecting phytoplankton migration in both satellite-based and mechanistic models are briefly discussed in Supplementary Section 4. In lakes, vertical migration by phytoplankton at high speed or over short distances was suggested earlier as an efficient strategy to enhance production, as summarized in Supplementary Section 3. Similarly, vertical migration observed for many crustaceans and other zooplankton taxa has been interpreted as an adaptive strategy in terms of both predator avoidance and resource exploitation 39 . A wide variety of organisms within the pelagic zone, such as cetaceans, employ migration to acquire multiple resources at distinct vertical positions. Even in ocean sediments, bacteria move up and down to obtain nitrate from above and sulfide from below 40 . Our study suggests that vertical migration by bulk marine phytoplankton is an equally widespread form of such migratory acquisition strategies that has not received sufficient attention so far, particularly considering its global-scale implications.

Conclusions
Our study shows that vertical migration of bulk phytoplankton can substantially fuel high NPP in the nutrient-depleted surface ocean. While it remains technically challenging to confirm this mechanism and its implications in situ, our model results may stimulate the development of new observational approaches for directly testing the presence of PVM. The results also suggest that the picture of phytoplankton as continually active drifters can explain observed vertical profiles as well as depth-integrated rates of NPP over a wide range of ocean regions much better than conventional models that represent phytoplankton as passive drifters. The vertical structure and functioning of pelagic ecosystems may therefore be controlled by phytoplankton behaviour to a much greater degree than previously thought and not solely by physical mechanisms such as nutrient supply via turbulent mixing. The substantial impact of this strategy found in our global NPP calculations points to a systematic underestimation of oceanic production, which probably exceeds our conservative estimate of 56 PgC yr −1 . Notably, biological nutrient pumping implies a counteracting NPP response to reduced physical mixing, which is relevant for assessing the response of the global ocean ecosystem to a warming future climate. Our results hence emphasize the need to revise classical grid-based model approaches and observational techniques by changing the perception of phytoplankton from passive to active drifters.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-022-01430-5.

Local and global model set-ups.
We applied the Lagrangian model for PVM in two settings: a 1D set-up for each of the four time-series observation sites (BATS, HOT, K2 and S1) and a 3D set-up covering the global ocean. In both set-ups, the previously described PVM model 16 used the documented parameterization, except for an additional co-limitation term (see the following). Their major dependencies on boundary forcing, especially the response of PVM to variations in chemocline depth, are illustrated in Supplementary Figs. 1 and Fig. 3.
Forcing and validation of the model in the 1D set-up relied on publicly available data for vertical profiles of chlorophyll, nutrients and NPP as described in the predecessor study 16 , which especially includes the NPP data from HOT 54 and BATS 55 . We also integrated downward short-wave radiation fluxes from the JRA55 re-analysis 56 , aggregated to daily (24 h) values from the originally provided 3-hourly data. In addition, estimates for MLD available for HOT were used to better reconstruct seasonal variations in mixing lengths, which otherwise were calculated using sea surface temperature data (Supplementary Section 5 16 ).
Forcing data for the 3D set-up originated from global datasets on a 0.5° grid. Nitrate profiles, covering the global ocean based on World Ocean Circulation Experiment and World Ocean Database 2009 57 , were available through the CSIRO Atlas of Regional Seas (CARS2009) database 58 . Surface values NO 3,0 were directly extracted ( Supplementary Fig. 4), whereas the profiles were further processed to z N using the algorithm documented before 16 (see annual average in Supplementary Fig. 7). Surface chlorophyll ( Supplementary Fig. 9), sea surface temperature and cloud cover were obtained from the European Space Agency Climate Change Initiative [59][60][61] . The datasets for the period 1999-2010 integrate observations from three satellites. All global boundary data were pooled into four seasonal time slices: December-February (DJF), March-May (MAM), June-August (JJA) and September-November (SON). Net short-wave radiation was calculated using an astronomical formula based on date, geographic position and cloud cover 62 .
The global set-up was run for the four seasonal climatologies in three simulation experiments: (1) a reference run, (2) with PVM switched off and (3) with an altered chemocline field emulating a projected reduction in MLD in a future climate. These future changes in z N were estimated using an RCP 8.5 climate change scenario run of the MPI-ESM-LR 36 . After re-gridding, absolute changes in MLD between the recent-past  and future (2071-2100) distributions were added to the existing distribution of chemocline depth. We chose MLD instead of z N , which were both calculated by the MPI-ESM, because changes in MLD were responsible for changes in z N , regardless of which specific BGC processes mediated those changes in the MPI-ESM. To avoid unrealistic adjustment, these differences in MLD were cut so that chemocline depths were altered by at most 50% ( Supplementary Fig. 8). The consistency of this treatment was qualitatively confirmed by the alternative 1D set-up, which explicitly accounts for most processes shaping the nutrient profile, hence z N , and for feedbacks between changed z N , migration behaviour and ecosystem dynamics ( Supplementary Fig. 6).
Optimization of migration traits. Active migration of unicellular autotrophs reflects tactic responses to external cues such as light as well as intracellular biochemical status 50,51,63 , which some 1D rule-based models already include for fast-moving phytoplankton [63][64][65][66] . Our model assumes near optimality of the migration pathway or strategy, respectively 16 . A variable migration strategy is defined by average positioning (z SCM ), migration amplitude (vertical distance 2δz, Fig. 1) and migration speed (v). A search algorithm identifies the optimal and near-optimal strategies in terms of net growth rate realized at a given environment (incident surface light, water temperature, surface CHL concentration, chemocline depth, surface nitrate concentration). The local optimization implies that the 3D global application neglects any lateral exchange. Local, close near-optimal strategies may follow from different migration pathways so that the averaged migration traits can already reflect an aggregation over diverse speeds or central positions.
Parameterization of co-limitation. The model has been modified from the original version 16 by adding a co-limitation factor f colim for immobile producers depending on surface nitrate values NO 3,0 (=NO 3 (z = 0), as part of the CARS2009 data 58 ) (see the preceding; Supplementary Fig. 4). High NO 3,0 can indicate limitation by iron or other nutrients, which will lower phytoplankton growth and production rates proportional to f colim . The factor is parameterized using the specific co-limitation variation f 0 colim : f colim = 1 − f 0 colim + f 0 colim exp(−NO 2 3,0 ).
Substantial residual surface nitrate has been observed at K2, where recurrent iron limitation has been reported. While the original model ( f 0 colim = 0) overestimates NPP Σ at K2 (not shown), with our setting f 0 colim = 0.5 it underestimates the observations. We decided to keep this parameterization because we aim at a conservative lower estimate of global NPP.

Data availability
All datasets used as model forcing are publicly available (see Methods).

Code availability
The MATLAB code required to produce all model results is available at https://doi. org/10.5281/zenodo.6608970. references acknowledgements ESA, CSIRO and JMA are acknowledged for making their data available. The work was supported by the Helmholtz society via the programme 'Changing Earth' . K.W. was supported by the BMBF project Multiple Stressors on North Sea Life-MuSSeL with project no. 03F0862A; M.M. was supported by Germany's Excellence Strategy EXC2037 (CLICCS-Climate, Climatic Change and Society) with project no 390683824. We thank P. Anugerahanti and C. Völcker for sharing their model set-ups for the HOT station. We also acknowledge the contributions of the many scientists and technicians who produced and made available the data from the time-series observations at stations BATS, HOT, K2 and S1.

Funding
Open access funding provided by Helmholtz-Zentrum hereon GmbH (4216).
Correspondence and requests for materials should be addressed to Kai Wirtz.
Peer review information Nature Climate Change thanks Qian Li, Tracy Villareal and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.
Reprints and permissions information is available at www.nature.com/reprints.