A meta-analysis on environmental drivers of marine phytoplankton C : N : P

The elemental stoichiometry of marine phytoplankton plays a critical role in global biogeochemical cycles through its impact on nutrient cycling, secondary production, and carbon export. Although extensive laboratory experiments have been carried out over the years to assess the influence of different environmental drivers on the elemental composition of phytoplankton, a comprehensive quantitative assessment of the processes is still lacking. Here, we synthesized the responses of P:C and N:C ratios of marine phytoplankton to five major drivers (inorganic phosphorus, inorganic nitrogen, inorganic iron, irradiance, and temperature) by a meta-analysis of experimental data across 366 experiments from 104 journal articles. Our results show that the response of these ratios to changes in macronutrients is consistent across all the studies, where the increase in nutrient availability is positively related to changes in P:C and N:C ratios. We found that eukaryotic phytoplankton are more sensitive to the changes in macronutrients compared to prokaryotes, possibly due to their larger cell size and their abilities to regulate their gene expression patterns quickly. The effect of irradiance was significant and constant across all studies, where an increase in irradiance decreased both P:C and N:C. The P:C ratio decreased significantly with warming, but the response to temperature changes was mixed depending on the culture growth mode and the growth phase at the time of harvest. Along with other oceanographic conditions of the subtropical gyres (e.g., low macronutrient availability), the elevated temperature may explain why P:C is consistently low in subtropical oceans. Iron addition did not systematically change either P:C or N:C. Overall, our findings highlight the high stoichiometric plasticity of eukaryotes and the importance of macronutrients in determining P:C and N:C ratios, which both provide us insights on how to understand and model plankton diversity and productivity.

1 Introduction

Elemental stoichiometry of biological production in the surface ocean plays a crucial role in the cycling of elements in the global ocean. The elemental ratio between carbon, nitrogen (N), and phosphorus (P) in exported organic matter expressed in terms of the $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio helps determine how much atmospheric carbon is sequestered in the deep ocean with respect to the availability of limiting nutrients. On geologic timescales, the N:P ratio reflects the relative availability of nitrate with respect to phosphate, both of which are externally supplied from the atmosphere via nitrogen fixation and/or continents via river supply and lost by denitrification and burial (Broecker, 1982; Lenton and Watson, 2000; Redfield, 1958; Tyrrell, 1999). On shorter timescales, the average stoichiometry of exported bulk particulate organic matter reflects the elemental stoichiometry of phytoplankton (Bonachela et al., 2016; Garcia et al., 2018; Martiny et al., 2013b), with additional influences from biological diversity and secondary processing of organic matter by zooplankton and heterotrophic bacteria. In the face of global change, understanding and quantifying the mechanisms that lead to variability in $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios are crucial in order to have an accurate projection of future climate change.

A key unresolved question is what determines $\mathrm{C}:\mathrm{N}:\mathrm{P}$ of individual phytoplankton. Phytoplankton grows in the upper light-lit layer of the ocean, where the amount of inorganic nutrients, light, and temperature vary spatially and temporally. Laboratory studies show that these fluctuations trigger responses at the cellular level, whereby cells modify resource allocation in order to adapt optimally to their ambient environment (Geider and La Roche, 2002). For example, phytoplankton may alter resource allocation between the P-rich biosynthetic apparatus, N-rich light-harvesting apparatus, and C-rich energy storage reserves (Moreno and Martiny, 2018). Under a typical future warming scenario, the global ocean is expected to undergo changes in nutrient availability, temperature, and irradiance (Boyd et al., 2010). These changes are likely to have profound effects on the physiology of phytoplankton (Finkel et al., 2010; van de Waal et al., 2010), and observations show that competitive phytoplankton species can acclimate and adapt to changes in temperature, irradiance, and nutrients on decadal timescales (Irwin et al., 2015). Numerous laboratory and field experiments have been conducted thus far to study the relationship between the $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio of phytoplankton and environmental drivers. It is, however, challenging to synthesize those studies and generalize the response of phytoplankton $\mathrm{C}:\mathrm{N}:\mathrm{P}$ to changes in environmental drivers. Individual studies employ different sets of statistical analyses to characterize the effects of the environmental driver(s) on elemental ratios, ranging from a simple t test to more complex mixed models, which makes interstudy comparisons challenging. In addition, since environmentally induced trait changes are driven by a combination of plasticity (acclimation), adaptation, and life history (Collins et al., 2020; Ward et al., 2019), stoichiometric responses of phytoplankton can be variable even amongst closely related species.

Meta-analysis/systematic review is a powerful statistical framework for synthesizing and integrating research results obtained from independent studies and for uncovering general trends (Gurevitch et al., 2018). The seminal synthesis by Geider and La Roche (2002), as well as the more recent work by Persson et al. (2010), has shown that C:P and N:P could vary by up to a factor of 20 between nutrient-replete and nutrient-limited cells. These studies have also shown that the C:N ratio can be modestly plastic due to nutrient limitation. A meta-analysis study by Hillebrand et al. (2013) highlighted the importance of growth rate in determining elemental stoichiometry and showed that both C:P and N:P ratios decrease with the increasing growth rate. Yvon-Durocher et al. (2015) investigated the role of temperature in modulating $\mathrm{C}:\mathrm{N}:\mathrm{P}$. Although their dataset was limited to studies conducted prior to 1996, they have shown a statistically significant relationship between C:P and temperature increase. MacIntyre et al. (2002) and Thrane et al. (2016) have shown that irradiance plays an important role in controlling optimal cellular C:N and N:P ratios. Most recently, Moreno and Martiny (2018) provided a comprehensive summary of how environmental conditions regulate cellular stoichiometry from a physiological perspective.

Here, we present results from a systematic literature review and subsequent meta-analysis to quantify how five key environmental drivers affect C:P and C:N ratios of marine phytoplankton. Unlike previous meta-analyses on the elemental stoichiometry of phytoplankton that strictly synthesized the effect of a single environmental driver, our study assessed the effects of five drivers, specifically for marine phytoplankton species. Importantly, we use a unique newly defined measure of effect size, a stoichiometry sensitivity factor (Tanioka and Matsumoto, 2017), which is a dimensionless parameter that relates a fractional change in P:C or N:C to a fractional change in a particular environmental driver. We compute the effect size for each driver–stoichiometry pair from independent studies and subsequently determine the weighted mean effect size for P:C and N:C ratios. Further, we compute the mean effect size within different subgroups of moderators such as plankton types and growth conditions to detect any systematic heterogeneity between those subgroups.

2 Materials and methods

2.1 Bibliographic search and screening

We systematically screened peer-reviewed publications on monoculture laboratory experiment studies that assessed the effects of dissolved inorganic phosphorus, dissolved inorganic nitrogen, dissolved iron, irradiance, and temperature on P:C and N:C ratios of marine phytoplankton. These five environmental drivers are considered to be the top drivers of the open-ocean phytoplankton group in studies (Boyd et al., 2010, 2015). Although CO₂ is another potentially important driver, we did not consider the effects of CO₂ on elemental ratios. The previous meta-analysis studies showed that no generalization could be made concerning the direction of trends in P:C or N:C ratios as a function of CO₂ concentration both in the laboratory-based experiments (Liu et al., 2010) and mesocosm/field-based experiments (Kim et al., 2018).

Firstly, we conducted a literature search using Web of Science (last accessed in February 2019) with the sequence of key terms (Table 1). This search yielded 4899 hits. We also closely inspected all the primary studies mentioned in the eight recent review papers on the elemental stoichiometry of phytoplankton in aquatic environments (Flynn et al., 2010; Geider and La Roche, 2002; Hillebrand et al., 2013; Moreno and Martiny, 2018; Persson et al., 2010; Thrane et al., 2016; Villar-Argaiz et al., 2018; Yvon-Durocher et al., 2015). The list is also augmented with six additional studies that did not appear in the literature search or the review papers but were cited elsewhere. Papers were further screened and selected to meet the following criteria. First, experiments must be carried out in controlled laboratory environments, where all the environmental factors, including temperature, photon flux density, salinity, and any other relevant conditions, are controlled. Second, all outdoor experiments, such as mesocosm or pond experiments, are excluded. Third, experiments must be conducted under unialgal/monoculture settings. However, we note that not all the experiments are carried out under strictly axenic conditions (i.e., not completely devoid of bacteria and viruses). Lastly, experiments must be conducted with replicates and must report either standard deviations or standard errors. Subsequent selection processes based on abstracts, graphs, tables, full text, and removal of duplicates led to a total of 104 journal articles (Fig. 1).

Table 1 Keyword search terms used for literature search (Web of Science, February 2019). In the search field, “TS” refers to a field tag for “topic” and “*” is a wildcard search operator.

Figure 1 Flow chart showing (1) the preliminary selection criteria and (2) the refined selection criteria used for determining s factors. Numbers (k values) correspond to the number of journal articles. See Supplement (Sect. S1) for a full list of studies included in the meta-analysis.

2.2 Data extraction

Data with means and standard deviations of P:C and N:C under varying environmental values provided by the original studies are used directly. GraphClick (Arizona Software Inc., 2010) was used to read off values from graphs when necessary. In cases where N:P and only one of either P:C or N:C is provided, the remaining ratio is determined by either multiplying or dividing by N:P. Similarly, elemental ratios are computed from the measurements of phytoplankton particulate organic carbon (POC), particulate organic nitrogen (PON), and particulate organic phosphorus (POP) when the ratios are not explicitly given in the original studies.

For nutrient (P, N, or Fe) manipulation studies, we selected two end-members (nutrient limited and nutrient replete) based on the definition given in the original studies. For batch and semicontinuous batch experiments, we compared the fractional change in initial concentrations between the nutrient-replete and nutrient-limited conditions when calculating the stoichiometry sensitivity factor (see Sect. 2.3.2). For continuous (chemostat or turbidostat) nutrient experiments, we used the difference in the inflow concentrations of the nutrient-replete and nutrient-limited cultures to determine the stoichiometry sensitivity factor. When multiple levels of concentrations are used, we selected two end-member points, one with the lowest growth rate and the other with the highest growth rate. When the growth rate was not provided in the original study, we selected two end-member values based on the highest and lowest nutrient uptake rate, chlorophyll concentration, or total concentration level with the underlying assumption that phytoplankton growth is nutrient limited within the range of nutrient levels considered.

For temperature and irradiance manipulations studies, we selected the lowest value and the optimal or saturating value that led to the maximum growth rate for phytoplankton. When the growth rate was not explicitly mentioned, we selected the lowest and the highest treatment values with the assumption that the phytoplankton is temperature or light limited within the range of values considered.

When more than two factors were manipulated in the same study, multiple experimental units are extracted if and only if each environmental driver was manipulated separately (i.e., conducted in a factorial manner). For example, we obtained a total of four experimental units from a two-by-two factorial study on temperature and nutrient: (1) comparing nutrient-limited vs. nutrient-replete treatment at low temperature; (2) same as in (1) at high temperature; (3) comparing low- vs. high-temperature response at nutrient-limited conditions; and (4) as in (3) at nutrient-replete conditions. An experimental unit refers to a controlled experiment of the same phytoplankton species between control and treatment groups, while all the other environmental factors are kept constant. If an experiment reported multiple measurements over time, only the final value was extracted.

We also extracted information on phytoplankton functional type (i.e., diatoms, coccolithophores, dinoflagellates, other eukaryotes, nondiazotrophic cyanobacteria, and diazotrophs; eukaryotes vs. prokaryotes; cold-water vs. temperate species), growth mode (i.e., batch vs. semicontinuous vs. continuous), growth phase at harvest (i.e., lag, exponential, stationary, decline), N form [${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}+{\mathrm{NH}}_{\mathrm{4}}^{+}$, N₂], and light regime (i.e., continuous vs. periodic light). Cold-water species is operationally defined if the control temperature (for P, N, Fe, or I manipulated experiments) or the maximum treatment temperature (for T manipulated experiments) was less than the threshold temperature of 10 ^∘C. Attempted but ultimately discarded moderators for subsequent analysis mainly due to the lack of sample size include salinity, the axenic nature of the culture, and the number of generations required for acclimation before the start of the experiment.

Our final dataset consists of 241 experimental units of P:C and 366 experimental units of N:C from 104 journal articles encompassing seven taxonomic phyla (Bacillariophyta, Chlorophyta, Cryptophyta, Cyanobacteria, Haptophyta, Miozoa, and Ochrophyta) and six plankton functional types (diatoms, coccolithophores, dinoflagellates, other eukaryotes, nondiazotrophic cyanobacteria, and diazotrophs), and they are available in the Zenodo data repository (https://doi.org/10.5281/zenodo.3723121; Tanioka and Matsumoto, 2020).

2.3 Statistical analysis

We used two different measures of effect size for this study. One is a commonly used natural-logarithm-transformed response ratio, ln (RR) (Hedges et al., 1999), and the other is the stoichiometry sensitivity factor (Tanioka and Matsumoto, 2017). By using two separate measures, we can give a more robust prediction on how elemental stoichiometry varies with a change in a given environmental driver. All statistical analyses were performed with R v3.5.2 (R Core Team, 2018).

2.4 Response ratio

The natural-logarithm-transformed response ratio ln (RR) of the individual experimental unit and its variance (v) were calculated following Lajeunesse (2015):

$$\begin{array}{}\text{(1)}& {\displaystyle }\ln \left(\mathrm{RR}\right)=\ln \left({\displaystyle \frac{Y_{\mathrm{t}}}{Y_{\mathrm{c}}}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{S_{\mathrm{t}}^{\mathrm{2}}}{N_{\mathrm{t}}\cdot Y_{\mathrm{t}}^{\mathrm{2}}}}-{\displaystyle \frac{S_{\mathrm{c}}^{\mathrm{2}}}{N_{\mathrm{c}}\cdot Y_{\mathrm{c}}^{\mathrm{2}}}}\right],\text{(2)}& {\displaystyle }v={\displaystyle \frac{S_{\mathrm{t}}^{\mathrm{2}}}{N_{\mathrm{t}}\cdot Y_{\mathrm{t}}^{\mathrm{2}}}}+{\displaystyle \frac{S_{\mathrm{c}}^{\mathrm{2}}}{N_{\mathrm{c}}\cdot Y_{\mathrm{c}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{S_{\mathrm{t}}^{\mathrm{4}}}{N_{\mathrm{t}}^{\mathrm{2}}\cdot Y_{\mathrm{t}}^{\mathrm{4}}}}+{\displaystyle \frac{S_{\mathrm{c}}^{\mathrm{2}}}{N_{\mathrm{c}}^{\mathrm{2}}\cdot Y_{\mathrm{c}}^{\mathrm{4}}}}\right].\end{array}$$

Y denotes mean P:C or N:C, S is the standard deviation of that mean, and N is the sample size for the treatment (subscript t) and the control (subscript c) groups. We removed any experimental unit with a studentized residual value of ln (RR) exceeding the absolute value of 3 as an outlier (Viechtbauer and Cheung, 2010).

2.4.1 Stoichiometry sensitivity factor

The second effect size is the newly defined stoichiometry sensitivity factor $s_X^Y$ (Tanioka and Matsumoto, 2017), which relates a fractional change in an elemental stoichiometry (response variable Y) to a fractional change in an environmental driver (variable X):

$$\begin{array}{}\text{(3)}& s_X^Y={\displaystyle \frac{(Y_{\mathrm{t}}-Y_{\mathrm{c}})/Y_{\mathrm{c}}}{(X_{\mathrm{t}}-X_{\mathrm{c}})/X_{\mathrm{c}}}}.\end{array}$$

We estimated the variance of $s_X^Y$ from the simple error propagation of Eq. (3) by assuming that the uncertainties associated with the environmental driver X are negligible compared to the errors associated with Y:

$$\begin{array}{}\text{(4)}& v_X^Y={\left({\displaystyle \frac{(Y_{\mathrm{t}}-Y_{\mathrm{c}})/Y_{\mathrm{c}}}{(X_{\mathrm{t}}-X_{\mathrm{c}})/X_{\mathrm{c}}}}\right)}^{\mathrm{2}}\left[{\displaystyle \frac{S_{\mathrm{t}}^{\mathrm{2}}/N_{\mathrm{t}}+S_{\mathrm{c}}^{\mathrm{2}}/N_{\mathrm{c}}}{{\left(Y_{\mathrm{t}}-Y_{\mathrm{c}}\right)}^{\mathrm{2}}}}+{\displaystyle \frac{S_{\mathrm{c}}^{\mathrm{2}}}{N_{\mathrm{c}}\cdot Y_{\mathrm{c}}^{\mathrm{2}}}}\right].\end{array}$$

In essence, the magnitude of the s factor is a measure of how sensitive Y (P:C or N:C) is to a change in stressor level X, and the sign indicates whether Y changes in the same direction as X (positive sign) or in the opposite direction to X (negative sign). The s factor allows for different kinds of response: a linear response of Y with respect to X ($s_X^Y=\mathrm{1}$), a near hyperbolic response that saturates at high X ($\mathrm{0}<s_X^Y<\mathrm{1}$), a logarithmic growth ($\mathrm{1}<s_X^Y$), a decay ($\mathrm{0}>s_X^Y$), and the null response ($s_X^Y=\mathrm{0}$). This s-factor metric is conceptually similar to the homeostasis coefficient H (Persson et al., 2010), which relates the fractional change in resource nutrient stoichiometry to the fractional change in the organism's nutrient stoichiometry.

Importantly, the advantage of using $s_X^Y$ as effect size is that its magnitude is a direct, quantitative measure of the strength of environmental drivers over the range of values examined. In contrast, ln (RR) only compares the effect of stressors without taking changes in the value of stressors into an account. Further, we can directly compare the strength of $s_X^Y$ across different pairs of X and Y as it is nondimensional. For convenience, we use the term “s factor” in the rest of this paper when describing $s_X^Y$ in a generic sense.

We used the same set of experimental units used in calculating ln (RR) to calculate s factors (i.e., any outliers are carried over). However, we did not calculate s factors for iron because the fractional change in dissolved iron concentration, often spanning multiple orders of magnitude, is substantially greater compared to the fractional change in P:C or N:C ratios, leading to an extremely low s factor. For temperature-manipulated experiments, we converted degrees Celsius into absolute temperature scale Kelvin. We used photon-flux density (PFD) measured in micromoles of photons per square meter per second (µmol photons m⁻² s⁻¹) for irradiance and micromolar (µM) for inorganic phosphorus and nitrogen experiments.

2.4.2 Meta-analysis and weighted mean responses

We calculated the weighted mean ln (RR) ($\stackrel{\mathrm{‾}}{\ln \left(\mathrm{RR}\right)}$) and s factor ($\stackrel{\mathrm{‾}}{s_X^Y}$) using the mixed-effects model with the R package metafor (Viechtbauer, 2010). The weighted mean (M) and its variance (V) were calculated as

$$\begin{array}{}\text{(5)}& {\displaystyle }M={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^k{W}_j{M}_j}{{\sum }_{j=\mathrm{1}}^k{W}_j}},\text{(6)}& {\displaystyle }V={\displaystyle \frac{\mathrm{1}}{{\sum }_{j=\mathrm{1}}^k{W}_j}},\end{array}$$

where k is the total number of experimental units, M_j is effect size (ln (RR) or $s_X^Y$) in experimental unit j, and W_j is the weighting factor, which is the inverse of the variance (Hedges et al., 1999). The 95 % confidence interval (CI) for the weighted mean was computed as

$$\begin{array}{}\text{(7)}& \mathrm{CI}=M\pm \mathrm{1.96}\times \sqrt{V}.\end{array}$$

In the subsequent sections of this paper, the values of $\stackrel{\mathrm{‾}}{\ln \left(\mathrm{RR}\right)}$ are back-transformed and represented as percent change:

$$\begin{array}{}\text{(8)}& \left(e^{\stackrel{\mathrm{‾}}{\ln \left(\mathrm{RR}\right)}}-\mathrm{1}\right)\times \mathrm{100}\mathit{\%},\end{array}$$

and they are considered statistically significant if the 95 % CIs do not overlap with zero.

2.4.3 Testing the effect of moderators

We determined the effects of moderators by the rma function of the metafor package, which is an omnibus test of between-moderator heterogeneity based on χ² distribution (Liang et al., 2020). The moderators we tested are plankton functional type (PFT), N form, growth mode, growth phase at extraction, and light regime (continuous vs. periodic). The effect of a moderator is considered significant when P value is less than 0.05. We use the weighted mean s factors in determining the effects of moderators except for iron experiments, where we used ln (RR) instead.

3 Results

Phosphate addition increases both the mean P:C (235 % [95 % CI: 169 %, 322 %]) and N:C (23 % [13 %, 34 %]) significantly (Fig. 2b). The mean stoichiometric sensitivity factor of P:C (${\mathrm{s}}_{\mathrm{P}}^{\mathrm{P}:\mathrm{C}}$) for change in phosphate is 0.21 [0.12, 0.29] (Table 2), which means that on average the P:C ratio of phytoplankton changes by 0.21 % for every 1 % increase in PO₄ concentration. The effect of phosphate on N:C is an order of magnitude smaller but also statistically significant and positively correlated (${\mathrm{s}}_{\mathrm{P}}^{\mathrm{N}:\mathrm{C}}=\mathrm{0.023}$ [0.004, 0.042]). Eukaryotic phytoplankton have significantly larger ${\mathrm{s}}_{\mathrm{P}}^{\mathrm{P}:\mathrm{C}}$ than prokaryotes (P<0.05, Fig. 3a), and the diatoms and coccolithophores especially have noticeably large $s_{\mathrm{P}}^{\mathrm{P}:\mathrm{C}}$ (Fig. S1a, Table S1 in the Supplement). In addition, phytoplankton grown under chemostat experiments have significantly larger stoichiometric sensitivity compared to those grown under batch or chemostat conditions (Fig. 3b, P<0.001). There was no between-moderator heterogeneity in ${\mathrm{s}}_{\mathrm{P}}^{\mathrm{N}:\mathrm{C}}$ (Table S1).

Table 2 Summary of the meta-analysis using the stoichiometry sensitivity factor and natural-logarithm-transformed response ratio. n, number of experimental units (numbers in bracket represent the number of outlier studies); $\stackrel{\mathrm{‾}}{s_X^Y}$, weighted mean stoichiometry sensitivity factor with environmental driver X and response variable Y; $\stackrel{\mathrm{‾}}{\ln \left(\mathrm{RR}\right)}$, weighted mean value of the natural-logarithm-transformed response ratio; ci.lb, lower boundary of 95 % CI; ci.ub, upper boundary of 95 % CI; sig., significance of the mean weighted effect size; ns, P>0.05; ^* P<0.05; ${}^{**}$ P<0.01; ${}^{***}$ P<0.001. Any experiment with a studentized residual value of ln (RR) exceeding 3 was removed as an outlier. Bold texts highlight statistically significant environmental drivers.

Figure 2 Summary plot showing weighted mean responses of P:C and N:C using the (a) stoichiometry sensitivity factor and (b) % changes between control and treatment. Numbers next to the plots in (b) correspond to the number of experimental units, and the numbers are identical in (a). Numbers in the outside column are the weighted means. ^* P<0.05; ${}^{**}$ P<0.01; ${}^{***}$ P<0.001; ns, not significant. Note that the x axis is different for temperature experiments in (a).

The response of N:C to changes in inorganic nitrogen is similar to the response of P:C to PO₄ changes where an increase in inorganic nitrogen raises N:C on average by 70 % [49 %, 93 %] (Fig. 2b), with the positive overall mean s factor $s_{\mathrm{N}}^{\mathrm{N}:\mathrm{C}}$ of 0.14 [0.08, 0.20] (Table 2). Again, eukaryotic phytoplankton have higher stoichiometric sensitivity than prokaryotes (Fig. 3a, P<0.05). Nitrogen addition does not affect the weighted mean P:C (Fig. 2). Surprisingly, however, phytoplankton grown with the culture made up of nitrate and ammonia have significantly larger $s_{\mathrm{N}}^{\mathrm{P}:\mathrm{C}}$ compared to those grown with nitrate only, with ammonia only, or under semidiazotrophic conditions (Fig. S2, Table S1). The small sample size, however, precludes us from making any firm conclusions.

Figure 3 Summary plot showing statistically significant effects of moderators. (a) Eukaryotes vs. prokaryotes, (b) growth mode, (c) growth phase at harvest, and (d) light regime. ^* P<0.05; ${}^{**}$ P<0.01; ${}^{***}$ P<0.001; ns, not significant.

An increase in iron availability does not lead to significant changes in both P:C and N:C (Fig. 2b). In addition, the effects of any moderators are not statistically significant (Table S1). Although diazotrophs that utilize N₂ as their nitrogen source display a large response compared to other PFTs (−20 % [−36 %, 1 %]) (Table S1), their stoichiometric response is not quite statistically significant.

Increase in light availability significantly decreases both P:C (−21 % [−38 %, −0.4 %]) and N:C (−18 % [−23 %, −12 %]), with overall negative s factors ($s_I^{\mathrm{P}:\mathrm{C}}=-\mathrm{0.034}$ [−0.062, −0.007], $s_I^{\mathrm{N}:\mathrm{C}}=-\mathrm{0.024}$ [−0.034, −0.013]). Although the magnitudes of both the response ratios and s factors are small compared to those of macronutrients, the responses across PFTs are consistent (Fig. S1c, f, Table S1). Phytoplankton grown under chemostat or batch conditions have significantly more negative $s_I^{\mathrm{N}:\mathrm{C}}$ compared to those grown under semicontinuous environments (Fig. 3b, P<0.01). Also, plankton grown under periodic light cycles have significantly lower ${\mathrm{s}}_I^{\mathrm{N}:\mathrm{C}}$ compared to those grown under continuous light (Fig. 3d, P<0.05).

The response of P:C to warming is significant, where on average P:C decreases by 15 % [−24 %, −5 %] with negative mean s factor of $s_T^{\mathrm{P}:\mathrm{C}}=-\mathrm{3.6}$ [−6.8, −0.4] (Fig. 2a, b). The large magnitude of the s factor compared to that of other drivers reflects the fact that the fractional change in temperature (measured in kelvins) is considerably smaller than the fractional change in P:C. There is a significant variability due to growth mode where batch culture and chemostat culture experiments respectively have more negative s factors for P:C and N:C (Fig. 3b, P<0.05). Further, phytoplankton extracted during the exponential phase have noticeably more negative s factors than those extracted during the stationary growth phase (Fig. 3c) for both P:C (P<0.001) and N:C (P<0.05). The difference in mean response s-factor ratio amongst PFTs and between cold vs. temperate species is not statistically significant (Fig. S1e, Table S1). The responses of N:C are mixed, and the weighted mean effect sizes are therefore not statistically significant.

4 Discussion

4.1 Basic framework

One of the fundamental tenets of chemical oceanography is the Redfield ratio, which implies that phytoplankton cells achieve a constant cellular $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio at the well-known molar ratio of 106 : 16 : 1 (Redfield et al., 1963). Constant $\mathrm{C}:\mathrm{N}:\mathrm{P}$ is achieved for algal cells growing under steady-state conditions, where the balance is achieved between uptake of elements and assimilation into a cellular functional pool (Berman-Frank and Dubinsky, 1999; Klausmeier et al., 2004). Under such conditions, the growth rate of all cellular constituents averaged over one generation is the same, whether it is the carbon-specific, nitrogen-specific, or phosphorus-specific growth rates (Falkowski and Raven, 2007). In the real ocean, however, balanced growth is not always achieved due to short-term and long-term changes in the physical conditions of the ocean (Moore et al., 2013; Moore and Doney, 2007). For example, the deficiency of essential nutrients limits the formation of building blocks of new cells (e.g., N for proteins, P for nucleic acids and ATP), light limitation slows carbon assimilation (i.e., making of carbohydrates and reductants), and low temperature slows down the essential cellular transport and enzymatic reactions for growth (Madigan et al., 2006). A good example of unbalanced growth is phytoplankton blooms in the spring, where the transient changes in surface temperature, irradiance, and nutrient supply rate alter the growth rate and the elemental stoichiometry of phytoplankton (Polimene et al., 2015; Talarmin et al., 2016). In addition, future environmental variabilities caused by climate change are expected to cause temporal shifts in phytoplankton $\mathrm{C}:\mathrm{N}:\mathrm{P}$ on longer timescales (Kwiatkowski et al., 2018, 2019; Tanioka and Matsumoto, 2017).

The degrees to which phytoplankton $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios are affected by environmental stresses depend both on the cellular stress response mechanisms and the magnitude of the environmental change as well as temporal variability of environmental drivers. Most types of stress responses can be divided into a stress-specific, primary response and a general secondary response (Brembu et al., 2017). The stress-specific responses are strong, robust, and consistently observed across photosynthetic organisms, while secondary responses are variable amongst different microorganisms. Primary and secondary responses are closely related to acclimation (plasticity response) and adaptation (evolutionary response), respectively. In essence, acclimation refers to the environmentally induced trait change of an organism in the absence of any genetic modification, while adaptation involves genetic changes driven by natural selection (Collins et al., 2020). Since primary responses do not involve genetic adjustment or natural selection, the responses are fast and often commonly shared amongst different marine phytoplankton. For example, changing the nutrient uptake affinity of a lineage within a generation in response to changing nutrient supply is a widely seen trait across all phytoplankton groups.

On the other hand, the secondary response depends both on the environmental condition and genotype (Brembu et al., 2017). The secondary responses take more time (usually up to a few hundred generations), and there is typically no single, unique response even when referring to a single species or functional group and a specific environmental driver (Collins et al., 2020). In the subsections below, we discuss any possible underlying cellular mechanisms responsible for producing changes in $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios (see Fig. 4 for schematic illustration).

Figure 4 Illustration of how the five environmental drivers under a typical future climate scenario affect the cellular allocation of volume between P-rich (red), N-rich (blue), and C-rich (orange) pools. The values for projected changes in C:P and C:N between 1981–2000 and 2081–2100 are given in Table 3.

4.2 Macronutrients (phosphate and nitrate)

Overall, we observe a consistent trend across all studies where P:C and N:C increases with an increase in the supply of dissolved inorganic phosphorus and nitrogen, respectively (Fig. 2). Since the changes in X:C and the supply of element X are positively related, ${\mathrm{s}}_{\mathrm{P}}^{\mathrm{P}:\mathrm{C}}$ and ${\mathrm{s}}_{\mathrm{N}}^{\mathrm{N}:\mathrm{C}}$ are both positive. Observations of phosphate (nitrate) against particulate organic matter P:C (N:C) across the global ocean indeed broadly follow this general trend (Galbraith and Martiny, 2015; Tanioka and Matsumoto, 2017).

Phytoplankton can temporally store excess nutrient intracellularly until the rate of carbon assimilation catches up to achieve steady-state balanced growth. Excess phosphorus, for example, can be stored mainly as polyphosphate (Dyhrman, 2016), and excess nitrate can be stored primarily as protein and free amino acids (Liefer et al., 2019; Sterner and Elser, 2002). Phytoplankton can consume these internal stores of nutrients (e.g., polyphosphates under P limitation) while maintaining the same level of carbon fixation, when the uptake of the nutrients does not meet their demand for growth (Cembella et al., 1984). Also, phytoplankton can reduce their number of ribosomes and RNA content under P limitation as RNA typically accounts for 50 % of nonstorage phosphorus (Hessen et al., 2017; Lin et al., 2016). Similarly, cells can reduce the synthesis of N-rich protein content under N limitation, resulting in a lower N:C ratio (Grosse et al., 2017; Liefer et al., 2019). These transient processes controlling the intracellular content of P or N (but not C content as much) likely result in a positive correlation between P:C and N:C with macronutrient concentrations.

Although $s_{\mathrm{P}}^{\mathrm{P}:\mathrm{C}}$ and $s_{\mathrm{N}}^{\mathrm{N}:\mathrm{C}}$ are consistently positive across all the studies, they are noticeably higher for eukaryotic phytoplankton than for prokaryotes (Fig. 3a). There are several hypotheses for explaining this trend. One of the most plausible explanations is related to the cell size and storage capacity difference amongst phytoplankton groups (Edwards et al., 2012; Lomas et al., 2014). Since eukaryotes are generally larger and possess more storage capacity, they are capable of greater luxury uptake and accumulation of internal P and N reserves when the nutrient is in excess (Talmy et al., 2014; Tozzi et al., 2004). When nutrients are scarce, the large cell size of eukaryotes allows them to increase their carbon content considerably by accumulating excess carbon as polysaccharides and lipids (Liefer et al., 2019; Lin et al., 2016). Another plausible hypothesis concerns variability in the acclimation/adaptation strategy at the genetic level (Dyhrman, 2016). Recent studies suggest that different phytoplankton groups exhibit different levels of transcriptional responsiveness and have different strategies for using nitrate (Lampe et al., 2019) and phosphate (Martiny et al., 2019). For example, diatoms have superior abilities to uptake and store nutrients by being able to quickly regulate their gene expression patterns required for nutrient uptake compared to other phytoplankton groups (Cáceres et al., 2019; Lampe et al., 2018, 2019). These hypotheses provide plausible explanations for why eukaryotes have elevated stoichiometry sensitivity to macronutrients compared to prokaryotes.

4.3 Iron

Iron is used in key biochemical processes such as electron transport, respiration, protein synthesis, and N fixation (Marchetti and Maldonado, 2016; Twining and Baines, 2013). Many of the iron-dependent processes are required for harvesting energy and for synthesizing biochemical intermediates (Price, 2005). As energy acquisition is equivalent to light acquisition in phototrophs, it makes sense that percent changes in stoichiometry for iron are similar in sign and magnitude to that of light (Fig. 2b). Although the effect of increasing iron on N:C is similar in sign and magnitude to that of light, increasing iron availability does not lead to a significant change in mean N:C (Fig. 2b). This suggests smaller-than-expected changes in the carbon or the nitrogen content (e.g., compounds such as porphyrin and phycobiliprotein that are essential for light harvesting) under Fe limitation (Falkowski and Raven, 2007; Twining and Baines, 2013). Alternatively, Fe availability may be affecting cellular C, N, and P more or less proportionally for all phytoplankton, leading to constant P:C and N:C (Greene et al., 1991; van Oijen et al., 2004; La Roche et al., 1993; Takeda, 1998). We also did not find noticeable heterogeneities in P:C and N:C amongst different moderators. Yet a number of laboratory studies, particularly those of picocyanobacteria (Prochlorococcus and Synechococcus), display significant effects of iron on $\mathrm{C}:\mathrm{N}:\mathrm{P}$ (e.g., Cunningham and John, 2017). Despite their ecological importance (Biller et al., 2015; Flombaum et al., 2013), these taxa are understudied compared to diazotrophic cyanobacteria and diatoms. Future studies could focus on these picocyanobacteria and combine cellular $\mathrm{C}:\mathrm{N}:\mathrm{P}$ information with other measures of phytoplankton physiology (e.g., chlorophyll fluorescence, F_v∕F_m ratio) to provide a more coherent, mechanistic picture of how changes in iron availability affect their physiology.

4.4 Irradiance

Light availability affects the photoacclimation strategy of phytoplankton and, subsequently, the cellular allocation of volume between the N-rich light-harvesting apparatus, P-rich biosynthetic apparatus, and C-rich energy storage reserves (Falkowski and LaRoche, 1991; Moreno and Martiny, 2018). At a fixed growth rate, high irradiance should downregulate the production of N-rich light-harvesting proteins and pigments to minimize the risk of photooxidative stress. Excess carbon fixed under high-irradiance conditions is stored as C-rich storage compounds such as lipids and polysaccharides (Berman-Frank and Dubinsky, 1999). As a result, N:C is expected to decrease under high light. In contrast, under low-light conditions, the macromolecular composition should favor the N-rich light-harvesting apparatus over C-rich storage reserves, thus elevating N:C. This line of reasoning would predict a negative relationship for the effect of irradiance increase on N:C, which is borne out in our meta-analysis (Fig. 2). Similarly, the P quota should be affected by a change in irradiance (Moreno and Martiny, 2018). P:C is expected to decrease at the increased light level because the total supply of inorganic phosphorus will not be able to keep up with the increase in photosynthetic carbon fixation, leading to a decoupled uptake of C and P (Hessen et al., 2002, 2008). Conversely, P:C is expected to increase at lower irradiance because carbon fixation decreases while phosphorus uptake remains constant (Urabe and Sterner, 1996).

The magnitudes of the weighted mean s factors for both P:C and N:C, however, are small, and the heterogeneity amongst PFTs is not discernible. This result agrees with a previous study that compiled experimental data prior to 1997 (MacIntyre et al., 2002). It is possible, however, that s factors obtained in our meta-analysis are underestimated as several factors may mute the effect of irradiance on the N:C ratio of phytoplankton. For example, an increase in nitrogen requirement for Rubisco (Li et al., 2015) and nutrient uptake machinery (Ågren, 2004) at high irradiance could partly offset the reduction in N content resulting from the downregulation of the light-harvesting apparatus. In addition, multiple studies have noted an increase in the protein demand (e.g., D1 protein) for repairing the damaged light-harvesting apparatus at high irradiance (Demmig-Adams and Adams, 1992; Li et al., 2015; Talmy et al., 2013), which also works in favor of stabilizing N content. Furthermore, we may have underestimated our s factor if the high end-member irradiance were above the optimal light level. This last reason is a fundamental limitation of s-factor determination as most studies we selected do not measure the actual optimal irradiance but simply report an arbitrary value that is either “high” or “light replete”.

Interestingly, we observed larger stoichiometric shifts in nutrient-replete batch and chemostat culture experiments compared to those conducted under semicontinuous settings (Fig. 3b). In addition, we found that experiments conducted under periodic daily light cycles have larger negative s factors compared to those experiments carried out under continuous light (Fig. 3d). These results are consistent with the global observation (Martiny et al., 2013a) and model studies (Arteaga et al., 2014; Talmy et al., 2014, 2016) which have shown that both the magnitude and temporal variability of N:C is higher in the nutrient-rich, light-limited polar regions than in the light-replete subtropics.

4.5 Temperature

We found that the P:C ratio decreases as temperature increases, while N:C remains relatively unchanged. Our result is consistent with a previous meta-analysis (Yvon-Durocher et al., 2015) that showed a decrease in phytoplankton P:C with temperature increase under laboratory and field settings. Moreover, our study and the study by Yvon-Durocher et al. (2015) support the idea that P:C is more flexible than N:C with respect to change in temperature, which suggests that intracellular P content is more sensitive to change in temperature than intracellular N content. Although the underlying mechanism for explaining lower P:C at higher temperature is not fully understood, there are currently three main hypotheses (Paul et al., 2015): (1) increase in metabolic stimulation of inorganic carbon uptake over phosphorus uptake; (2) increase in nutrient use efficiency which enables greater carbon fixation for given nutrient availability; and (3) translation compensation theory, which predicts that less P-rich ribosomes are required for protein synthesis and growth as the translation process becomes kinetically more efficient (McKew et al., 2015; Toseland et al., 2013; Woods et al., 2003; Xu et al., 2014; Zhu et al., 2017).

Differences in s factors amongst PFTs were not statistically significant, and none of the PFT displayed a statistically significant response in isolation. In other words, we did not see any PFT-specific adaptive/evolutionary response to warming (Schaum et al., 2018; Taucher et al., 2015). However, we observed noticeable variability due to the difference in culture growth mode (Fig. 3b) and the growth phase at extraction (Fig. 3c). The latter factor is particularly noticeable for P:C, where phytoplankton extracted during the nutrient-replete exponential growth phase have significantly more negative stoichiometric flexibility with a larger magnitude compared to those extracted during the nutrient-deplete stationary phase. This is consistent with multiple recent studies which suggest that the effect of temperature on growth and metabolic rates is greater when plankton are not nutrient or light limited (Aranguren-Gassis et al., 2019; Marañón et al., 2018; Roleda et al., 2013). This leads us to hypothesize that change in the P:C ratio due to ongoing warming will be more noticeable in the nutrient-rich polar regions, especially given the fact that temperature is already increasing at a startling rate due to polar amplification (Post et al., 2019).

4.6 Limitations and caveats

In the real ocean, none of the environmental changes discussed will likely occur in isolation because changes in irradiance, temperature, and nutrient availability are often linked. For example, an increase in sea surface temperature enhances the vertical stratification of the water column, which leads to greater levels of irradiance and nutrient limitation for phytoplankton trapped in a more shallow mixed layer (Boyd et al., 2015; Hutchins and Fu, 2017). Indeed, a meta-analysis on the pair-wise effects of environmental drivers on the elemental stoichiometry of phytoplankton has shown that the interactions of two environmental stressors can impose predominantly nonadditive effects to $\mathrm{C}:\mathrm{N}:\mathrm{P}$ of phytoplankton so that the overall effect of multiple stressors is more than simply the sum of its parts (Villar-Argaiz et al., 2018). In addition to the individual phytoplankton stoichiometry, the bulk organic matter stoichiometry also reflects the phytoplankton community composition (Bonachela et al., 2016; Weber and Deutsch, 2010) as well as the stoichiometry of detrital material. Processes such as decomposition (Karl and Dobbs, 1998; Verity et al., 2000; Zakem and Levine, 2019), viral shunt (Jover et al., 2014), and preferential remineralization of phytoplankton macromolecules (Frigstad et al., 2011; Grabowski et al., 2019; Kreus et al., 2015) can also decouple phytoplankton $\mathrm{C}:\mathrm{N}:\mathrm{P}$ from the bulk organic matter $\mathrm{C}:\mathrm{N}:\mathrm{P}$.

Table 3 Projected changes in C:P (molar) and C:N (molar) between 1981–2000 and 2081–2100 given model-based projected changes in environmental drivers from Boyd et al. (2015). Changes in C:N and C:P are calculated separately for each driver, with s factors from Table 2 combined with reference $\mathrm{C}:\mathrm{N}:\mathrm{P}$ of 146 : 20 : 1, a global biomass-weighted mean ratio of particulate organic matter (Martiny et al., 2013b). Ranges are derived from propagating uncertainties for the weighted mean s factors in Table 2. We used Eq. (9) in the main text for estimating the combined effect of multiple drivers.

4.7 Implications for global ocean biogeochemistry

Recent global biogeochemical models are starting to incorporate a more realistic representation of plankton physiology, which includes flexible phytoplankton $\mathrm{C}:\mathrm{N}:\mathrm{P}$ (e.g., Buchanan et al., 2018). Modeling studies with flexible phytoplankton stoichiometry have demonstrated that proliferation of C-rich phytoplankton under future climate scenarios has the potential to buffer expected future decline in carbon export and net primary productivity caused by increased stratification (Kwiatkowski et al., 2018; Moreno et al., 2018; Tanioka and Matsumoto, 2017). This buffering effect cannot be simulated by biogeochemical models with fixed phytoplankton $\mathrm{C}:\mathrm{N}:\mathrm{P}$.

One way to model the dependencies of multiple environmental drivers (e.g., P, N, irradiance, and temperature) on $\mathrm{C}:\mathrm{N}:\mathrm{P}$ of marine phytoplankton is the power-law formulation by Tanioka and Matsumoto (2017):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\left[X:\mathrm{C}\right]& ={\left[X:\mathrm{C}\right]}_{\mathrm{0}}{\left({{\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}\right]}{\left[{\mathrm{PO}}_{\mathrm{4}}\right]}}}_{\mathrm{0}}\right)}^{s_{{\mathrm{PO}}_{\mathrm{4}}}^{X:\mathrm{C}}}{\left({{\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}}}_{\mathrm{0}}\right)}^{s_{{\mathrm{NO}}_{\mathrm{3}}}^{X:\mathrm{C}}}\\ & \cdot {\left({{\displaystyle \frac{I}{I}}}_{\mathrm{0}}\right)}^{{\mathrm{s}}_I^{X:\mathrm{C}}}{\left({{\displaystyle \frac{T}{T}}}_{\mathrm{0}}\right)}^{s_T^{X:\mathrm{C}}}(X=\text{P or N})\end{array},\end{array}$$

where subscript “0” indicates reference values. The s factors obtained from this meta-analysis are the exponents of Eq. (9) for different environmental drivers. Within the context of the power-law formulation, our results would indicate, for example, that eukaryotic phytoplankton would have the largest plasticity in P:C and N:C compared to prokaryotes with respect to the change in nutrient availability. Under future warming, high s factors of eukaryotes may thus play an important role in buffering the expected future decline in carbon export and net primary productivity (Kemp and Villareal, 2013).

We can give a first-order estimate of how much the elemental stoichiometry of marine phytoplankton may change in the future using Eq. (9) given a typical projection of the change in the key environmental drivers (Table 3; Fig. 4). Global climate models generally predict a decline in macronutrients and an increase in temperature and irradiance as a result of surface warming, increased vertical stratification, and reduced mixed layer depth (Bopp et al., 2013; Boyd et al., 2015). With large projected declines in macronutrients (−28.0 % for phosphate, −18.7 % for nitrate) we can predict an increase in C:P and C:N by ∼10 units (molar ratio) and ∼0.2 units, respectively, assuming the mean biomass-weighted particulate organic matter $\mathrm{C}:\mathrm{N}:\mathrm{P}$ of $\mathrm{146}:\mathrm{20}:\mathrm{1}$ as the present-day value (Martiny et al., 2013b). Further increase in C:P is expected due to the temperature increase of around 1 % (∼3 K). The total C:P change ranges from +6 to +25, considering all the uncertainties associated with the s factors. For C:N, we estimate an overall increase by 0.1–0.4 units largely driven by a decrease in nitrogen availability. The effect of change in irradiance is noticeably smaller (Table 3). In summary, this simple calculation highlights a potentially large shift for $\mathrm{C}:\mathrm{N}:\mathrm{P}$, whose change is predominantly driven by a reduction in macronutrients and temperature increase.

5 Conclusions

Our meta-analysis represents an important bottom-up approach in predicting how elemental stoichiometry of phytoplankton may evolve with climate change. We conclude that macronutrient availability is the most significant and shared environmental driver of $\mathrm{C}:\mathrm{N}:\mathrm{P}$. Changes in $\mathrm{C}:\mathrm{N}:\mathrm{P}$ by macronutrients are driven by primary/plasticity responses commonly shared across phytoplankton. Our analysis shows that eukaryotic phytoplankton have higher stoichiometric plasticity compared to prokaryotes. Eukaryotes' large stoichiometric flexibility and high intrinsic growth rate can explain their unexpectedly high diversity (Malviya et al., 2016) and large contribution to carbon export globally, even in oligotrophic regions (Agusti et al., 2015; Nelson and Brzezinski, 1997). The effects of temperature on C:P are also significant, suggesting that a future ocean with elevated temperature and increased stratification will favor the production of carbon-rich organic matter. Future laboratory-based studies exploring how the multiple environmental drivers interactively alter the elemental composition of phytoplankton would be needed for a complete understanding. In addition, a further investigation on how a change in environmental drivers affects the stoichiometry of heterotrophs and zooplankton will be useful in filling the gaps to gain more mechanistic views on how these drivers affect the whole marine ecosystem.