Interaction matters: Bottom‐up driver interdependencies alter the projected response of phytoplankton communities to climate change

Phytoplankton growth is controlled by multiple environmental drivers, which are all modified by climate change. While numerous experimental studies identify interactive effects between drivers, large‐scale ocean biogeochemistry models mostly account for growth responses to each driver separately and leave the results of these experimental multiple‐driver studies largely unused. Here, we amend phytoplankton growth functions in a biogeochemical model by dual‐driver interactions (CO2 and temperature, CO2 and light), based on data of a published meta‐analysis on multiple‐driver laboratory experiments. The effect of this parametrization on phytoplankton biomass and community composition is tested using present‐day and future high‐emission (SSP5‐8.5) climate forcing. While the projected decrease in future total global phytoplankton biomass in simulations with driver interactions is similar to that in control simulations without driver interactions (5%–6%), interactive driver effects are group‐specific. Globally, diatom biomass decreases more with interactive effects compared with the control simulation (−8.1% with interactions vs. no change without interactions). Small‐phytoplankton biomass, by contrast, decreases less with on‐going climate change when the model accounts for driver interactions (−5.0% vs. −9.0%). The response of global coccolithophore biomass to future climate conditions is even reversed when interactions are considered (+33.2% instead of −10.8%). Regionally, the largest difference in the future phytoplankton community composition between the simulations with and without driver interactions is detected in the Southern Ocean, where diatom biomass decreases (−7.5%) instead of increases (+14.5%), raising the share of small phytoplankton and coccolithophores of total phytoplankton biomass. Hence, interactive effects impact the phytoplankton community structure and related biogeochemical fluxes in a future ocean. Our approach is a first step to integrate the mechanistic understanding of interacting driver effects on phytoplankton growth gained by numerous laboratory experiments into a global ocean biogeochemistry model, aiming toward more realistic future projections of phytoplankton biomass and community composition.


| INTRODUC TI ON
Temperature, nutrient availability, light, and the concentration of inorganic carbon species are the main bottom-up drivers that control phytoplankton growth rates and, thus, marine primary production and phytoplankton community structure (Cooley et al., 2022).
Ongoing climate change modifies environmental conditions, leading to the emergence of new mean states of these drivers, with potential impacts on phytoplankton growth and ocean productivity (e.g., Winder & Sommer, 2012). Global marine biogeochemistry models including multiple phytoplankton functional types (PFTs) are used to project future marine primary production under a range of different climate change scenarios (e.g., Kwiatkowski et al., 2020;Laufkötter et al., 2015;Nakamura & Oka, 2019;Stock et al., 2020;Tagliabue et al., 2021). The majority of these models projects a decrease in future net primary production (NPP) of 2%-13% by the end of the century in a high-emission climate scenario (e.g., Bopp et al., 2013;Cooley et al., 2022;Kwiatkowski et al., 2019Kwiatkowski et al., , 2020Laufkötter et al., 2015;Steinacher et al., 2010;Tagliabue et al., 2021). This was explained by increasing nutrient limitation with increased upper ocean stratification and stronger grazing pressure under future conditions (e.g., Laufkötter et al., 2015;Nakamura & Oka, 2019).
Interactive effects of multiple co-occurring environmental changes on phytoplankton growth are, however, to a large extent not considered in these large-scale ocean models that usually apply PFT models (i.e., based on the trophic structure and taxonomic affiliation of organisms) or trait-based models (i.e., based on the ecological function and traits of organisms, Kwiatkowski et al., 2014) to simulate biological and biogeochemical processes (e.g., Kwiatkowski et al., 2020;Laufkötter et al., 2015).
For varying single bottom-up environmental drivers, phytoplankton growth responses in models are commonly described by mechanistic approaches such as the Arrhenius curve for temperature (Arrhenius, 1889), the Michaelis-Menten, Monod, or Droop functions for different nutrients (Droop, 1973;Michaelis & Menten, 1913;Monod, 1942Monod, , 1949, photosynthesis-irradiance curves for light (Geider et al., 1996;Geider & Osborne, 1992), and functions that describe growth responses to the carbonate system Gafar et al., 2018;Paul & Bach, 2020;Seifert et al., 2022). Subsequently, the applied growth descriptions in models usually consist of the multiplication of these single driver effects; for different nutrient limitations, often a minimum-based approach is chosen, justified by Liebig's law of the minimum. Physiologybased models like the model by Geider et al. (1998), which forms the base of the model employed here, take interactions between nutrient status and photoacclimation into account. In recent decades, the importance of so-called "interactive driver effects" became increasingly evident from laboratory experiments using full-factorial manipulations (e.g., Boyd et al., 2018;Folt et al., 1999;Gao et al., 2020;Häder & Gao, 2015;Rost et al., 2008), implying that the use of mechanistic functions for single drivers is not sufficient to describe growth responses of phytoplankton to jointly changing drivers. This means that a shift in one driver does not only affect the growth of a species directly but can also modulate its sensitivity to other drivers. For instance, ocean acidification decreases the tolerance of Southern Ocean diatoms to high light conditions (Trimborn, Thoms, et al., 2017) and reduces the temperature above which the growth rate of the coccolithophore Emiliania huxleyi becomes negative by 1-2°C (Listmann et al., 2016). Mutual modifications of driver responses may arise from energy re-allocations within the phytoplankton cells, that is, energetic and elemental costs are reduced at the new level of one driver liberating resources for other processes (Rokitta & Rost, 2012;Van de Waal & Litchman, 2020).
Thereby, such interactions can obviously shape the outcome of environmental changes, both in strength and direction of the response (Boyd et al., 2014;Vinebrooke et al., 2004), and may furthermore induce shifts in the phytoplankton community composition (Boyd & Hutchins, 2012;Vinebrooke et al., 2004).

Important attempts have been made to develop mechanistic
functions of bottom-up driver interactions for models. However, they either do not yet take into account the wealth of information collected in recent laboratory multiple-driver studies or are not thus far applicable to large-scale ocean models. For instance , Geider et al. (1996, Geider et al. ( , 1998) developed a light limitation term that depends on the nitrate limitation of a cell. This function is widely applied in biogeochemical models (e.g., Aumont et al., 2015;Hauck et al., 2013;Moore et al., 2004) but was developed from a relatively small set of experiments (most importantly, Laws & Bannister, 1980) and is limited to interactions between light and nitrogen only, although there exist extensions for silicon (Hohn, 2009) and iron (Buitenhuis & Geider, 2010). The chain model by Pahlow and Oschlies (2009), which takes into account interactions between phosphorus, nitrogen, and light limitations and which was applied in Arteaga et al. (2016) in a global model, starts from a different, optimalitybased approach. In another early study on parameterizing driver interactions, Smith (2011) developed temperature dependencies of Michaelis-Menten terms for nutrient limitation, which are, however, only based on a confined collection of experiments all performed in waters of the North Atlantic gyre over a 3-year period (Harrison et al., 1996). Similarly, Wrightson et al. (2022) developed temperature-dependent iron and phosphorus use efficiencies for diazotrophs in a global biogeochemistry model based on confined laboratory experiments by Jiang et al. (2018) and Yang et al. (2022). Thomas et al. (2017) describe a function for interactive effects of temperature and nutrient concentration on phytoplankton growth, which was developed based on conceptual understanding, and K E Y W O R D S biogeochemical modeling, bottom-up effects, coccolithophores, diatoms, interactive effects, multiple driver, ocean acidification, warming experiments were only used to assess the application of the function in a global species distribution model. Finally, Taherzadeh et al. (2019) describe a theoretical framework on interactive driver sensitivities of community growth rates and cell sizes to changing nitrate availability and grazing pressure for a zero-dimensional box model. This framework, however, cannot yet be accommodated in large-scale ocean models that use PFT and trait-based models without significant modifications of the typical model structure.
In summary, the cumulated knowledge from laboratory studies is not yet fully exploited to parameterize interactive driver effects in large-scale ocean biogeochemistry models.
In the present study, we develop a parametrization of interactive effects between the carbonate system and both temperature and light on phytoplankton growth for a global biogeochemical model by applying the findings from our recently published meta-analysis (Seifert et al., 2020). The meta-analysis assesses growth-rate measurements of published multiple-driver laboratory studies. While not sufficient data was available from other driver interactions to obtain robust results, increasing partial pressure of carbon dioxide (CO 2 ) was found to profusely dampen the growth-enhancing effects of high temperature and high light (Seifert et al., 2020). Here, we develop and implement a new parameterization of the dual driver interactions of CO 2 with temperature and light into the global ocean biogeochemistry model FESOM-REcoM (Hauck et al., 2013;Karakuş et al., 2021;Seifert et al., 2022). In the first step, both the initial model setup without driver interactions and the newly developed model setup with driver interactions are used to assess changes in biomass, NPP, and individual driver limitations between present-day and future conditions. In a second step, we compare our findings between the simulations with and without driver interactions at present-day and future conditions. This allows us to test the hypotheses that (i) projections for future global and regional marine phytoplankton biomass are altered by driver interactions and that (ii) phytoplankton groups are affected differentially, resulting in an altered community composition compared with control simulations without driver interactions.

| Ocean biogeochemical model
We use the global Regulated Ecosystem Model version 2 (REcoM-2-M) coupled to the Finite Element Sea Ice-Ocean Model (FESOM 1.4, Schourup-Kristensen et al., 2014;Wang et al., 2014). We briefly introduce REcoM and the relevant model equations here but refer the reader to the Supporting Information (Text S.1) for a more detailed description of the model. The structure of REcoM mainly follows Hauck et al. (2013), Karakuş et al. (2021), and Seifert et al. (2022) including two detritus groups (slow-sinking and fast-sinking particles), two zooplankton groups (small, fast-growing zooplankton and slow-growing polar macrozooplankton), and three phytoplankton groups (diatoms, coccolithophores, and small-sized phytoplankton).
Phytoplankton growth, that is, the biomass build-up of each phytoplankton group p, depends on a group-specific maximum growth rate (μ max , Table S.1) as well as sensitivities to nutrient availability (f N ), temperature (f T ), light (f L ), and the carbonate system (f CO 2 ):

| Parametrizations of dual-driver interactions
We base the new parametrization of dual-driver interactions on the meta-analysis by Seifert et al. (2020), which highlights that dualdriver interactions between the carbonate system (called "CO 2 " for simplicity) and temperature as well as CO 2 and light have statistically significant effects on phytoplankton growth rates. As both driver interactions (in the following called C × T and C × L, respectively) include the effect of CO 2 on growth, we add interactive effects with temperature and light to the model term describing carbonate system dependencies (f CO 2 , Bach et al., 2015;Seifert et al., 2022) which reads in the original formulation: The term f CO 2 depends on water column concentrations of bicarbonate (HCO − 3 ) and CO 2(aq) , as well as H + (derived from pH) and is limited to be ≤3 (i.e., can triple phytoplankton growth if no other driver limits growth) for numerical reasons. It accounts for growthenhancing effects of the substrates for photosynthesis (HCO − 3 and CO 2(aq) ) as well as growth-inhibiting effects by increasing stress on cell-internal homeostasis at decreasing pH. The parameters a, b, c, and d (with d = d initial in the original formulation) for each phytoplankton group were derived from least-square fits to laboratory data ( Table 1; following Seifert et al., 2022). Parameters a and b determine the response to substrate availability, c determines the minimum CO 2(aq) requirement, and d the response to decreasing pH levels (i.e., the effect of ocean acidification). For the parameterization of interactive effects, we only modify the part of the function, (2)  (Schoemann et al., 2005).
We here use the collection of growth rates from the analysed laboratory studies in the supporting information of Seifert et al. (2020) ( Table 2; Figure 1a). The studies were previously filtered for a sufficient number of replicates as well as controlled experimental conditions including stable carbonate system parameters (Seifert et al., 2020). In addition, cultures must have been fully acclimated to experimental conditions (minimum 1 week or seven generations; Table S.2). While the laboratory studies usually consider two drivers at two driver levels (Boyd et al., 2019;Kreyling et al., 2018), experiments accounting for gradients of driver levels are rarely done, mostly due to logistic reasons (Boyd et al., 2019;Collins et al., 2022).
Nonetheless, it is possible to infer the directionality and magnitude of the interaction between two drivers. The laboratory studies were conducted in a multifactorial design: Phytoplankton growth was measured at (i) low control levels of drivers A (either temperature or light) and B (CO 2 ; A1B1), (ii) high control levels of drivers A and B (A2B2), (iii) low control level of driver A, high control level of driver B (A1B2), and (iv) high control level of driver A, low control level of driver B (A2B1). The interaction of drivers is synergistic if the effect of simultaneous changes in multiple drivers on phytoplankton growth (A2B2) is larger than the sum of single-driver responses (i.e., A2B2 > A1B2 + A2B1 − A1B1). The interaction is antagonistic if drivers mutually dampen the growth response compared with the sum of individual driver responses (A2B2 < A1B2 + A2B1 − A1B1; Boyd & Hutchins, 2012;Folt et al., 1999;Przeslawski et al., 2015).
For each selected laboratory study, we compute a theoretical growth rate A2B2 add from the measured growth rates (A1B1, A1B2, A2B1; Table 2; Figure 1a), which represents growth at high driver levels under the assumption that drivers have no interactive effects, that is, are additive: To assess how much the interaction between drivers affects growth at high driver levels, we compute the relative difference g rel between the measured growth rate A2B2 and the hypothetical additive growth rate A2B2 add (Table 2; Figure 1a): The interaction is antagonistic for g rel < 0 and synergistic for g rel > 0. We then calculate the mean g rel for each driver interaction and phytoplankton group (Table 2; Figure 1b). The values of g rel yield synergistic effects of C × T for diatoms (g rel = 1 % ± 13 %) and coccolithophores (g rel = 11 % ± 30 %), and antagonistic effects of C × T for small phytoplankton (g rel = − 6 % ± 3 %). Interactive effects of C × L are antagonistic for diatoms (g rel = − 11 % ± 15 %) and coccolithophores (g rel = − 2 % ± 34 %), and synergistic for small phytoplankton (g rel = 12 % ± 27 %; Table 2). We acknowledge that the spread is large across studies (standard deviation of up to ±34% for the coccolithophore C × L interaction). We trace this variability back to the diversity of species, origins of isolates, and driver ranges within the collection of studies (Table S.2). In addition, also the propagated standard deviations of A2B2 add in the individual studies (i.e., the square-rooted sums of the squared standard deviations of A1B1, A1B2, and A2B1) are rather high (up to ±0.43 day −1 ; Table S.2). Therefore, we performed additional sensitivity simulations to assess the robustness of our parameterization (see Section 2.3 and Text S.2).
For the subsequent development of model terms for driver interactions, we also calculate the mean temperature, light, and CO 2 values at low and high control levels from all laboratory studies. This calculation yields 380 ± 17 atm, 14 ± 3°C, and 52 ± 24 mol photons m −2 s −1 for the low CO 2 (LC), temperature (LT), and light (LL) levels, respectively, and 880 ± 125 atm, 28 ± 2°C, and 295 ± 175 mol photons m −2 s −1 for the high CO 2 (HC), temperature (HT), and light (HL) levels, respectively. Laboratory light levels are converted to levels of photosynthetically active radiation (PAR), that is, from mol photons m −2 s −1 to W m −2 following Brock (1981) for a wave length distribution corresponding to PAR at the sea surface, resulting in 11.3 W m −2 (LL) and 64.1 W m −2 (HL).

| Interaction between CO 2 and temperature
In the following, we describe how we modify the growth rates GR p resulting from the original (initial) model Equation (1) at high control levels of CO 2 and temperature (HCHT), further on called GR HCHT initial , according to g rel (hence, the relative change caused by driver interactions). For that, we use a temperature-dependent response to ocean acidification, hence, a term d inter,C×T that substitutes d in Equation (2) by considering the interaction (inter) between CO 2 and temperature (C × T; Figure 1c; Table 2). The growth rates that result from this modified model equation and that take interactive effects between CO 2 and temperature into account (4) g rel = A2B2 − A2B2 add ∕ A2B2 add . TA B L E 2 Growth rates (day −1 ) resulting from CO 2 and temperature as well as CO 2 and light interactions based on experimental data in the literature collection by Seifert et al. (2020). Note: Driver A is either temperature or light, driver B is always CO 2 ; 1 and 2 denote low and high driver levels, respectively. In the column A2B2 add theoretical growth rates at high driver levels assuming additive driver interactions are listed (Equation 3), and g rel denotes the relative difference between A2B2 add and the realized growth rates of the experiments A2B2 in % (Equation 4). Positive values of g rel indicate synergistic, and negative values of g rel indicate antagonistic interactions. Mean relative changes for each interaction and phytoplankton group are printed in bold (standard deviation in parentheses). For additional details on the individual studies, see Table S.2. are called GR HCHT inter . That means, if g rel = 0.1, d inter,C×T should increase GR HCHT initial by 10% yielding GR HCHT inter , which would represent a synergistic interaction of C × T. Note that even though not noted in the equations for simplicity, all calculations were performed for each phytoplankton group separately.
We first calculate the growth rates at LC and LT levels (i.e., 380 μatm and 14°C), GR LCLT . To derive concentrations of HCO − 3 , CO 2(aq) , and H + to use in f CO 2 that comply with the CO 2 levels ap-  Table 2). Likewise, we calculate growth rates at high CO 2 and temperature control levels (i.e., 880 μatm and 28°C; GR HCHT initial ): GR HCHT initial (derived from the model equations) corresponds to the values of A2B2 add (computed as described above, Table 2), and, thus, represents growth rates at high control levels without driver interaction.
We then assume a relative change in GR HCHT initial according to g rel ( Table 2), which is caused by the interaction of C × T. The resulting interactive growth rate, GR HCHT inter , is computed from GR HCHT initial by: In the next step, we compute d HCHT for each phytoplankton group, which yields values of d (Equation 2) that result in the interactive growth rates GR HCHT inter by solving the following equation for d HCHT : We assume that for each phytoplankton group, values of d initial (Table 1) are valid at our low control levels (i.e., 380 μatm and 14°C), and values of d HCHT are valid at our high control values (i.e., 880 μatm and 28°C). In the absence of information about the functional form, we additionally assume that d changes linearly and temperaturedependent between d initial and d HCHT (Figure 2). The slope m inter,C×T of that linear change of d is computed as: We use the resulting values of m inter,C×T to derive a temperaturedependent formulation for d, which is d inter,C×T : (9) m inter,C×T = d HCHT − d initial HT − LT .

F I G U R E 1
Stepwise description of the parameterization of interactive effects between CO 2 and temperature as well as light on phytoplankton growth. Abbreviations are given at the right side of the figure, and more details on parameters and calculations are given in the text. Note the different annotations for real growth rates that are based on the literature compilation of laboratory data ("Laboratory data") and computed growth rates based on the model equations of phytoplankton growth ("Model equations").
with T deg being temperature in °C (Figures 2, 3). Values of d HCHT and m inter,C×T are listed in Table 3, and slopes are illustrated in Figure S.1.
We assume that the linear change of d inter,C×T and thus, the strength of interaction between CO 2 and temperature, changes linearly with the two variables also outside the calibration interval ( Figure 2). This means that the mode of interaction changes below 14°C (for coccolithophores, for instance, from synergistic to antagonistic interaction), and the magnitude of the interactive effect gets closer to additive interaction below 380 μatm. Above 28°C and 880 μatm, the magnitude of interactive effects becomes larger.
We acknowledge that a linear interpolation of interactive effects between low and high control levels and an extrapolation to even lower/higher driver levels is a simplification of interactions across possible driver ranges, which can theoretically be described as a 3-dimensional landscape (Boyd et al., 2018;Collins et al., 2022;Cooley et al., 2022). However, response patterns over multiple levels of two continuous drivers, which would be required to develop such landscape functions (Kreyling et al., 2018), have not yet been determined in laboratory studies. Further, fitting complex regressions to the small number of driver levels that are available (usually only two driver levels are investigated in multiple-driver laboratory studies) would possibly add bias through overfitting (Babyak, 2004). We, therefore, consider a linear regression as the most simple approach until enough data are available for three-dimensional landscape functions.

| Interaction between CO 2 and light
The effects of C × L interactions are computed analogously to the interaction with temperature ( Figure 1d). Here, we use the low and high control levels of light (11.3 and 64.1 W m −2 , respectively) and the light function f L (Equation S.5) to compute growth at low and high CO 2 and light control levels (i.e., GR LCLL and GR HCHL initial following Equations 5 and 6). Furthermore, we use the relative growth rate changes g rel caused by C × L interaction as listed in Table 2 to compute GR HCHL inter (following Equation 7). For the light function f L in the present computation, we assume fixed, but physiologically meaningful values for the variable chlorophyll-to-carbon ratio q Chl in each group (Nielsen, 1997;Sathyendranath et al., 2009;Thomalla et al., 2017) derived from earlier model simulations [0.01 g Chl a (g C) −1 for diatoms and 0.005 g Chl a (g C) −1 for coccolithophores and small phytoplankton]. From GR HCHL inter we compute d HCHL and m inter,C×T following Equations (8) and (9), and derive d inter,C×L : with PAR in W m −2 (Figure 3b,e). Altogether, this results in a set of values for d HCHT , d HCHL , and m inter for each phytoplankton group (diatoms, coccolithophores, small phytoplankton) and interaction (C × T and C × L, respectively), which are listed in Table 3.

| Model simulations
We perform two sets of climatological simulations, one with presentday and one with future repeated year atmospheric forcing including CO 2 to assess the possible range of interactive effects between two climate states ( Table 4). All runs are performed on the mesh that was also used in FESOM for the CORE-II model intercompari- F I G U R E 2 Conceptual representation of the relation between driver levels (LCLT, HCHT = low and high CO 2 and temperature control levels, respectively), d inter,C×T including the slope of the linear relationship m inter,C×T , and the relative change in the growth rate (GR). For instance, to derive a growth rate GR HCHT inter at high control levels of both drivers that is 11% higher than the respective non-interactive growth rate GR HCHT initial (as given by g rel in Table 2), d HCHT (Table 3) is lower than d initial ( Table 1). The driver space between d initial and d HCHT , d inter,C×T , is assumed to be linearly dependent on temperature (Equation 10). For a detailed overview of temperature-and light dependent changes in d inter of the separate phytoplankton groups, see While it was criticized for its underlying assumption on maybe unrealistic high emissions (Hausfather & Peters, 2020), the resulting climate state may nevertheless be reached even with lower emissions due to carbon cycle feedbacks that are currently not accounted for (Schwalm et al., 2020). Acknowledging this uncertainty, we use it to test maximal future effects of driver interactions on phytoplankton and biogeochemistry. To reduce the impact of climate variability, we identify "normal" atmospheric conditions by jointly assessing the state of the Southern Annual Mode (SAM index), the El-Niño Southern Oscillation (Oceanic Niño index), and the North Atlantic Oscillation (NAO index) for the present day (2000-2020) and the future (2080-2100) period. In particular, we compute seasonal anomalies of all indices over the respective periods and rank the years according to the annual means of the anomalies and their standard deviations, which yields 2014 and 2088 as years with the smallest F I G U R E 3 (a-c) Coccolithophore growth rates (arbitrary scale) for an example CO 2 system at high temperature and light control value with (solid lines) and without driver interaction (original parametrization, dashed lines). Arrows at 880 μatm mark the growth rate increase (+11%) for the CO 2 and temperature interaction, and the growth rate decrease (−2%) for the CO 2 and light interaction compared with the original parametrization. The growth rate resulting from the full interaction [(c), +9%] is caused by the individual interactions between CO 2 and temperature (a) and between CO 2 and light (b). (d-f) Display the respective f CO 2 functions at increasing pCO 2 levels. In addition to the interactive (solid lines) and the original growth rate formulation (dashed lines) at high driver control levels, we display f CO 2 for low temperature and light levels (dotted lines). At high pCO 2 with d inter and with temperature and light values below 14°C and 11.3 W m −2 (2°C and 5 W m −2 in the example) the direction of change in the slope of f CO 2 reverses. Abbreviations: C × L: CO 2 and light interaction; C × T: CO 2 and temperature interaction; HC, high CO 2 ; HL, high light; HT, high temperature.

TA B L E 3
Parameter values for d HCHT and d HCHL , respectively (in kg mol −1 ) and m inter (in °C mol kg −1 and W mol m −2 kg −1 , respectively) of each phytoplankton group and driver interaction.

TA B L E 4 List of model simulations performed in this study.
Forcing Interaction None Abbreviations: C × L = CO 2 and light interaction; C × T = CO 2 and temperature interaction; C × T × L = both CO 2 and temperature and CO 2 and light interaction.

C × T C × L C × T × L
anomalies for the present-day and the future period, respectively.
We use annual mean atmospheric CO 2 levels of 400 ppm for 2014 and of 971 ppm for 2088 corresponding to the global atmospheric CO 2 levels of the respective years used in the AWI-CM simulations.
The difference between the present-day and the future climate under the chosen scenario outweighs any interannual variability, and the choice of any particular year within the present-day and future periods, therefore, has a negligible effect on the analysis performed here. Temperature and salinity are initialized with hydrographic data from the Polar Science Center Hydrographic Climatology (Steele et al., 2001). The nutrients dissolved inorganic nitrogen and dissolved silicic acid are initialized from the World Ocean Atlas (Garcia et al., 2013), and dissolved iron from PISCES output (Aumont et al., 2003), which was corrected using observed profiles (Boye et al., 2001;de Baar et al., 1999). In both the present-day (PR) and the future (FU) set, we perform four simulations (Table 4). In a control simulation, we use the original CO 2 function with d = d initial (Table 1) in Equation (2)  Phytoplankton biomass is, hence, relatively robust against variations of g rel .

| Model analysis
We evaluate total biomass, NPP, limitation terms, and biogeochemical fluxes for each simulation globally and in three larger regions, namely the northern and southern high latitudes as well as the tropics/subtropics (Fay & McKinley, 2014) Total biomass and NPP were integrated vertically over the whole water column and horizontally over the respective region. Limitation terms (f N , f T , f L , and f CO 2 ) were calculated online and averaged over the depth layer from surface to the euphotic zone depth, that is, where PAR is 1% of its surface value ( Figure S.3). We note that thereby, we mix different ecological niches with diverging environmental conditions (e.g., low nutrient concentrations and high light in the mixed layer, high nutrient concentrations and low light below the mixed layer), but we consider limitation terms in the euphotic zone to be most representa-

| Projections of end-of-century biomass and attribution to drivers in the control simulations
Global total biomass and NPP decrease by 68.0 Tg C and 1.6 Pg C year −1 between the FU_CTRL and the PR_CTRL simulation Globally, decreasing coccolithophore and small-phytoplankton biomass (−10.8% and −9.0%, respectively; Figure 5) are primarily driven by stronger growth-limitation by CO 2 and light (Table 6).
Stronger light limitation is caused by a deepening of the euphotic zone in most parts of the ocean (Figure S.7). Light limitation terms are additionally modified by alterations in nutrient availability, temperature, and the chlorophyll-to-carbon ratio (Equation S.5).
Warming-induced growth enhancement (up to 9.9%) dampens the negative effect of changes in the CO 2 and light limitation terms for small phytoplankton and coccolithophores, and outbalances completely the growth-decreasing effects of these drivers on diatoms ( Figure 5; Table 6).
In the northern high latitudes of the control simulations, enhanced growth limitations by CO 2 and light (up to −7.0%; Table S.3) result in lower small-phytoplankton biomass compared with present-day conditions (−9.1%, Figure 5). Diatom biomass decreases

F I G U R E 4 Depth-integrated phytoplankton biomass at present-day (PR; a-d) and future (FU; e-h) conditions in the interaction simulations (PR_INTER and FU_INTER). Biomass in the FU_INTER simulation is shown as difference to the PR_INTER simulation.
F I G U R E 5 Relative differences of spatially and depth-integrated phytoplankton biomass between the control simulations at future (FU) and present-day (PR) conditions (FU_CTRL and PR_CTRL; hatched bars) as well as between the interaction simulations at future and present-day conditions (FU_INTER and PR_INTER; filled bars) globally, in the northern and southern high latitudes, as well as in the tropics/ subtropics. Regions as indicated on the map. Total biomass differences are denoted at the respective bars.
by −0.8% because of stronger growth limitation by CO 2 and nutrients (up to −3.4%; Table S.3). The decrease in biomass of both groups is reinforced by a warming-induced stronger grazing pressure (up to +34%). Coccolithophores benefit from a strong warming-related growth enhancement under future compared with present-day conditions (+27.6% biomass increase).
In the tropics/subtropics, where the deepening of the euphotic zone is most pronounced (up to 300 m; Figure  In the southern high latitudes, small-phytoplankton biomass decreases from present-day to future conditions (−13.1%; Figure 5) due to stronger limitations by CO 2 , nutrients, and light (up to −8.5%) as well as increasing grazing pressure (+6.8%; Table S.3). Diatoms and coccolithophores benefit (+14.5% and +49.7%, respectively) from the growth-increasing changes in the temperature limitation terms (up to 24%; Table S.3). In addition, the grazing pressure on diatoms becomes smaller (−41.7%).

| Effects of driver interactions on phytoplankton biomass in the present-day simulations
The consideration of driver interactions affects total present-day biomass only little in all regions, with relative changes between the PR_CTRL and the PR_INTER simulation of total global phytoplankton biomass being in the range of ±2.5% (Figure 6a). However, biomass differences between both simulations are larger on regional scales and within phytoplankton functional groups. Regionally, biomass differs most in the southern high latitudes, with lower diatom biomass (−20.9% or 32.3 Tg C) and higher small-phytoplankton biomass (+10.6% or 23.1 Tg C) in the PR_INTER compared with the PR_CTRL simulation (Figure 6a). High relative differences between the PR_CTRL and the PR_INTER simulation in coccolithophore biomass of +35.9% in the northern high latitudes and −10.5% in the southern high latitudes translate into small total biomass differences (1.3 and 0.4 Tg C, respectively) due to the generally low coccolithophore biomass in these regions (Figures 4b and 6a). Limitation terms in the PR_INTER simulation are in a similar range as in the PR_CTRL simulation (Tables S.3 and S.4).
In the tropics/subtropics and partly in the subpolar regions, interactive effects between CO 2 and light or temperature change phytoplankton growth rates by more than 0.001 day −1 between the PR_INTER and the PR_CTRL simulation (which corresponds to a relative growth rate change of more than 1-2%), while driver interactions barely alter the growth rates in the polar regions compared with the control simulation (difference in growth rates <0.001 day −1 ; compare white and colored regions in Figure 7a-c). Diatoms are predominantly affected by the C × L interaction (Figure 7a), which lowers their growth rates compared with the CO 2 term without the driver interactions ( Table 2). Regions with high coccolithophore biomass ( Figure 4b) are mainly influenced by the growth-enhancing C × T interaction ( Figure 7b; Table 2). For small phytoplankton, the growthstimulating C × L interaction dominates in most parts of the ocean ( Figure 7c; Table 2). In subpolar regions, the growth-dampening C × T interaction is more prevalent but likely outbalanced by the C × L interaction indicated by a higher small-phytoplankton biomass in the PR_INTER compared with the PR_CTRL simulation ( Figure 7c;

| Effects of driver interactions in the future simulations
While the difference in total global phytoplankton biomass between the FU_CTRL and the FU_INTER simulation remains small (1.5%, 19.1 Tg C), biomass differences of the individual phytoplankton groups are three to four times more pronounced than in the presentday simulations when accounting for driver interactions, with globally 15.4% fewer diatoms, 65.2% more coccolithophores, and 9.5% more small phytoplankton (Figure 6b). Regionally, the southern high latitudes exhibit the largest relative and total differences of diatom and small-phytoplankton biomass between the FU_CTRL and the FU_INTER simulation, respectively, with 36.1% lower diatom biomass (63.9 Tg C) and 26.3% higher small-phytoplankton biomass (49.7 Tg C; Figure 6b). Here, coccolithophore biomass barely changes between both simulations, and total phytoplankton biomass is 3.8% lower in the interaction simulation. In the northern high latitudes, driver interactions have the strongest relative effect on total phytoplankton biomass (+5.3% or 12.6 Tg C compared with the FU_CTRL simulation; Figure 6b). Here, diatom biomass is lowered by 10.5% TA B L E 6 Relative differences of globally averaged limitation terms and growth rates and total differences in global depth-integrated biomass between the present-day and future simulations (control simulations: PR_CTRL and FU_ CTRL, interaction simulations: PR_INTER and FU_INTER) as well as total biomass and relative biomass contributions of the three phytoplankton groups within each simulation.
and small-phytoplankton biomass is enhanced by 18.7%. Like in the present-day simulations, large relative interactive effects on coccolithophore biomass in the polar regions translate into only small total biomass differences. In the tropics/subtropics, higher phytoplankton biomass in the FU_INTER compared with the FU_CTRL simulation (diatoms: +5.6%, coccolithophores: +60.9%, small phytoplankton: +0.8%) results in a higher total phytoplankton biomass of 3.2% (20.3 Tg C; Figure 6b).
Generally, the CO 2 term (Equation 2), which includes the parametrization of driver interactions in the FU_INTER simulation and which F I G U R E 6 Total differences of spatially and depth-integrated phytoplankton biomass (a) between the interaction (INTER) and the control (CTRL) simulation at present-day conditions (PR_INTER and PR_CTRL), and (b) between the interaction and the control simulation at future conditions (FU_INTER and FU_CTRL) globally, in the northern and southern high latitudes, as well as the tropics/subtropics. Regions as indicated on the maps. The bars reveal the effect of the interactive driver parametrization at present-day and future conditions, respectively. Percentages above/below each bar specify the relative biomass changes. (c) Relative changes in limitation terms (CO 2 , nutrients, light, temperature) and growth rates between the FU_INTER and the FU_CTRL simulation for the same regions. Limitation terms and growth rates were averaged over the euphotic zone (1% PAR to +8% between both simulations but is not always the dominant driver that causes differences in phytoplankton biomass (Figure 6c). Globally, differences in diatom and small-phytoplankton biomass between the two future simulations (Figure 6b) are mainly caused by the light limitation term (−6.5% and +6.4%, respectively, in the FU_INTER compared with the FU_CTRL simulation; Figure 6c).

Instead, variations in the CO
Coccolithophores benefit from a more growth-enhancing CO 2 term (+3.0%; Figure 6c), which is likely the result of the interaction between increasing CO 2 and warming ( Figure 7e; Table 2).
In the northern and southern high latitudes, biomass differences of diatom and small phytoplankton between the future simulations are mainly caused by the light limitation term (up to −7.8% and +9.6% in the FU_INTER compared with the FU_CTRL simulation; Figure 6b,c), similar to the global ocean. Small-phytoplankton growth is additionally enhanced by changes in the CO 2 term (up to +2.6%; Figure 6c). In contrast to the global ocean, the CO 2 term of coccolithophores is more growth-limiting in the FU_INTER compared with the FU_CTRL simulation (up to −5.9%; Figure 6c). This dampens the growth-enhancing effect of alleviated nutrient limitation (+0.9%;

| Altered phytoplankton response to climate change with driver interactions
The decrease in global total biomass and NPP from present to future conditions is about the same in the simulation with driver interactions (−79.3 Tg C and −2.3 Pg C year −1 ; FU_INTER minus PR_INTER) F I G U R E 7 Spatial distribution of the dominating interaction for present-day and future climate in the control and the interaction simulations. For the present-day (PR) simulations (a-c), mesh nodes with a growth rate difference of ≥0.001 day −1 between the control simulation (PR_CTRL) and the simulation with CO 2 and temperature interaction (PR_CT) or CO 2 and light interaction (PR_CL) were extracted and categorized according to the larger impact on growth rate changes into CO 2 and temperature interaction (C × T) and CO 2 and light interaction (C × L). Subsequently, these categories at each node were summarized over depth by weighting according to the size of the depth level and by selecting the most common interaction. No interaction dominates when the differences in growth rates between the PR_CTRL and the PR_CT or PR_CL simulations are <0.001 day −1 . A similar computation was performed for the future (FU) simulations (d-f; i.e., by using the FU_CTRL, FU_CT, and FU_CL simulations). and the one without (−68.0 Tg C and −1.6 Pg C year −1 ; FU_CTRL minus PR_CTRL), that is, about 5%-6% ( Figure 5; Figure S.4; Table 5), but we identify differences on a regional and group-specific level.
From present to future conditions, global diatom biomass decreases (−8.1%) in the interaction simulation while remaining constant in the control simulation (Figures 4e and 8; Figure S.8; Table 6). This can mainly be attributed to a stronger decrease of the light term (−7.9% in the interaction simulations vs. −3.4% in the control simulations; Table 6). Small-phytoplankton biomass decreases less (−5.0% instead of −9.0%) from present to future conditions in the interaction simulation (Figures 4g and 8; Figure S.8; Table 6) because light limitation intensifies less severely than in the control simulation (−5.5% instead of −8.4%; Table 6). The biomass of coccolithophores even increases (+33.2%) instead of decreases (−10.8%) once driver interactions are considered ( Figure 4f; Table 6), mainly because the CO 2 term dampens future growth less without driver interactions (−6.3% instead of −8.4%; Table 6). A reduced grazing pressure on coccolithophores (−7.1%) despite increasing coccolithophore biomass, likely because of the higher availability of small phytoplankton as preferred prey, contributes to a stronger biomass build-up (  Table 6).
Similar to the global ocean, the phytoplankton community in the northern high latitudes and in the southern high latitudes is composed of a higher share of coccolithophores and small phyto- F I G U R E 8 Schematic figure summarizing phytoplankton biomass changes from present-day (PR) to future (FU) climate (regions as indicated on the maps). Dashed arrows represent the biomass changes between the control simulations (PR_CTRL and FU_CTRL; upwarddirected = higher biomass, downward-directed = lower biomass) for each phytoplankton group separately and together, without specifying the magnitude of change. For each group, the dominant drivers responsible for the biomass change are indicated on the arrows. Solid arrows represent the biomass changes between the interaction simulations (PR_INTER and FU_INTER) and are scaled according to the magnitude of change compared with the control simulations. For instance, diatom biomass in the northern high latitudes decreases stronger from the PR_ INTER to the FU_INTER simulation than from the PR_CTRL to the FU_CTRL simulation, while the decrease in small-phytoplankton biomass is weaker. In case both arrows have the same size, changes were not affected by accounting for interactive effects. The delta signs at each arrow indicate the main drivers responsible for the difference to the control simulations (only indicated when the biomass change differs between the control and the interaction simulations). , opal export at 100 m depth increases by +4% in the control simulations from present-day to future but decreases by −13% in the interaction simulations (Table 5). Calcium carbonate export at 100 m depth increases (+24%) instead of decreases (−20%) between present-day and future conditions once driver interactions are considered (Table 5). Similar to the global ocean, POC export at 100 m depth of the southern high latitudes decreases by a similar amount from present-day to future conditions in the control and the interaction simulation (by 0.12 and 0.10 Pg C year −1 , respectively; Table 5).
Alongside with a decreasing instead of increasing diatom biomass in the interaction simulation ( Figure 5), opal export at 100 m becomes lower (−20%) and not higher (+12%) under future conditions ( Table 5). The stronger increase in coccolithophore biomass in the interaction compared with the control simulation ( Figure 5) results in a stronger increase of the calcium carbonate export at 100 m depth (+90% instead of +27%; Table 5).

| DISCUSS ION
Our study introduces a parametrization that accounts for interacting effects between CO 2 , temperature, and light on phytoplankton growth in a global ocean biogeochemistry model. We show that the overall decrease of global phytoplankton biomass and productivity under future conditions remains unchanged when accounting for driver interactions, whereas the projected changes on a regional scale are altered, with implications for biogeochemical fluxes.

| Drivers of changes in future phytoplankton biomass without interaction
Under future climate conditions assuming persistently high emissions, our simulations with and without driver interactions reveal 5%-6% lower phytoplankton biomass and 6%-8% lower NPP compared with present-day levels. This is in the range of previous modeling studies that apply similar scenarios [e.g., −2% to −13% NPP in the models analysed by Steinacher et al. (2010), −8.6% NPP in the models analysed by Bopp et al. (2013), −6.1% biomass in the models analysed by Kwiatkowski et al. (2019), −3.0% NPP in the models analysed by Cooley et al. (2022)]. The patterns of change are consistent with previous findings, with declines in NPP and biomass mainly in the temperate, subtropical, and tropical regions, and increases in parts of the Southern Ocean, the North Atlantic, and the equatorial Pacific (Figure 4h; Bopp et al., 2013;Cabré et al., 2015;Cooley et al., 2022;Lotze et al., 2019;Moore et al., 2018;Vichi et al., 2011).
In our study, decreasing growth rates and shifts in the phytoplankton community under future conditions are not driven by one  Laufkötter et al. (2015), causing unchanged NPP under future conditions in some of the analysed models (Laufkötter et al., 2015).
Yet, the identification of dominant drivers critically depends on the depth layer considered for the computation of driver strengths. In the low latitudes, light limitation barely changes from present to future conditions in the study of Laufkötter et al. (2015), who investigate surface limitation terms, and nutrient limitation becomes significantly stronger in the study of Marinov et al. (2010), who compute limitation terms over the ML depth, which is projected to become shallower in the future tropics and subtropics. Limitation terms in our study were computed over the euphotic zone, which deepens under future conditions and, hence, decreases light availability but increases nutrient availability. Moreover, the growthincreasing effect of warming in Laufkötter et al. (2015) may be an upper limit as the increase in temperature at the sea surface is higher than at depth (Fu et al., 2016). Our study shows, however, that warming within the euphotic zone can in most cases not fully compensate, but dampen the growth-decreasing effects of other drivers. Finally, growth responses to future environmental conditions in our model are additionally modified by changes in the carbonate system (Seifert et al., 2022), which can even dominate the response for some regions and phytoplankton groups and is not considered in other models. While the decreasing phytoplankton biomass in most models analysed by Laufkötter et al. (2015) was mainly caused by warming-induced increasing grazing pressure on a global scale, this effect appears more important on a regional and group-specific scale in our model (small phytoplankton in the northern high latitudes and only diatoms globally). Top-down effects will be discussed in more detail in Section 4.3. In summary, the identification of predominant drivers responsible for changes in phytoplankton biomass depends on which depth layer is considered.

| Future changes in phytoplankton biomass and biogeochemical fluxes
Based on our model projections, coccolithophores and small phytoplankton are 'winners' and diatoms are 'losers' in the simulations that include driver interactions on a global scale (Figure 8).
In the northern high latitudes, the increase of coccolithophore biomass from present-day to future conditions in the simulations without driver interactions (Figure 8; Figure S.8) is in line with observed recent coccolithophore biomass increase in the North Atlantic (Beaugrand et al., 2013;Rivero-Calle et al., 2015).
We show that this trend intensifies with dual-driver interactions ( Figure 8; Figure S.8). Coccolithophore biomass also increased in the future tropics/subtropics when interactive effects were considered, while it decreased in the future tropics/subtropics of the control simulation ( Figure 8). Laboratory experiments show that coccolithophores are often negatively affected by ocean acidification (Hoppe et al., 2011;Meyer & Riebesell, 2015;Seifert et al., 2020). While this is mirrored in the growth-decreasing CO 2 term in the model equations and in the result of our control simulation, interactions with other environmental drivers as well as physiological and ecological feedback can even reverse the sign of the biomass change (Seifert et al., 2022).
Notably, in the southern high latitudes, diatom biomass is projected to increase in the control simulations but to decrease in the interaction simulations. Previously, future Southern Ocean diatom biomass and NPP were projected to increase in the PlankTOM5 Biogeochemical Model (Laufkötter et al., 2015) and the Biogeochemical Elemental Cycle model (BEC, Krumhardt et al., 2022) because of warming and enhanced nutrient supply due to increasing upwelling rates and nutrient trapping (Krumhardt et al., 2022;Laufkötter et al., 2015;Moore et al., 2018). In our control simulation, diatom biomass in the future Southern Ocean increases because of warming, improved light availability, and reduced grazing pressure, but driver interactions counteract these biomass-enhancing effects.
In addition to warming, reviews of modeling, experimental, and observational results identified increasing iron supply as a likely cause for increasing diatom biomass in the Southern Ocean (Deppeler & Davidson, 2017;Henley et al., 2020). Further, experiments have shown that ocean acidification dampens the positive effect of iron enrichment on Southern Ocean diatoms . Thus, including interactions involving the nutrient limitation term may yield a more complete picture of the future of Southern Ocean diatoms.
Our results show that the future POC flux is reduced by a similar amount in both the control and the interaction simulation compared with present-day fluxes (Table 5). However, while a community shift from large to small phytoplankton was detected in Bopp et al. (2005) without driver interaction, this shift is only seen in the simulation with driver interaction. As diatoms form silicious cell walls (frustules) and coccolithophores build calcium carbonate structures around the cells (coccospheres), the altered projection of the future phytoplankton community composition affects the export of these minerals to depth. More specifically, the export of opal is reduced (−13.3%) instead of increased (+4.1%) and the export of calcium carbonate increases (+23.5%) instead of decreases (−20.0%, Table 5) at 100 m depth from present-day to future conditions. In the southern high latitudes where the consideration of driver interactions has the largest effects on phytoplankton community composition, the opal export is decreasing (−20%) instead of increasing (+12%), and the calcium carbonate export increases more strongly (+90% instead of +27%, Table 5). The contribution of these minerals to the total sinking flux at >1000 m depth was observed to correlate with the POC flux, leading to the assumption that minerals ballast sinking particles and thereby increase the transfer of carbon to the deep ocean (Francois et al., 2002;Klaas & Archer, 2002). To date, our model does not account for the ballasting of minerals, and we can only speculate how these shifts in mineral export can affect carbon transfer to deeper water layers. While the ballasting effect of frustules is still a matter of discussion and presumably depends on the diatom species composition and the degree of silicification (Tréguer et al., 2018), there is consensus about the ballasting effect of coccospheres (Bach et al., 2016;Francois et al., 2002). In addition, zooplankton fecal pellets were found to sink faster when they contain frustules or coccospheres (Iversen & Ploug, 2010;Ploug et al., 2008). Because of the smaller share of coccolithophores compared with diatoms, we assume that an enhanced calcium carbonate flux cannot balance the reduced ballasting effect of frustules both on global scale and in the southern high latitudes, which may lead to a further reduction of the future POC flux. This could be amplified by a reduction of a cell's silicification or calcification caused by acidification (Findlay et al., 2011;Petrou et al., 2019;Raven & Beardall, 2021), and by a diet-related shift toward a microzooplankton-dominated grazer community producing fecal pellets with a lower sinking velocity (e.g., Steinberg & Landry, 2017;Wilson et al., 2008), which is due to the fact that small phytoplankton cells are the usual diet of microzooplankton, while copepodes and heterotrophic dinoflagellates usually graze on large diatom cells (e.g., Calbet, 2008;Landry & Calbet, 2004;Sommer et al., 2005;Sommer & Sommer, 2006).
A reduction in the dissolution of sinking opal by ocean acidification (Taucher et al., 2022) could furthermore strengthen silicic acid deficiencies in surface layers and lead to an ever stronger decrease in diatom biomass and silicification, especially in regions where silicic acid is limiting. The reduced usage of silicic acid in the southern high latitudes due to a lower diatom biomass, however, could reduce the role of the Southern Ocean as a "silicon trap" (Hauck et al., 2018;Holzer et al., 2014;Nissen et al., 2021) as less silicic acid is exported and trapped in the Southern Ocean's interior. Consequently, more silicic acid could be available for diatoms outside the Southern Ocean (Hauck et al., 2018). Finally, other modeling studies show that not only the transfer of carbon (Gehlen et al., 2006;Lima et al., 2014), but also the transfer of nutrients (Oka et al., 2008) is largely determined by the mineral load of sinking particles globally, which is directly linked to the plankton community structure (Lima et al., 2014;Nissen et al., 2021).

| Do interactions matter for projections of future phytoplankton?
Our results show that interactive effects between bottom-up drivers can alter the magnitude and even the direction of future phytoplankton biomass projections and are, hence, important amendments to the traditional temperature-nutrient-light-driven system in PFT models. The quantitative effects of modified bottom-up drivers can yet be obscured by the grazing function applied and by the number of grazer groups considered in the model (e.g., Karakuş et al., 2022;Le Quéré et al., 2016;Prowe et al., 2012;Vallina et al., 2014). For instance, solely the separation of the generic small zooplankton group into micro-and mesozooplankton in a different version of our model lead to a 25% increase in NPP (Karakuş et al., 2022). Furthermore, driver interactions can also affect the response of zooplankton and, hence, grazing to changing environmental conditions, although the database is still very limited, for example, due to the longer generation time of zooplankton (Otto et al., 2020). For example, the acclimation of the copepod Acartia tonsa to warming was shown to be impaired by ocean acidification (Dam et al., 2021), but responses to concurrent warming and acidification are highly variable depending on the trait measured in the experiment (e.g., growth, size, mortality; Garzke et al., 2020). Nonetheless, despite the regulating effect of grazing on phytoplankton biomass, the regional and seasonal effect of altering environmental conditions on phytoplankton is controlled by a close interplay between bottom-up and top-down effects (Hashioka & Yamanaka, 2007;Ward et al., 2014) and an improved consideration of bottom-up processes in a model is a vital step toward more realistic future projections of phytoplankton biomass.

| Limitations and caveats
Given the complexity of driver interactions and knowledge gaps, there are some aspects that could not be considered in our study.
For example, only the knowledge about interactions between CO 2 and temperature or light were considered robust enough in Seifert et al. (2020) to be implemented into our biogeochemical model.
Other dual-driver interactions, e.g., those including nutrients, are also important, although they could not be underpinned by a comprehensive statistical analysis. For instance, it was shown that the nutritional status of a phytoplankton cell can change its response to other environmental drivers (e.g., Van de Waal & Litchman, 2020). Stronger nutrient limitation suppresses the temperature dependence of metabolic rates in several species (Marañón et al., 2018) and modulates the response to ocean acidification (Eberlein et al., 2016) and high light conditions (Alderkamp et al., 2012) with likely implications on phytoplankton community composition (Alvarez-Fernandez et al., 2018;Flynn et al., 2015). Because changes in nutrient limitation are one of the most important drivers for changes from present-day to future phytoplankton communities (Dutkiewicz et al., 2013;Laufkötter et al., 2015;Marinov et al., 2010;Nakamura & Oka, 2019) it is important to gain a better understanding on the interaction between nutrients and other drivers to allow for the implementation into models (Van de Waal & Litchman, 2020). Besides, processes other than phytoplankton growth, such as calcification, are affected by driver interactions (e.g., Sett et al., 2014), but the number of studies investigating interactive driver effects on these processes is lower than on growth rates (Harvey et al., 2013). On top of that, driver interactions are not restricted to two drivers but can include several drivers that all interactively affect any phytoplankton's growth sensitivity to varying environmental conditions (Griffen et al., 2016;Ye et al., 2021). Following up on the assessment of mean climate states (i.e., present-day vs. future), transient simulations could furthermore provide additional information about decadal changes in phytoplankton responses to interacting drivers.
We simplified the mode of interaction by a linear interpolation between growth responses at low and high control levels, despite acknowledging that interactions are likely non-linear and should rather be described by a three-dimensional landscape (Boyd et al., 2018;Van de Waal & Litchman, 2020). However, given that most multipledriver laboratory studies assess only two levels per driver (Kreyling et al., 2018), linear interpolation between two points is the most prudent approach (Babyak, 2004). More laboratory experiments investigating multiple levels of each driver would be required to develop three-dimensional reaction norms (i.e., response surfaces) for models based on mechanistic understanding (Collins et al., 2022). These elaborate and time-consuming experiments should focus on key species of large-scale ocean biomes (e.g., Fay & McKinley, 2014) to allow for group-and region-specific considerations of driver interactions.
The increasing availability of omics-data (i.e., genetic data and data on gene expression; e.g., Chust et al., 2017) may constitute a solution to the limitations by laboratory experiments. Implications of interactive drivers on the genetic and transcriptomic level could then be used by cell-based models (e.g., Clark et al., 2011;Toseland et al., 2013), which in turn inform large-scale models about the mechanistic basis of plankton community shifts on a global scale . In addition, the modification of phytoplankton responses to bottom-up effects by evolutionary processes could be simulated by adaptive models (Le Gland et al., 2021). Although the process understanding to translate omics-data into changes in biological rates, e.g., growth, is still limited , complementing PFT models by a combination of different model types (e.g., cell-based, trait-based and adaptive models; D'Alelio et al., 2019) may be an emerging and promising field of biogeochemical modeling.

| Conclusion
In summary, our study shows that (i) by the end of the century, global phytoplankton biomass decreases by almost the same amount (5%-6%) in simulations with and without driver interactions in a high-emission scenario.
However, considering driver interactions alters the future projections of phytoplankton biomass on a regional and groupspecific level (e.g., by up to 65% for coccolithophores; Figures 6b and 8).
(ii) the biomass shifts between phytoplankton groups in turn lead to an altered community composition ( Figure 8). Specifically, diatoms will belong to the 'losers', whereas coccolithophores and small phytoplankton will be 'winners' of the anticipated changes within the future ocean when considering interactive effects.
Essentially, processes that are identified to be important in field and laboratory studies (e.g., interactive driver effects) should be implemented into models, and modellers, in turn, are in need of a database that is comprehensive enough to simulate processes for an entire phytoplankton community on the global scale (e.g., Chust et al., 2017;Davidson, 2014;Everett et al., 2017;Lombard et al., 2019;Moore, 2022). In addition, experimental set-ups using shipboardbased mesocosms on natural plankton communities  can test our finding that phytoplankton communities in the southern high latitudes are more susceptible to be affected by driver interdependencies than communities in the tropics and subtropics. This is especially important in light of the role of Southern Ocean diatoms in sustaining an efficient export flux of carbon and silicate. A continued collaboration between both fields will facilitate the improvement of models used for future projections of phytoplankton.

ACK N OWLED G M ENTS
MS and JH acknowledge funding from the Initiative and Networking The work reflects only the authors' views; the European Commission and their executive agency are not responsible for any use that may be made of the information the work contains. The authors thank two anonymous reviewers for their constructive comments on this work. Open Access funding enabled and organized by Projekt DEAL.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have declared no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available