We constructed a model of a spherical diatom cell to quantify its acquisition rates of CO₂ and/or HCO₃^– and to solve for the optimal allocation of metabolic energy between HCO₃^– uptake and biomass synthesis (i.e., growth) over a day (Figure 1). Given cell radius R (μm) and CO₂ concentration in bulk water [CO₂]_bulk (μM), the model simulates a daily growth rate g /day:

Note that although the model uses the cell radius to simplify calculations, the results are conventionally reported using the cell diameter.

The model first calculates the cell-specific acquisition rate of CO₂ under given R and [CO₂]_bulk. The model then counted the CO₂ leakage from the cell. The model reaches a maximal growth rate if the net acquisition of CO₂ is sufficient; otherwise, it lowers the growth rate to save metabolic energy for actively acquiring HCO₃^– as a carbon complement. The lowered growth also reduces the total carbon demand, partly relieving the requirement for HCO₃^–. Note that we set a minimum HCO₃^– acquisition at 20% of the total carbon demand even if the CO₂ supply is sufficient, reflecting findings that even very small diatoms also take HCO₃^– (e.g., Burkhardt et al., 2001; Shi et al., 2019). Therefore, there exists a solution for the modeled growth rate at which cellular energy is optimally allocated between HCO₃^– acquisition and growth.

The parameter values of the model (Supplementary Table 1) were set up based on the literature as described below. Our model only represents fundamental metabolic processes of cell growth and focuses on the responses of diatom growth to increasing CO₂, particularly the dynamics of HCO₃^– acquisition as the secondary inorganic carbon source, while ignoring other complex physiological processes.

Our model scheme of the CO₂ supply rate to the cell (F_T, μmol/cell/day) is established by improving those of previous studies (Gavis and Ferguson, 1975; Wolf-Gladrow and Riebesell, 1997; Reinfelder, 2011) in which both CO₂ diffusion and CO₂ dehydrated from HCO₃^– near the cell surface are included

where 3.456 × 10^–10 = 4 × 10^–15 L ⋅μm^–3 × 86,400 s/day, in which the constant of 4 is from the original equation (Gavis and Ferguson, 1975), 10^–15 is the factor to convert μm³ to L for consistency with the unit used in CO₂ concentrations, and 86,400 is the factor to convert day to second; D (μm²/s) is the diffusivity of CO₂ in water k′ (s^–1) is a combined rate constant for the hydration of CO₂ with H₂O and OH^–, with its value estimated at 15°C, pH of 8.1 and salinity of 35 PSU (Supplementary Table 1; Wolf-Gladrow and Riebesell, 1997; Reinfelder, 2011), and [CO₂]_sfc (μM = μmol L^–1) is the seawater CO₂ concentration near the cell surface. Although [CO₂]_sfc can be effectively eliminated from the equation by assuming that it was 1/3 of [CO₂]_bulk, as done in a previous study (Reinfelder, 2011), [CO₂]_sfc is likely to decline with increasing cell size (Flynn et al., 2012). In this study, we further solve [CO₂]_sfc by assuming a balance between CO₂ supply to the cell surface (F_T) and CO₂ permeation across the cell membrane (F_P, μmol/cell/day):

where P (μm/s) is the CO₂ permeability coefficient of the cell membrane, Δ_CO2 = 10% is the relative reduction in CO₂ concentration from near the cell surface to the cytoplasm (Hopkinson et al., 2016), that is, Δ_CO2[CO₂]_sfc is the difference in CO₂ concentration between the two sides of the cell membrane, and 3.456 × 10^–10 = 4 × 10^–15 L ⋅μm^–3 × 86,400 s/day. From Eqs. 2 and 3, [CO₂]_sfc can be solved as a function of input variables R and[CO₂]_bulk:

where D′=D⋅(1+R⁢k′/D). This equation generates a declining [CO₂]_sfc with increasing cell size, as expected (Supplementary Figure 1). F_T can then be obtained by substituting Eq. 2 into either Eq. 3 or Eq. 4.

This scheme introduces a new parameter P, for which the maximal value of 560 μm/s found in experiments of several diatom species (Hopkinson et al., 2011) is used in the model (Supplementary Table 1). The assumption of F_T=F_P used in the scheme (Eq. 3) implies that CO₂ supplied to the cell surface is completely acquired by the cell, and the CO₂ supply rate is solved at its maximal potential, which we annotate as F_T,max.

The energetic cost of growth is estimated by adding the theoretical requirements of ATP when synthesizing particulate organic carbon (POC) from CO₂ [3 mol ATP (mol C) ^–1] (Raven, 1991) and synthesizing particulate organic nitrogen (PON) from nitrate [2 mol ATP (mol N) ^–1] (Eichner et al., 2014), which is equivalent to ∼ 0.30 ATP (mol C) ^–1 by using a Redfieldian molar C:N ratio of 106:16. By further assuming small additional energy of 0.2 mol ATP (mol C) ^–1 used in other cellular processes, the energy cost rate for growth, e_g, is set at 3.5 ATP (mol C) ^–1 (Supplementary Table 1).

The intracellular energy production rate E₀ [ATP/mol C/day] is set at:

where the maximal daily average growth rate g₀ = 3.0/day represents those maximal rates found in previous experiments culturing diatoms under saturating light and nutrients (Sarthou et al., 2005). The instantaneous maximal growth rate (and consequently, E₀) is further set to vary in a 12:12 h light–dark cycle in which the rate peaks for 6 h around noon to represent light saturation, while its daily average is kept at g₀ (Supplementary Figure 2). This setup can improve the model performance to mimic the diel variation of inorganic carbon demand. For example, a model cell that would have sufficient CO₂ supply, should its maximal growth rate be set constant at g₀, may become, after implementing the diel cycle, CO₂-limiting and have to acquire HCO₃^– around noon when the instantaneous maximal growth rate is 8/3 times higher than g₀ (Supplementary Figure 2).

ATP expenditure by CCM is generally high in diatoms (Raven, 1991), which may be partially satisfied by the Mehler reaction (Behrenfeld et al., 2008). The energy cost for diatoms to transport 1 HCO₃^– molecule has been estimated at 0.5 ATP to cross the cytoplasmic membrane (Liu et al., 2017) and 1 ATP to cross the chloroplast membrane to reach pyrenoids (Raven et al., 2000). Therefore, we set the energy cost rate for HCO₃^– acquisition, e_bc, at 1.5 ATP (mol C) ^–1 (Supplementary Table 1).

The net inorganic carbon acquisition rate F (μmol/cell/day) is

where F_bc (μmol/cell/day) is the HCO₃^– acquisition rate, and l_k = 30% is the fraction of acquired carbon leaking from the cell based on a previously reported CO₂ efflux for diatoms (Burkhardt et al., 2001). Note that the leakage in the model does not count the exudation, which, however, could be small (∼5% of primary production) in exponentially growing phytoplankton (Nagata, 2000). The model sensitivity to l_k was also tested (see below).

If both F_T and F_bc are determined (discussed below), F can be solved (Eq. 6), and the growth rate is then

where q_c is the cell carbon quota (μmol cell^–1) estimated by adopting an empirical relationship with the cell volume of diatoms (Menden-Deuer and Lessard, 2000):

F_T and F_bc are optimally solved under the constraint of intracellular energy production E₀. The intracellular energy requirement is:

in which

A trial rate of HCO₃^– acquisition (F_bc) is first set to be the minimal fraction (fm⁢i⁢nb⁢c) of total carbon uptake (F_T,max + F_bc) so that:

Then, a trial energy requirement E is calculated using Eqs. 10–12 and setting F_T = F_T,max (see section “Model Scheme of the CO₂ Supply”). When E > E₀, energy production limits growth and the total carbon supply is more than sufficient. Therefore, the effective F_T and F_bc must be smaller than their trial rates while their ratio holds. Therefore, a factor of E₀/E is multiplied by F_T,max and the trial F_bc to obtain the effective F_T and F_bc, which ensures E = E₀.

However, a trial rate of E larger than E₀ indicates that the inorganic carbon supply is insufficient, additional acquisition of HCO₃^– is needed, and the effective F_T is at its maximal allowed rate F_T,max. The effective F_bc is calculated from Eqs. 9–11 by setting E = E₀ (Eq. 9) and F_T = F_T,max (Eq. 10).

Once the effective F_T and F_bc are determined in either case, the growth rate g can be calculated from Eqs. 6, 7. Noting that the intracellular energy produced is fully used, the scheme described here solves an optimal allocation of intracellular energy to estimate g at its maximal potential.

The cell growth rate was simulated at two levels of bulk-water CO₂ concentrations of 10 μM (low carbon condition, LC) and 20 μM (high carbon condition, HC) over a diameter range of 2–200 μm. LC was selected at 10 μM to represent the typical condition of subtropical surface ocean waters. GRR is calculated as the ratio of the simulated growth rates under HC and LC. We also calculated two key indicators, the value of the peak GRR over the simulated diameter range and the corresponding cell diameter at the peak GRR (namely, the optimal cell diameter).

Our model does not represent other factors, such as light and nutrient concentrations, that can limit the diatom growth rate. Instead, the model parameter g₀ can be considered to represent collective effects from multiple limiting factors. We then conducted model experiments by setting g₀ at 0.1, 0.5, 1.0, and 2.0/day.

We also compared the model sensitivity, in terms of the peak GRR and the optimal cell diameter, to 6 model parameters, including P, l_k, e_bc, g₀, △_CO2, and fm⁢i⁢nb⁢c. These parameters were tested in a range of ± 50% of their default values. The model was also tested over 12–30 μM of the input variable [CO₂]_bulk in HC.

We collected 85 pairs of growth rate data from published CO₂ enrichment experiments culturing diatoms, among which 5 pairs of data were removed due to extremely high culturing CO₂ concentrations (Supplementary Table 2). Due to different magnitudes of CO₂ enrichment, we linearly adjusted the originally reported GRR (GRR_orig) cultured under high ([CO₂]_H) and low ([CO₂]_L) CO₂ concentrations to GRR_adj at a CO₂-enrichment factor (f_CO2) of 200%:

The cell diameters of the diatoms were either reported in the experiments or estimated using the average diameters of the same species collected in a global diatom database (Leblanc et al., 2012a). Considering the purpose of this study, for nonspherical diatom species, we simply used the length of their shortest dimension as the cell diameter because it was the shortest distance that inorganic carbon was transported to the cell center, admitting that this simplification slightly overestimated the CO₂ supply rate.

We projected the peak GRR and optimal cell size of diatoms to CO₂ increase in global surface seawater over the 21st century. We first obtained projected results from the Community Earth System Model (CESM) 1.0-BGC module under the RCP8.5 scenario, including the annual climatology of CO₂ partial pressure, pH, temperature, diatom carbon fixation rate, and diatom carbon biomass during the historical (1990–2000) and future (2090–2100) periods¹. The CO₂ partial pressure, pH, and temperature were used to calculate the CO₂ concentration using CO2SYS (Lewis and Wallace, 1998). The diatom growth rates were calculated as the ratio of the diatom carbon fixation rate to the diatom carbon biomass.

The cellular model ran using historical and future CO₂ concentrations to estimate monthly peak GRR and optimal cell size at a resolution of 1°× 1° in the global ocean, in which historical diatom growth rates were used to set location specificity g₀. Because CESM did not simulate diatom sizes, our results can only be regarded as a demonstration of the spatial variations in the responses of diatoms impacted by the magnitude of CO₂ increase and environmental favorability (such as nutrients, temperature, and light), with the latter collectively represented here by the CESM-projected diatom growth rate.

The scientific color map batlow was used in the plot of the global prediction to prevent visual distortion of the data and exclusion of readers with color-vision deficiencies (Crameri et al., 2020).

The modeled cell growth rate declines with increasing cell diameter under LC, while the modeled growth under HC remains at a maximal rate of 2.67/day before it starts to decline when its diameter is larger than 7 μm (Figure 2A). As the ratio of the modeled growth rate under HC to that under LC, the GRR increases with cell diameter, reaches a peak value of 1.14 at an optimal cell diameter of 7 μm, and then gradually decreases to 1.02 at the modeled maximal cell diameter (Figure 2B).

We then analyze the inorganic carbon (C_i) budget of the model cell in two cell diameter ranges. As expected, the model simulates a decreasing CO₂ supply rate with increasing cell diameter (Figure 3). First, for cells with a diameter < 7 μm, CO₂ and minimum HCO₃^– supply under LC are insufficient to meet the requirements of both growth and leakage, causing HCO₃^– acquisition to increase with cell size (Figure 3A). However, under HC, CO₂ and minimum HCO₃^– supply are sufficient, and cells do not acquire additional HCO₃^– (Figure 3B). Second, for cells with a diameter > 7 μm, CO₂ and minimum HCO₃^– supply become insufficient even under HC (Figure 3B), and the modeled HCO₃^– acquisition increases with cell size under both HC and LC (Figure 3).

The modeled pattern between GRR and cell size (Figure 2B) is mainly determined by the magnitude of energy savings on HCO₃^– acquisition from LC to HC (Figure 4). In the cell diameter range < 7 μm, the energy expenditure on HCO₃^– increases greatly with increasing cell size under LC but is unchanged under HC, resulting in substantial energy savings that increase with cell size (Figure 4). In the cell diameter range > 7 μm, however, the energy used in HCO₃^– acquisition increases with cell size under both HC and LC, leading to a narrowed difference between the two (Figure 4). This can be explained by the decrease in the relative contribution of CO₂ to the total C_i acquisition with increasing cell size (Figure 3): The fraction of HCO₃^– acquisition that can be substituted by the same magnitude of the elevated CO₂ supply decreases in larger diatoms. In other words, larger cells become less sensitive to the CO₂ increase.

The model reveals a general pattern in which (1) the smallest diatoms respond weakly to the increase in seawater CO₂ because CO₂ supply fulfills most of the cell demand even without that increase, (2) the response intensifies with a moderate enlargement in cell size, while (3) large cells are insensitive to the CO₂ increase because CO₂ supply is a small contribution to total C_i acquisition. There exists an optimal, intermediate cell size in which the modeled cell growth rate is enhanced the most by the same magnitude of CO₂ increase. Our model shows that the optimal cell size corresponds to a point where the model cell under HC just starts to use HCO₃^–.

The maximal growth rate (g₀) has been set at a high level assuming that the modeled diatom grows at optimal conditions. We then conduct the model experiments with lowered g₀, representing that the growth of the modeled diatom is limited by other factors. The results show that the modeled diatom starts to benefit from the CO₂ increase at larger cell sizes when g₀ decreases, so that the optimal cell size also increases (Figure 5A). This result is expected because the modeled diatom needs less inorganic carbon at a lowered g₀, and the CO₂ supply alone becomes sufficient for larger diatoms. At very low g₀ values, such as 0.1/day, the CO₂ supply is even sufficient at the maximal cell diameter simulated in our model (200 μm), and diatoms do not respond to the CO₂ increase in the whole model domain (Figure 5A). Additionally, g₀ does not impact the modeled peak GRR (Figure 5A), the reason for which is discussed below. These results indicate that the optimal cell size of diatoms in response to the CO₂ increase tends to be small when the environmental conditions favor their growth but can be larger when other factors are also limiting.

We also compared the sensitivity of the modeled peak GRR (Figure 5B) and optimal cell size (Figure 5C) to various model parameters. The CO₂ supply rate, mainly determined by the production of the CO₂ permeability (P) and CO₂ gradient across the cytoplasmic membrane (△_CO2) (Eqs. 3, 4), has a positive relationship with the optimal cell size but does not impact the GRR. This can be more clearly understood from a model experiment with a 50% higher P (Supplementary Figure 3A): The cell needs to use additional HCO₃^– at large cell sizes with increased CO₂ supply. Hence, the curves of HCO₃^– acquisition energy under both HC and LC move toward larger cell sizes without noticeable changes in the magnitude of the energy cost, leading to a moderately elevated optimal cell size but unchanged peak GRR (Figures 5B,C; Supplementary Figure 3A). Similarly, a smaller maximal daily growth rate (g₀) also leads to a larger optimal cell size but no effect on peak GRR (Figure 5 and Supplementary Figure 3B) because the lowered carbon demand also allows larger cells to only use CO₂. In contrast, reduced CO₂ leakage (l_k) lowers the requirement of HCO₃^– acquisition and shrinks the difference in energy cost between HC and LC, causing a decreased GRR and enlarged optimal cell size (Figure 5 and Supplementary Figure 3C).

The elevated energy cost rate of HCO₃^– acquisition (e_bc) substantially increases the peak GRR but has little effect on the optimal cell size (Figures 5B,C) because it mainly determines the energy allocation and not the cell size at which HCO₃^– acquisition is needed. The high bulk-water CO₂ concentration of HC (HC[CO₂]_bulk) greatly increases both the peak GRR and the optimal cell size (Figures 5B,C), which is expected because the cell under HC benefits more from higher CO₂ while starting to use HCO₃^– at a larger size. The minimum HCO₃^– uptake proportion (fm⁢i⁢nb⁢c) slightly impacts both the peak GRR and the optimal cell size: a higher fm⁢i⁢nb⁢c leads to a reduced model sensitivity to the change in CO₂ and results in a lower peak GRR and a larger optimal cell size.

Our model sensitivity tests reveal that the peak GRR is largely determined by the magnitude of the seawater CO₂ increase and the energy consumption rate in HCO₃^– acquisition. Meanwhile, the optimal cell diameter increases with increasing CO₂ supply (higher membrane permeability to CO₂ and larger increase in seawater CO₂) but decreases with cellular carbon demand, such as higher growth and leakage rates.

Our theoretical model reveals a unimodal relationship between the GRR and diatom cell size. Although the increased CO₂ concentration and acidified seawater can have multiple physiological effects on diatoms, our model considers the most direct effect: the variations in species and rates of inorganic carbon acquired by diatoms. The absolute rates of growth do not solely determine the competition among phytoplankton because they can also be controlled by other factors, such as zooplankton grazing and viral lysis, which is also suggested by the fact that phytoplankton with different growth rates coexist in the ocean. A phytoplankton species with a stronger enhancement in the growth rate from the CO₂ increase, can have a better advantage and reach higher biomass when competing with other phytoplankton, which is also confirmed by a modeling study (Dutkiewicz et al., 2015). Therefore, the unimodal relationship between GRR and diatom cell size shown by our model suggests a competitive advantage for middle-sized diatoms of ∼7 μm in the future high-CO₂ ocean.

It is generally true that organisms limited more strongly by a resource should benefit more from its repletion, but only if there are no substitute resources. CO₂, for example, is not the only carbon resource for phytoplankton and can be complemented by, although not preferentially, HCO₃^– when CO₂ is insufficient. This is the important reason why our model reveals that middle-sized diatoms can have a competitive advantage over larger diatoms from elevated CO₂, even though the latter is limited more strongly by CO₂.

We also explored the GRR–cell size relationship using historical data from CO₂ manipulation experiments of diatoms. The originally reported GRR in those experiments was first interpolated to the same degree of CO₂ enrichment (Supplementary Figure 4 and Supplementary Table 2, see section “Materials and Methods”). The adjusted GRRs varied greatly (0.7–1.4) and did not support the previous speculation that the GRR should increase with cell size (Spearman’s correlation, p = 0.59) (Figure 6). The GRRs also did not directly show the same pattern as that proposed by our model. However, all GRRs higher than 1.2 were only found in cell diameters of 3–7 μm, while the GRRs at other cell diameters were mostly lower than 1.1 (Figure 6). This may indicate that the growth rate of diatoms in this middle cell size range had the greatest potential to respond to increasing seawater CO₂. Additionally, the GRRs in the same size range also varied the most (Figure 6), suggesting that the GRR of middle-sized diatoms can be impacted by different factors and mechanisms and/or be species specific. Interestingly, considering the species in the genus Thalassiosira that had the largest cell diameter range (3.9–47 μm) in our dataset, their GRRs decreased with cell diameter (Spearman’s correlation = −0.45, p < 0.05), with GRRs in cell diameter ≤ 7 μm (1.14 ± 0.13, mean ± s.d.; n = 11) significantly higher than those in cells > 7 μm in diameter (1.00 ± 0.06, mean ± s.d.; n = 10) (t-test, p < 0.01). Nevertheless, these historical data showed some consistency with our model, while the GRR can be a synergistic result impacted not only by cell size but also by other physiological and ecological characteristics and culturing conditions. More studies are needed to reveal complex mechanisms controlling the GRR.

Other studies may directly or indirectly support our model. Hancock et al. (2018) cultured a near-shore Antarctic community over a gradient of CO₂ levels, showing that the abundance of nanosized diatoms (2–20 μm in diameter) increased with CO₂, while larger microplanktonic diatoms (>20 μm in diameter) and smaller discoid centric diatoms of 1–2 μm in diameter had no significant response to CO₂. Other studies also showed that middle-sized diatom species (∼8 μm) had an increased contribution to the community under increased CO₂ (Hoppe et al., 2013; Eggers et al., 2014).

We also constructed a new scheme of the maximal seawater CO₂ supply rate to diatom cells as a function of cell size. Compared to other studies (Milligan et al., 2009; Reinfelder, 2011), the new scheme solves a cell size–dependent CO₂ concentration near the cell surface (Supplementary Figure 1). Our scheme introduces a new constant parameter, the CO₂ permeability of the cell membrane. Although elevated seawater CO₂ can lead to a saturated CO₂ supply and reduce the effective CO₂ permeability (e.g., Sultemeyer and Rinast, 1996), our scheme is constructed to study the maximal effect of seawater CO₂ increase and therefore uses the maximal potential of CO₂ permeability, which theoretically is determined by the number of channels such as pores and aquaporins (Matsui et al., 2018; Blanco-Ameijeiras et al., 2020) allowed per unit area of the cell membrane and should be independent of cell size. Limited data on diatom CO₂ permeability (Hopkinson et al., 2011) also do not show a clear relationship to cell size. Our model scheme therefore can be more suitable to quantify the relationship between phytoplankton cell size and their maximal potential for CO₂ uptake.

The growth rates of diatoms are generally inversely related to cell sizes (Cosper, 1982; Sarthou et al., 2005; Marañón, 2015), which could result from a lower CO₂ supply rate to larger diatoms, as observed in this study, and/or other size-dependent physiological characteristics of diatoms. This implies that the constant model parameter of maximal growth rate g₀ can in reality decline with increasing cell size. We then conducted a model experiment by using an empirical function for diatoms (Sarthou et al., 2005):

where V was the cell volume (μm³). In the results, the unimodal pattern between GRR and cell size remained (Supplementary Figure 5), although the optimal cell size increased to 45 μm, which was expected because of the negative relationship between the optimal cell size and g₀ (Figure 5B).

As eCA can be potentially important in particularly large phytoplankton in converting HCO₃^– to CO₂ extracellularly (Martin and Tortell, 2008; Tortell et al., 2008; Chrachri et al., 2018) and therefore saving energy on HCO₃^– acquisition, we tested its role by increasing the conversion rate of HCO₃^– to CO₂ in the model by an enhancement factor (f_eCA) of 2–10 (Tortell et al., 2008). As the ratio of the chemical conversion between HCO₃^– and CO₂ to the CO₂ diffusivity is R⁢k′/D (Eq. 2) (Reinfelder, 2011), the CO₂ supply (F_T) was then enhanced to a new rate:

In the model experiments, although the GRR was elevated in the large cells, the unimodal pattern between the GRR and cell size, the optimal cell diameter, and the magnitude of peak GRR were barely changed (Supplementary Figure 6A). We also conducted another model experiment by linearly increasing f_eCA from 1 (i.e., no enhancement) at the smallest cell diameter (2 μm) to 10 at a cell diameter of 200 μm to mimic a hypothetical scenario of stronger eCA of larger cells, and still obtained a similar pattern between the GRR and cell size, except that the GRR remained at a constant low level when the cell diameter was >100 μm (Supplementary Figure 6B).

The two parameters D and k′ used in the CO₂ reaction–diffusion kinetics (Eq. 2) are fixed in our model, while in reality, they change with temperature (Wolf-Gladrow and Riebesell, 1997). A sensitivity test of temperature-dependent D and k′, however, showed no obvious effect on the model results (Supplementary Figure 7).

The above experiments can partly support the robustness of our model results, suggesting that the negative correlation between diatom growth rates and cell size, the extracellular conversion of HCO₃^– to CO₂, and temperature do not change the unimodal relationship between the GRR and the cell size of diatoms.

Our model sensitivity experiments reveal that the modeled optimal cell size enlarges with stronger seawater CO₂ supplies associated with its higher concentration and permeability across the cell membrane, but diminishes with elevating cellular carbon demand determined by higher rates of cell growth potential (Figure 5). Among these factors, the magnitude of future CO₂ increase and the maximal growth rate of diatoms can vary substantially in the global ocean. Particularly, as already discussed, our model does not represent the effects of other factors that can limit diatom growth. We then ran our model in the global ocean using the projected levels of these two variables in this century under the RCP8.5 scenario from an Earth system model (see section “Materials and Methods”), in which multiple limiting factors on diatom growth were simulated. The spatial variation in the modeled peak GRR was generally small in the global ocean, with the lowest value in the tropical Pacific upwelling zone (Figure 7A) associated with its smallest magnitude of CO₂ increase (Figure 7C). The modeled optimal cell size, however, varied greatly from < 20 μm in tropical oceans to the highest modeled cell diameter (i.e., GRR increases monotonically with cell size) in the subpolar regions (Figure 7B). The modeled optimal cell size was systematically smaller than what was predicted in our theoretical model (Figures 2, 5B) because the CESM-simulated diatom growth rates incorporated other limiting factors and were mostly below 1.0/day (Figure 7D), much lower than those found in laboratories with ideal culturing conditions (Sarthou et al., 2005). The lowest optimal cell size was also in the tropical Pacific upwelling zone (Figure 7B), contributed synergistically by the high growth rates of diatoms and low magnitude of CO₂ increase (Figures 7C,D).

Our projection for diatoms in the global ocean should be regarded at most as the first-order estimation. However, some interesting implications emerged from the projection. Diatoms are one of the important contributors to carbon export to the deep ocean (Falkowski et al., 2004), particularly in tropical oceans (Siegel et al., 2014, 2016). The relatively small optimal cell size projected in these regions (Figure 7B) implies that, without considering other effects, the size structure of the diatom community would not shift to larger species in the future high-CO₂ ocean, which consequently would not lead to an increase in the sinking speed of organic particles and carbon export.

Carbonic anhydrase in diatoms can convert CO₂ to HCO₃^– in the cytoplasm and result in much-increased HCO₃^– transport further into chloroplasts, as found for diatom Phaeodactylum tricornutum (Hopkinson et al., 2011; Hopkinson, 2014). Additionally, a high leakage of CO₂ from the chloroplast to the cytoplasm has been suggested (Hopkinson et al., 2011), although it is unclear how much of the leaked CO₂ further leaks to the environment or is converted to HCO₃^– in the cytoplasm and transported back to the chloroplast. The CO₂ leakage from the cell to the environment was already counted in the model. However, if there is substantial recycling of the leaked CO₂ to HCO₃^– in the cytoplasm, the rate of HCO₃^– transfer across the chloroplast membrane can be further elevated. In other words, e_bc can possibly be higher for at least some diatoms, leading to higher GRRs (Figure 5A) but small changes in the optimal cell size (Figure 5B).

There are some limitations in our model that can be investigated in future studies. First, our model does not include the negative effects of the lower pH concomitant with elevated CO₂, which can be important to some species or functional groups of phytoplankton (e.g., Taylor et al., 2012; McMinn et al., 2014; Cyronak et al., 2016; Hong et al., 2017; Gafar et al., 2018; Luo et al., 2019; Shi et al., 2019). If there is no relationship between the negative effects of lower pH and diatom cell sizes, which, however, is unclear, our modeled pattern between GRR and diatom cell sizes can still hold.

Second, RuBisCO in different diatoms shows greater variation in the Michaelis constant for CO₂ (23–68 μM) (Young et al., 2016), indicating diverse patterns of CCM, RuBisCO kinetics, and catalytic features in diatoms. Our model does not separate diatom species but focuses on the general response of diatoms of different sizes to seawater CO₂ increases. Nevertheless, the variation of these factors can change model parameters of the cell membrane permeability to CO₂ (P) and the energy cost rate in HCO₃^– uptake (e_bc), among which a ± 50% change in the former can vary the optimal cell size in a range of 2–13 μm while that in the latter has no substantial effect (Figure 5C). That is, if the CO₂ permeability is substantially lower than that used in our model, the GRR can be highest in the smallest diatoms and then decreases with increasing cell size.

Third, light appears to interact with elevated CO₂ in determining the collective effect in cell size–dependent responses of diatoms, such as that we already discussed for the stronger response of larger diatoms to CO₂ increase under very high light. Another study also proposed a light-based maximal benefit of CCM in middle-sized phytoplankton: The depletion of light near the center of large cells depresses photosynthesis and reduces the energy available for CCM, making CCM economically less efficient in large cells (Wirtz, 2011). This effect and our proposed mechanism, if valid, can synergistically intensify the competitive advantage of middle-sized diatoms under the CO₂ increase.

Last, cultured diatoms also often respond differently to CO₂ increases depending on the availability of nutrients such as nitrogen (Li et al., 2012), phosphorus (Sun et al., 2011), and silicon (Tatters et al., 2012), while these nutrients are not simulated in our model. Instead, our model experiments implicitly represent the limitation of these nutrients by applying a lower model parameter g₀, showing that middle-sized diatoms are more likely to have better competition from the CO₂ increase when these nutrients are not severely limiting (Figures 5, 7). Nevertheless, our model does not simulate the cell size–dependent limitation of nutrients on diatoms. In the future warmer and more stratified oceans, nutrients can become more limiting, and smaller diatoms may gain another competitive advantage because of their lower surface-to-volume ratios. Considering this effect, the increasing trend of GRR from small to middle-sized diatoms can be weakened or even reversed. However, the decreasing trend of GRR from middle-sized to large-sized diatoms may be even stronger.

Overall, despite the limitations of our model, the discussions above suggest that the GRR in the future high-CO₂ ocean, in general, likely decreases from middle-sized to large-sized diatoms, except in the regions where the light is substantially high or nutrients are substantially limiting on diatoms. The GRR may also increase from small to middle-sized diatoms unless the cell membrane permeability to CO₂ is low or nutrients become severely limiting, particularly on middle-sized diatoms.

In summary, this study constructed a theoretical model and revealed a competitive advantage of middle-sized diatoms of ∼7 μm over both small and large diatoms when seawater CO₂ increases. Although the model is mostly based only on the basic principles of the carbonate equilibrium system and inorganic carbon diffusion with constraints on cellular energy allocations, it appears to be a robust pattern, although the exact cell size of which diatoms benefit most from the CO₂ increase can change substantially. If the mechanism proposed in this study dominates, the diatom community in the future high-CO₂ ocean may not shift toward larger cells, and carbon export by sinking diatoms may not be substantially stimulated, while we admit that other known and unknown factors can also influence the cell size–dependent responses of marine diatoms to the CO₂ increase. Our simplified physiological model processes, or even more simply the size-dependent pattern revealed by our model, can be integrated into marine ecosystem models to improve the predictions of the size composition of diatoms, phytoplankton community structure, and productivity for future high-CO₂ oceans.