Photoacclimation by phytoplankton determines the distribution of global subsurface chlorophyll maxima in the ocean

Subsurface chlorophyll maxima are widely observed in the ocean, and they often occur at greater depths than maximum phytoplankton biomass. However, a consistent mechanistic explanation for their distribution in the global ocean remains lacking. One possible mechanism is photoacclimation, whereby phytoplankton adjust their cellular chlorophyll content in response to environmental conditions. Here, we incorporate optimality-based photoacclimation theory based on resource allocation trade-off between nutrient uptake and light harvesting capacity into a 3D biogeochemical ocean circulation model to determine the influence of resource allocation strategy on phytoplankton chlorophyll to carbon ratio distributions. We find that photoacclimation is a common driving mechanism that consistently explains observed global scale patterns in the depth and intensity of subsurface chlorophyll maxima across ocean regions. This mechanistic link between cellular-scale physiological responses and the global scale chlorophyll distribution can inform interpretation of ocean observations and projections of phytoplankton responses to climate change. Trade-offs in the cellular allocation of resources in response to environmental conditions consistently explain the depths of subsurface chlorophyll maxima across the global ocean, according to simulations with a biogeochemical ocean circulation model.

P hotoacclimation is a dynamic physiological response to light availability, manifested as variations of the intracellular concentrations of light-harvesting pigments, which are most commonly observed as chlorophyll. Under low light conditions, phytoplankton increase their chlorophyll to carbon biomass ratio in a cell (cellular chl:phyC) in order to harvest more photons [1][2][3] . The cellular chl:phyC ratio is also constrained by nutrient availability and temperature [1][2][3][4] . Data 5 compiled from incubation experiments reveal that cellular chl:phyC ratios vary by a factor of 60 6 . In field data, cellular chl:phyC varies by a factor of 10 7,8 . Geider 3 proposed an empirical formula using differential equations for the dynamics of cellular chl:phyC, based on the results of incubation experiments. Other empirically based formulae for cellular chl:phyC have been developed to estimate net oceanic primary production from satellite observations 6,9 . Photoacclimation, via changes in cellular chl:phyC, substantially impacts bulk chlorophyll concentrations and the resultant global distribution of oceanic chlorophyll. Vertical chlorophyll profiles generally exhibit maxima below the ocean surface, which are termed subsurface chlorophyll maxima (SCM) or deep chlorophyll maxima. SCM are observed nearly ubiquitously, in subtropical 10 , equatorial [11][12][13] , subpolar 14 , polar [15][16][17][18] , and upwelling regions 19 . Observations reveal substantial regional differences in typical SCM depth, which varies from 100-160 m in the subtropics 10 , to around 40 m in subpolar regions 14 . SCM depth is closely associated with the depth of the nutricline 20 , at which nutrient concentrations sharply increase downward. The mechanisms underlying SCM formation have long been discussed and debated 20 . In subtropical regions, photoacclimation is recognized as an important mechanism underlying SCM formation [21][22][23][24] , where it is attributed mostly to increases in cellular chl:phyC with depth. For equatorial SCM, some studies have highlighted the importance of photoacclimation 25 , while others regard SCM as peaks of biomass 26 . In the subpolar, polar, and upwelling regions, the role of photoacclimation in SCM formation remains unclear.
The purpose of the study is to evaluate to what degree photoacclimation, based on strategic resource optimization, determines global distributions of chlorophyll and hence SCM. We implemented the 0D Flexible Phytoplankton Functional Type (FlexPFT) model of Smith et al. 27 within a global biogeochemical ocean circulation model. The FlexPFT model combines theories for phytoplankton physiology as proposed by Pahlow 28 and others 29,30 . This theory accounts for intracellular resource allocation among structural material, nutrient uptake, and light harvesting associated with the chloroplast ( Supplementary  Fig. 1a). Phytoplankton are assumed to optimize their resource allocation to nutrient uptake versus light-harvesting biomolecules, depending on light, nutrient, and temperature conditions, in order to maximize their growth rate. Under low light and nutrient-replete conditions, phytoplankton increase their resource allocation for light harvesting (chloroplast) while decreasing that for nutrient uptake ( Supplementary Fig. 1b). This resource allocation theory has recently provided the first theoretical derivation 31 of the widely applied empirical Droop cell quota model 32 . Pahlow's model 28 has also explained the results of incubation experiments 33 .

Results and discussion
Role of photoacclimation in SCM formation. The simulation reproduces the regional differences in SCM depth among subpolar, subtropical, and equatorial regions as observed in the North Pacific and Atlantic, with a slightly shallower SCM in the Atlantic subtropical regions (Fig. 1a-d). Simulated SCM is frequently located around the nutricline, where vertical nutrient gradients are steepest (Fig. 1e, f). Since chlorophyll concentration is given by the product of phytoplankton carbon biomass concentration (phyC) (Fig. 1g, h) and cellular chl:phyC (Fig. 1i, j), SCM depth varies with latitude as a consequence of regional differences in the vertical profiles of phyC and cellular chl:phyC. phyC is maximal at the surface in most oceanic regions, (Fig. 1g, h), consistent with observations 34 . Maximal values of cellular chl: phyC (Fig. 1i, j), on the other hand, tend to occur at greater depths (60-130 m). This result highlights the need to account for photoacclimation when modeling and interpreting observations of chlorophyll, which despite being the most commonly observed metric of phytoplankton, is not a good proxy for their biomass 4,6 .
By comparing the vertical profiles of chlorophyll concentration, cellular chl:phyC, and phyC in the subtropical, equatorial, subpolar, Antarctic, and upwelling regions (Fig. 2), respectively, we show that photoacclimation mainly determines the SCM depth. In our model results, SCM do not correspond to the maxima of phyC, which is consistent with previous laboratory results 4 . In subtropical regions, SCM appear at 118 m depth, at which phyC is one-third of its surface value, while cellular chl:phyC is 50 times its value at the surface (Fig. 2f). Since the cellular chl:phyC varies more over depth than does phyC, the SCM depth resides near the depth of the maximum cellular chl:phyC (Fig. 2a), i.e., deep SCM are formed. In subpolar regions, where the SCM are shallow (28 m depth), cellular chl:phyC at the surface is about 1/5 of its maximum at 82 m depth, while phyC at the surface is about 16 times its value at 82 m (Fig. 2h). Therefore, the SCM depth occurs near the depth of the maximum phyC, relatively close to the surface (Fig. 2c). Hence in these regions, the SCM depth is primarily determined by the extent to which cellular chl:phyC at the surface is less than its maximum value, which occurs in the subsurface. This SCM formation mechanism is common to the equatorial (Fig. 2b, g), Antarctic (Fig. 2d, i), and upwelling regions (Fig. 2e, j).
In the resource allocation theory, cellular chl:phyC is the product of the chlorophyll to carbon biomass ratio within the chloroplast (chloroplast chl:phyC) (Fig. 1k, l), and the resource allocation ratio for (i.e., the 'size' of) the chloroplast (Fig. 1m, n). For example, if chloroplast chl:phyC is 0.3 g chl mol C −1 and the resource allocation ratio for the chloroplast is 0.5 (nondimensional), cellular chl:phyC is 0.15 g chl mol C −1 . Chloroplast chl:phyC depends on only light intensity 27 . As light intensity decreases with depth, chloroplast chl:phyC increases from the surface to its maximal depth, which is about 90 m in equatorial regions and 50 m in subpolar regions. Below the maximal depth, chloroplast chl:phyC decreases with depth. The balance between benefits and respiration costs of maintaining chlorophyll 30 determines the depth of maximal chloroplast chl:phyC. The depth of maximal chloroplast chl:phyC is, in turn, a major determinant of the depth of maximal cellular chl:phyC.
Values of the intracellular resource allocation ratio to the chloroplast near 0 (or 1) mean that resources are mainly allocated to nutrient uptake (or light harvesting), with minimal allocation to the competing use. At depths below about 130 m, where light intensity is consistently low, this resource allocation is heavily skewed toward the chloroplast in all oceanic regions. By contrast, above about 130 m in low nutrient regions such as the subtropics, this resource allocation is heavily skewed toward nutrient uptake, and therefore few resources are allocated to the chloroplast. On the contrary, near the surface in high nutrient, e.g., subpolar, regions, phytoplankton resource allocation is more evenly balanced between chloroplast and nutrient uptake. Regional differences in the fractional resource allocation to the chloroplast caused by different nutrient conditions mainly determine regional patterns of cellular chl:phyC. This explains , e, f nitrate (NO 3 ) concentration (μmol N L −1 ), g, h phytoplankton carbon biomass concentration (phyC, μmol C L −1 ), i, j cellular chlorophyll to carbon ratio (chl:phyC in a cell, g chl molC −1 ), k, l chloroplast chlorophyll to carbon ratio (chl:phyC in the chloroplast, g chl molC −1 ), m, n fractional resource allocation to chloroplast. why cellular chl:phyC tends to be very low near the surface in the subtropics, and more generally, why SCM tend to be deeper in low nutrient regions than in high nutrient regions.
Previous observational studies have highlighted the important role of nutricline depth as a determinant of SCM depth. Resource allocation theory provides a more detailed mechanistic basis for the dependence of SCM depth on nutricline depth as well as light. Because nutrient concentrations increase steeply with depth near the nutricline, fractional resource allocation to the chloroplast also increases steeply downward, resulting in sharp increases in cellular chl:phyC. As a result, maximal bulk chlorophyll concentration, i.e., SCM formation, tends to occur around the nutricline. Significant relations among the strong vertical gradients of cellular chl:phyC, nutricline, and SCM obtained in our model are consistent with a previous study for the California Current System 35 .
Global SCM distribution. The simulated global patterns in SCM depth (Fig. 3) are consistent with observed results, as described below. In our simulation, SCM are ubiquitous across the global ocean as in a SCM estimation from observed surface chlorophyll distributions 36 . Global features, such as a remarkable change in SCM depth across the subtropical-subpolar boundary and SCM shallowing in upwelling regions, agree with patterns of the satellite-estimated global SCM distribution 36 . The simulation also successfully describes observed shallowing of SCM depths toward the coasts 20,21 . In summer in the Arctic and Antarctic Ocean, the simulated SCM depth of 50-65 m is in line with observations [15][16][17][18] . In February and August, SCM occur at similar depths in most of oceanic regions. The successful simulation of the global patterns in SCM depth was achieved by explicitly incorporating phytoplankton's resource allocation strategy into the model. Limitations of the model. Our model employs only one generic phytoplankton species, and therefore does not explicitly resolve community composition. In case studies in which the chlorophyll-specific initial slope of growth versus irradiance, a I , is changed ( Supplementary Fig. 2), increasing a I causes deeper SCM and increase in chlorophyll concentration around the SCM, which improves reproducibility of chlorophyll around SCM, but at the expense of making surface chlorophyll concentration deviate far from the observations. In the standard case, our model is tuned in order to simulate both surface and deep chlorophyll concentration to a reasonable degree with the single generic phytoplankton species. However, different photoacclimation responses have been observed for different phytoplankton lineages 37 . In the future, by introducing surface and deep species having different traits for irradiance, our model can be expected to reproduce more realistic SCM depth and chlorophyll concentration around the SCM. Another limitation of our model is that it represents chlorophyll concentration only at scales larger than about 10 m due to its vertical resolution. For finer-scale vertical variation of chlorophyll concentration, aggregation related to swimming behavior or buoyancy control 20 may play some roles. A mesoscale eddy has been observed to alter SCM depth substantially, from 80 to 120 m 38 , which is beyond the scope of this study due to the coarse resolution of our physical model.
Relevance to future climate change. Climate change is expected to alter the oceanic environment, especially temperature, stratification, and nutrient supply rates. Its impact on phytoplankton 39 , which are the sustaining basis of marine food chains, will be crucial for humankind. Our research has proceeded on the premise that improving models' ability to reproduce the observed chlorophyll distribution will yield more reliable future projections. By establishing a mechanistic link between the individual-level photoacclimation response of phytoplankton and chlorophyll distributions at the global scale, our results provide a more sound basis for interpreting chlorophyll observations and predicting how oceanic phytoplankton in particular, and hence marine ecosystems more broadly, will respond to future climate change.  In accordance with Pahlow's resource allocation theory 28 , the FlexPFT model assumes that resources are allocated among structural material, nutrient uptake and, light harvesting ( Supplementary Fig. 1a). The fraction of structural material is assumed to be Q s /Q, where Q is the nitrogen cell quota, which is the intracellular nitrogen to carbon ratio (mol N mol C −1 ), and Q s is the structural cell quota (mol N mol C −1 ) given as a fixed parameter. The fraction of nutrient uptake is defined as f V (non-dimensional), so that the residual fraction available for light-harvesting is equal to ð1 À Q s Q À f v Þ. Optimal uptake kinetics further sub-divides the resources allocated to nutrient uptake between surface uptake sites (affinity) and enzymes for assimilation (maximum uptake rate), the fraction of which is given by f A and  (1 − f A ), respectively. Under nutrient-deficient conditions, the number of surface uptake sites (and hence affinity) increases, while enzyme concentration (hence, maximum uptake rate) decreases. The FlexPFT model assumes instantaneous resource allocation, which means that resource allocation tracks temporal environmental change with no lag time. It has elsewhere been demonstrated that an instantaneous acclimation model provides an accurate approximation of a fully dynamic acclimation model 44 .

Methods
We assume that acclimation responds to daily-averaged environmental conditions, which are used to calculate the optimal values of f V , f A , and Q as Phytoplankton growth rate per unit carbon biomass (day −1 ), μ, is given by whereμ I is the potential carbon fixation rate per unit carbon biomass (day −1 ), ζ N is the energetic respiratory cost of assimilating inorganic nitrogen (0.6 mol C mol N −1 ), andV N is the potential nitrogen uptake rate per unit carbon biomass (mol N mol C −1 day −1 ). Equation (1) represents the balance of net carbon fixation and respiration costs of nitrogen uptake, which are proportional to the fraction of resource allocation.V N ð½ N; TÞ iŝ whereÂ 0 andV 0 are the maximum value of affinity and maximum nitrogen uptake rate. From here, we will explain how the optimized values such as f o V , f o A , and Q o are calculated. The optimal fraction of resource allocation to affinity, f o A , is given by which is derived by substituting Eqs. (18) and (19) in Smith2016 into Eq. (17). Fð TÞ is temperature dependence, defined as where E a is the parameter of the activation energy of 4.8 × 10 4 J mol −1 , R is the gas constant of 8.3145 J (mol K) −1 , and T ref is the reference temperature of 20°C. Optimization for light-harvesting is described below. The potential carbon fixation rate per unit carbon biomass (day −1 ),μ I (day −1 ), in Eq. (1) iŝ whereμ 0 and k Fe are the maximum carbon fixation rate and half saturation constant for iron, respectively. S specifies the dependence of light. Defininĝ Iron limitation is imposed by substitutingμ 0 toμ limFe 0 in all equations in Smith2016. S is defined as where a I is the chlorophyll-specific initial slope of growth versus irradiance.Θ o , optimal chloroplast chl:phyC (g chl (mol C) −1 ), iŝ where constant parameters ζ chl and R chl M are the respiratory cost of photosynthesis (mol C (g chl) −1 ) and the loss rate of chlorophyll (day −1 ), respectively. L d is the fractional day length in 24 h. W 0 is the zero-branch of Lambert's W function. I 0 is the threshold irradiance below which the respiratory costs overweight the benefits of producing chlorophyll: The optimal fraction of resource allocation to nutrient uptake, f o V , is The optimal nitrogen cell quota, Q o is Optimal cellular chl:phyC (g chl (mol C) −1 ), Θ o , is which is the multiplication of the fraction of resource allocation to light-harvesting and optimal chloroplast chl:phyC. The cellular chl:phyC and chloroplast chl:phyC in Figs. 1 and 2 are optimal cellular chl:phyC, Θ o , and optimal chloroplast chl: phyC,Θ o , respectively. The relation in Eq. (12) is displayed in Fig. 1i-n. If we artificially turn off the optimization of resource allocation by applying the constant Q o and f o V to the all grid points, optimal cellular chl:phyC (Fig. 1i,j) only depends on optimal chloroplast chl:phyC (Fig. 1k, l), and therefore significant variation of SCM depth across the equatorial, subtropical, and subpolar regions is not reproduced.
In the above equations, Eqs. (3), (8), (10), (11), and (12), optimized values related to acclimation processes are obtained and then used in calculating the phytoplankton growth rate. Phytoplankton growth rate per unit carbon biomass (day −1 ), μ, in Eq. (1) is calculated at each time step: whereμ I ðI; T; ½Fe d Þ is obtained by substituting I, T, and [Fe d ] for I, T, and ½Fe d in Eq. (5), respectively. Note that the model calculates circadian variation in solar irradiance, I, and therefore the phytoplankton growth rate, μ, reaches its maximum at noon local time and is zero during night. On the other hand, phytoplankton optimization is assumed to respond to daily-averaged conditions. The FlexPFT model introduces phytoplankton respiration proportional to chlorophyll content, which is another important originality of Pahlow's resource allocation theory 30,33 .
The carbon biomass-specific respiratory costs of maintaining chlorophyll, R chl , is The growth rate per unit nitrogen biomass, μ N , is equal to that per unit carbon biomass, μ. Instantaneous acclimation assumes that the quota of nitrogen to carbon biomass obtained by phytoplankton growth is equal to the nitrogen quota in a cell: When temporal Q o change occurs, to satisfy the mass conservation, carbon or nitrogen biomass is adjusted with the other fixed. The FlexPFT fixes carbon biomass, while the FlexPFT-3D fixes nitrogen biomass to the temporal Q o change.
The rate of change in the phytoplankton nitrogen concentration, [p N ], except for the advection and diffusion terms is given by the following equation: where R cnst and M p are the coefficient of respiration not related to chlorophyll concentration and mortality rate coefficient, respectively. The extracellular excretion is where γ ex is the coefficient of extracellular excretion. The grazing term is represented by where G 20deg is the maximum grazing rate at 20°C, and [z N ] is zooplankton concentration. Temperature dependency, F(T), is obtained by substituting T for T in Eq. (4). a H is the parameter controlling Holling-type grazing, which takes a value from 1 to 2. k H is the grazing coefficient in Holling-type grazing.
Once [p N ] is calculated, phytoplankton carbon concentration (mol C L −1 ), and chlorophyll concentration (g chl L −1 ) are uniquely determined in an environmental condition, without prognostic calculation. Therefore, an instantaneous acclimation model can represent stoichiometric flexibility with lower computational costs compared with a dynamic acclimation model 44 .
Model validation. The spatial pattern of simulated annually mean chlorophyll at the ocean surface agrees with that of satellite observation 45 (Supplementary Fig. 3). The model reproduced the contrast of the surface chlorophyll concentration between subtropical and subpolar regions, although simulated surface chlorophyll concentration in subtropical regions is lower than that of the observation partly due to the lack of nitrogen fixers. Nitrogen fixation is estimated to support about 30-50% of carbon export in subtropical regions 46,47 . Simulated surface chlorophyll distribution in the Pacific equatorial region is close to the observed.
Our model properly simulates the meridional distribution of nitrate compared with that of observations 48 (Supplementary Fig. 4). The simulated horizontal distribution of primary production is consistent with that estimated by satellite data 9,49 (Supplementary Fig. 5), although simulated primary production is underestimated in subtropical regions, associated with the underestimation of surface chlorophyll in these regions ( Supplementary Fig. 3).

Code availability
MRI.COM3 is available on reasonable request from Meteorological Research Institute. The website (https://mri-ocean.github.io/) describing the details of use is currently only available in Japanese, so if you have any questions, please contact the corresponding author. The core code of the FlexPFT-3D is available on reasonable request from the corresponding author. The NCAR Command Language (NCL) code used to create the figures is available at https://www.ncl.ucar.edu/.