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

The elemental stoichiometry of marine phytoplankton plays a critical role in the global biogeochemical cycle through its impacts 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 10 N:C ratios of marine phytoplankton to five major drivers (inorganic phosphorus, inorganic nitrogen, inorganic iron, irradiance, and temperature) by meta-analysis of laboratory experimental data across 366 experiments from 104 journal articles. Our results show that the response of the 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 15 the changes in macronutrients compared to prokaryotes, possibly due to their larger cell size and their abilities to quickly regulate their gene expression patterns required for nutrient uptake. The effect of irradiance was significant and constant across all studies where an increase in irradiance decreased both P:C and N:C. The response to temperature changes was mixed depending on the culture growth mode and the growth phase of phytoplankton at the time of harvest but the weighted mean P:C ratio decreased 20 significantly with warming. Along with other oceanographic conditions of the subtropical gyres (e.g., low macronutrient availability), elevated temperature may explain why P:C is consistently low in subtropical oceans. Iron addition did not systematically change neither 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 25 productivity.


Introduction
Elemental stoichiometry of biological production in the surface ocean plays a crucial role in cycling of elements in the global ocean. The elemental ratio between carbon and the key limiting macronutrients, nitrogen (N) and phosphorus (P), in exported organic matter expressed in terms of C:N:P ratio helps 30 determine how much atmospheric carbon is sequestered in the deep ocean with respect to the availability of limiting nutrients. On geologic timescale, N:P ratio reflects the relative availability of nitrate with respect to phosphate, both of which are externally supplied from atmosphere via nitrogen-fixation and/or continents via river supply and lost by denitrification and burial (Broecker, 1982;Lenton and Watson, environmental drivers. Individual studies employ different sets of statistical analyses to characterize effects of 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. 60 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 the more recent work by Persson et al. (2010) have shown that C:P and N:P could vary up to a factor of 20 between nutrient-replete and nutrient-limited cells. These studies have also shown that C:N ratio is plastic due to 65 nutrient limitation. A meta-analysis study by Hillebrand et al. (2013) highlighted the importance of growth rate in determining elemental stoichiometry showed that both C:P and N:P ratios decrease with increasing growth rate. Yvon-Durocher et al. (2015) investigated the role of temperature in modulating C:N: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 70 moderators such as plankton types and growth conditions for detecting any systematic heterogeneity between those subgroups.

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 open-ocean phytoplankton group in studies (Boyd et al., 90 2010(Boyd et al., 90 , 2015. Although CO2 is another potentially important driver, we did not consider the effects of CO2 on elemental ratios as previous meta-analysis studies showed that no generalization can be made with respect to the direction of trends in P:C or N:C ratios as a function of CO2 concentration both in the laboratory-bases 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) 95 with the sequence of key terms (Table 1). This search yielded 4899 hits. We also closely inspected all the primary studies mentioned in the 8 recent review papers including meta-analyses studies on elemental stoichiometry of phytoplankton in aquatic environment (Flynn et al., 2010;Geider and La Roche, 2002;Hillebrand et al., 2013;Persson et al., 2010;Thrane et al., 2016;Villar-Argaiz et al., 2018;Yvon-Durocher et al., 2015). The list is also augmented with data from additional six studies 100 that did not appear in the literature search or in the review papers but were cited elsewhere. Papers were further screened and selected to meet the following criteria: (1) experiments must be carried out in the controlled laboratory environments, where all the environmental factors including temperature, photon flux density, salinity, and any other relevant conditions are controlled; (2) all outdoor experiments such as mesocosm or pond experiments are excluded; (3) experiments must be conducted under 105 unialgal/monoculture settings. However, we note that not all the experiments are carried under strictly axenic condition (i.e., not completely devoid of bacteria and virus); and (4) experiments must be conducted with replicates and must report either standard deviation or standard error. Subsequent selection processes based on abstracts, graphs, tables, and full text, and removal of duplicates led to a total of 104 journal articles (Fig. 1). 110

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, 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 115 computed from the measurements of phytoplankton POC, PON, and 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 semi-continuous batch experiments, we compared fractional change in initial concentrations between nutrient replete and 120 limited conditions when calculating stoichiometry sensitivity factor (see section 2.3.2). For continuous (chemostat or turbidostat) nutrient experiments, we used difference in the inflow concentrations of the nutrient replete and limited cultures to determine 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 highest growth rate. When the growth rate was not provided in the original study, we selected 125 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 growth rate 130 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 extracted total of 4 experimental units from a 2-by-2 factorial study 6 on temperature and nutrient: (1) comparing nutrient limited vs. replete treatment at low temperature; (2) same as in (1) at high temperature; (3) comparing low vs. high temperature response at nutrient limited condition; and (4) as in (3)  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, axenic nature 150 of the culture, and the number of generations required for acclimation before the start of the experiment. 7 prediction on how elemental stoichiometry varies with a change in given environmental driver. All statistical analyses were performed with R v3.5.2 (R Core Team, 2018).

Response ratios
The natural logarithm-transformed response ratios ln(RR) of individual experimental unit and its variance (v) was calculated following Lajeunesse (2015): (1) (2) 170 Y denotes mean P:C or N:C, S 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).

Stoichiometry sensitivity factor
The second effect size is the newly defined stoichiometry sensitivity factor s 9 : ( In essence, the magnitude of s-factor is a measure of how sensitive Y (P:C or N:C) is to a change 185 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 9 : = 1), a near hyperbolic response that saturates at high X (0 < s 9 : < 1), a logarithmic growth (1 < s 9 : ), a decay (0 > s 9 : ), and the null response (s 9 : = 0). This s-factor metric is conceptually similar to the homeostasis coefficient H (Persson et al., 2010), which relates 190 fractional change in resource nutrient stoichiometry to fractional change in organism's nutrient stoichiometry.
Importantly, an advantage of using s 9 : as effect size is that its magnitude is a direct, quantitative measure of the strength of environmental driver over the range of values examined. In contrast, ln(RR) only compares the effect of stressor on two end point values (control and treatment) without taking 195 changes in the stressor into an account. Further, we can directly compare the strength of s 9 : across different pairs of X and Y as it is non-dimensional. For convenience, we use the term "s-factor" in the rest of this paper when describing s 9 : 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 200 change in dissolved iron concentration, often spanning multiple orders of magnitude, are substantially larger compared to the fractional change in P:C or N:C ratios leading to extremely low s-factor. For temperature-manipulated experiments, we converted degrees Celsius into absolute temperature scale Kelvin. We used photon-flux density (PSD) measured in µmol photons m -2 s -1 for irradiance and µM for inorganic phosphorus and nitrogen experiments. 205

Meta-analysis and weighted mean responses
We calculated weighted mean ln(RR) (ln ( ) GGGGGGGGGG and s-factor (s 9 = ± 1.96 × √ (7) 215 In the subsequent sections of this paper, the values of ln ( ) GGGGGGGGGG are back-transformed and represented as percent change: and considered statistically significant if 95% CIs do not overlap with zero. 220

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

Results
Phosphate addition increases both the mean P:C (235% [95% CI: 169%, 322%]) and N:C (23% [13%, 34%]) significantly (Fig. 2b). Mean stoichiometric sensitivity factor of P:C (s P P:C ) with respect to change 230 in phosphate is 0.21 [0.12, 0.29] ( Table 2) which means that on average P:C ratio of phytoplankton changes by 0.21% for every 1% increase in PO4 concentration. The effect of phosphate on N:C is an order of magnitude smaller but also statistically significant and positively correlated (s P N:C = 0.023 [0.004, 0.042]). Eukaryotic phytoplankton have significantly larger s P P:C than prokaryotes (P < 0.05, Fig. 3a) and the diatoms and coccolithophores especially have noticeably large s P P:C (Fig. S1a, Table S1). In addition, 235 phytoplankton grown under chemostat experiments have significantly larger stoichiometric sensitivity compared to those grown under batch or chemostat condition ( Fig. 3b, P < 0.001). There was no betweenmoderator heterogeneity in s P N:C (Table S1).
The response of N:C to changes in inorganic nitrogen is similar to the response of P:C to PO4 changes where an increase in inorganic nitrogen raises N:C on average by 70% [49%, 93%] ( Fig. 2b) with 240 the positive overall mean s-factor s N N:C of 0.14 [0.08, 0.20] (Table 2). Again, eukaryotic phytoplankton have higher stoichiometric sensitivity than prokaryotes ( Fig. 3a, P < 0.01). Nitrogen addition does not affect the weighted mean P:C (Fig. 2). Surprisingly however, phytoplankton grown with the culture made of up nitrate and ammonia have significantly larger s N P:C compared to those grown with nitrate only, ammonia only, or those under semi-diazotrophic condition (Fig. S2, Table S1). The small sample size 245 however precludes us from making any firm conclusions.
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 moderatos are not statistically significant (Table S1). Although diazotrophs that utilize N2 as its nitrogen source have significantly large response compared to other PFTs (-20% [-36%, 1%]) (Table S1), their stoichiometric response is not quite statistically significant. 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 P:C = -3.6 [-6.8, -0.4]) ( Fig. 2a, b). The large magnitude of s-factor compared to that of other drivers reflects the fact that the fractional change in temperature (measured in 260 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 sfactors for P:C and N:C (Fig. 3b, P < 0.05). In addition, phytoplankton extracted during exponential have noticeably more negative s-factors than those extracted during 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 265 between cold vs. temperate species is not statistically significant (Fig. S1e, Table S1). Response of N:C is mixed and the weighted mean effect sizes are therefore not statistically significant.

Basic framework
One of the fundamental tenets of chemical oceanography is the Redfield Ratio, which implies that 270 phytoplankton cells achieve a constant cellular C:N:P ratio at the well-known molar ratio of 106:16:1 (Redfield et al., 1963). Constant C:N:P is achieved for algal cells growing under steady state conditions where the balance is achieved between uptake of elements and assimilation into 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, 275 nitrogen (protein)-specific, or phosphorus (DNA)-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 physical conditions of 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 280 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 bloom in the spring where the transient changes in surface temperature, irradiance and nutrient supply rate alter the growth rate and 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 285 shift in phytoplankton C:N:P on longer timescales (Kwiatkowski et al., 2018b(Kwiatkowski et al., , 2019Tanioka and Matsumoto, 2017).
The degrees to which phytoplankton C:N:P ratios are affected by 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, 290 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 organisms. Primary and secondary responses are closely related to acclimation (plasticity response) and adaptation (evolutionary response) respectively. In essence, acclimation refers to environmentally induced trait change of an organism in the absence of any genetic 12 change, 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 commonly seen trait across all phytoplankton groups. On the other hand, secondary response depends both on the 300 environmental condition and genotype (Brembu et al., 2017). The secondary responses take longer time (usually up to 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 C:N:P ratios (see Fig. 4 for schematic illustration). 305

Macronutrients (Phosphate and Nitrate)
Overall, we observe a consistent trend across all studies where P:C and N:C increases with 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, s P P:C and s N N:C are both positive. Observations of 310 phosphate/nitrate against particulate organic matter P:C and 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 315 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 its demand for growth (Cembella et al., 1984). In addition, phytoplankton can reduce their number of ribosomes and RNA content under P limitation as RNA typically accounts for 50% of non-storage phosphorus (Hessen et al., 2017;Lin et al., 320 2016) which would conserve phosphorus for other uses in a cell, resulting in lower P:C ratios. Similarly, cells can reduce synthesis of N-rich protein content under N limitation resulting in 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 positive correlation between P:C and N:C with macronutrient concentrations. 325 Although s P P:C and s N N: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 hypotheses is related to the 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 330 internal P and N reserves when the nutrient is in excess (Talmy et al., 2014;Tozzi et al., 2004). When 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., 2018Lampe et al., , 2019. These hypotheses provide plausible explanations for why eukaryotes have elevated stoichiometry sensitivity to macronutrients 340 compared to prokaryotes.

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 345 processes are required for harvesting energy and biochemical intermediates. As energy acquisition is equivalent to light acquisition in phototrophs, it makes sense that % changes in stoichiometry for iron are similar in sign and magnitude as for light (Fig. 2b). Although the effect of increasing iron on N:C is similar in sign and magnitude to that of light, we found unlike irradiance increase that increasing iron availability does not lead to a significant change in mean N:C (Figure 2b). This suggests smaller than 350 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 355 and N:C amongst different moderators. In the future study, we could combine cellular C:N:P information with other measures of phytoplankton physiology (e.g., chlorophyll fluorescence, Fv/Fm ratio) in order to provide a more coherent, mechanistic picture on how changes in iron availability affect their physiology.

Irradiance
Light availability affects the photoacclimation strategy of phytoplankton and subsequently the cellular allocation of volume between N-rich light-harvesting apparatus, P-rich biosynthetic apparatus, and C-rich energy storage reserves (Falkowski and LaRoche, 1991; . At a fixed growth rate, high irradiance should downregulate production of N-rich light harvesting proteins and pigments in 365 order to minimize the risk of photooxidative stress. Excess carbon fixed under high irradiance condition 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 condition, macromolecular composition should favor N-rich light harvesting apparatus over C-rich storage reserves, thus elevating N:C. This line of reasoning would predict negative relationship for the effect of irradiance 370 increase on N:C, which is borne out in our meta-analysis (Fig. 2).
Similarly, P quota should be affected by change in irradiance if P is the main limiting nutrient . Under P limitation, P:C is expected to decrease at 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 decoupled uptake of C and P (Hessen et al., 2002(Hessen et al., , 2008. 375 Conversely, P:C is expected to increase at lower irradiance because carbon fixation decreases while phosphorus uptake remains constant (Urabe and Sterner, 1996). As we did observe such P:C responses with statistically significant negative s-factor (Fig. 2), we can infer that most of the experiments were likely to have been P-limited, although such information is not necessarily given in the original studies.
The magnitude of the weighted mean s-factors for both P:C and N:C however are small and the 380 heterogeneity amongst PFTs are not discernible. This result agrees with a previous study which 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 there are several factors that may mute the effect of irradiance on N:C ratio of phytoplankton. For example, increase in nitrogen requirement for Rubisco (Li et al., 2015) and nutrient uptake machinery (Ågren, 2004) at high irradiance could be partly offset the reduction in N 385 content resulting from the down regulation of light harvesting apparatus. In addition, multiple studies have noted increase in the protein demand (e.g., D1 protein) for repairing 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 is a fundamental limitation of s-390 factor determination as the original studies do not measure the true optimal irradiance across the range of irradiance values 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 compared to those cultures conducted under semi-continuous setting (Fig. 3b). In addition, we found that experiments conducted under periodic daily light cycle have larger negative s-factors compared to 395 those experiments carried under continuous light (Fig. 3d). This is consistent with the global observation (Martiny et al., 2013a) and model studies (Arteaga et al., 2014;Talmy et al., 2014Talmy et al., , 2016 which have shown that both the magnitude and temporal variability of N:C of phytoplankton are higher in the nutrientrich, light-limited polar/subpolar regions than in the light-replete subtropics.

Temperature
We found that 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 decrease in phytoplankton P:C under both laboratory and field settings. Moreover, our study and the study by Yvon-temperature, which suggest that intracellular P content is more sensitive to change in temperature than intracellular N content.
Although the underling 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 410 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). 415 Differences in s-factors amongst PFTs was not statistically significant and none of the PFT displayed 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 growth phase at extraction (Fig. 3c). The latter factor is particularly noticeable for P:C, where phytoplankton extracted 420 during nutrient-replete exponential growth phase have significantly more negative stoichiometric flexibility with larger magnitude compared to those extracted during nutrient-deplete stationary phase. This is consistent with multiple recent studies which suggest that the effect of temperature on growth and metabolic rates are greater when plankton are not nutrient and/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 P:C ratio due 425 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).

Limitations and caveats
In the real ocean, none of the environmental changes discussed will likely occur in isolation because 430 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 elemental stoichiometry of phytoplankton has shown that the interactions of two environmental 435 stressors can impose predominantly non-additive effects to C:N: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., 440 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 C:N:P from the bulk organic matter C:N:P.

Implications for global ocean biogeochemistry 445
Recent global biogeochemical models are starting to incorporate a more realistic representation of plankton physiology, which includes flexible phytoplankton C:N:P (e.g., Buchanan et al., 2018).
Modeling studies with flexible phytoplankton stoichiometry have demonstrated that proliferation of Crich phytoplankton under future climate scenario has the potential to buffer expected future decline in carbon export and net primary productivity caused by increased stratification (Kwiatkowski et al., 2018a;450 Moreno et al., 2018;Tanioka and Matsumoto, 2017). This buffering effect cannot be simulated by biogeochemical models with fixed phytoplankton C:N:P. One way to model the dependencies of multiple environmental drivers (e.g., P, N, irradiance, and temperature) on C:N:P of marine phytoplankton is the power-law formulation by Tanioka and Matsumoto where subscript "0" indicates reference values. The s-factors obtained from this meta-analysis are the exponents of equation (9) for different PFTs. 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 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 equation (9) given a typical projection of the change in the key 465 environmental drivers and the estimates of the s-factors (Table 3; Fig. 4). Global climate models generally predict a decline in macronutrients and 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 increase in C:P and C:N by ~10 units (molar ratio) and ~0.2 units, respectively, assuming the 470 mean biomass-weighted particulate organic matter C:N:P of 146:20:1 as the present-day value (Martiny et al., 2013b). Further increase in C:P is expected due to temperature increase of around 1% (~3K). The total C:P change ranges from +6 ~ +25 taking into account all the uncertainties associated with the sfactors. For C:N, we estimate an overall increase by 0.1~0.4 units largely driven by decrease in nitrogen availability. The effect of change in irradiance is noticeably smaller (Table 3). In summary, this simple 475 calculation highlights potentially a large shift for C:N:P, whose change is predominantly driven by reduction in macronutrients and temperature increase.

Conclusions
Our meta-analysis represents an important bottom-up approach in predicting how elemental stoichiometry 480 of phytoplankton may evolve with the climate change. We conclude that macronutrient availability is the most significant and shared environmental driver of C:N:P. Changes in C:N: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 485 (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 is also significant suggesting that future ocean with elevated temperature and increased stratification will favor production of carbon-rich organic matter. Future laboratory-based studies focused on exploring the effects of multiple environmental drivers would interactively alter the elemental composition of phytoplankton would be 490 needed for a complete understanding. In addition, a further investigation on how change in environmental drivers affect 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.

495
Data availability: All the data and codes used in the meta-analysis are available in Zenodo data repository (http://doi.org/10.5281/zenodo.3723121).
Author contributions: TT and KM designed the study. TT carried out the literature review, data selection, analysis, and created figures. Both TT and KM wrote the manuscript. 500 Competing interests: The authors declare no conflict of interest.
Acknowledgements: This research was supported by a grant from the US National Science Foundation (OCE-1827948). TT acknowledges support from University of Minnesota Doctoral Dissertation 505 Fellowship. KM acknowledges sabbatical support by the Leverhulme Trust Visiting Professorship and the University of Oxford. We thank Carolyn Bishoff, Julia Kelly, and Amy Riegelman from University of Minnesota Library for helping out literature search and data selection. We also thank James Cotner for providing us feedback on the manuscript. gradients reveal diverse acclimation responses in phytoplankton phosphate uptake, ISME J., 13(11), 2834-2845, doi:10.1038/s41396-019-0473-1, 2019 The Utilization of Inorganic and Organic Phosphorous Compounds as Nutrients by Eukaryotic Microalgae: A Multidisciplinary Perspective: Part 2, CRC Crit. Rev. Microbiol., 11(1), 13-81, doi:10.3109/10408418409105902, 1984. 550 Collins, S., Boyd,P. W. and Doblin,M. A.: Evolution,Microbes,and Changing Ocean Conditions,Ann. Rev. Mar. Sci.,12 (1) Table 3.

Tables
Key search terms (TS=(phytoplankton OR algae OR microalgae OR diatom OR coccolithophore* OR cyanobacteri* OR diazotroph*) AND TS=(stoichiometr* OR "chemical composition" OR "element* composition" OR "nutritional quality" OR "nutrient composition" OR "nutrient content" OR "nutrient ratio*" OR C:N OR C:P OR N:P OR P:C OR N:C OR "cellular stoichiometr*" OR C:N:P OR "element* ratio*" OR "food qualit*" OR "nutrient concentration" OR "carbon budget") AND TS = (phosph* OR "phosph* limit*" OR nitr* OR "nitr* limit*" OR iron OR "iron limit*" OR nutrient OR "nutrient limit*" OR "nutrient supply" OR "nutrient availabilit*" OR "supply ratio*" OR eutrophication OR fertili* OR enrichment OR temperature OR warming OR light OR irradiance OR "light limit*") AND TS = (marine or sea or ocean OR seawater OR aquatic)). Table 1. Key word 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. 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 = number of outlier studies). ; ? ( GGG , weighted mean stoichiometry sensitivity factor with environmental driver X and response variable Y; ln ( ) GGGGGGGGGG , weighted mean value of the natural logarithm-845 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 experiments with studentized residual value of ln(RR) exceeding 3 was removed as outliers. Red bold texts highlight statistically significant environmental driver using both effect sizes.  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 C:N:P of 148: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 855 Table 2. We used Equation (9) in the main text for estimating the combined effect of multiple drivers.