Effect of nutrient supply on cell size evolution of marine phytoplankton

1.   Introduction

Phytoplankton is an important functional group in the ocean and plays a central role in marine ecosystems. As a kind of photoautotrophic microbes, phytoplankton is responsible for about half of global primary productivity ^[1]. Over the past several decades, global warming has asserted significant effects on the coastal environment and marine organisms ^[2]. Recently, ecologists found that marine phytoplankton meets Bergmann's Rule, that is, increasing temperature can lead phytoplankton cells to become smaller ^[3]. Extensive studies have been carried out to explore why phytoplankton becomes smaller with increasing temperature. It was found that an increase in physiological and metabolic rates directly caused by increasing temperature may lead to the reduction of phytoplankton size. However, the increased physiological and metabolic rates are not necessarily the main reason for the reduction in phytoplankton size ^[4]. It has been found that reduced nutrient availability indirectly caused by climate warming is a dominant factor for the decrease in phytoplankton cell size ^[5], as well as phytoplankton biomass ^[4].

The availability of nutrients plays a crucial role in determining the biomass and productivity of phytoplankton. Recently, some studies have found that the global phytoplankton biomass is declining at a rate of approximately 1% per year in numerous ocean areas, especially in the south and equatorial Atlantic regions ^[6]. Many pieces of evidence indicate that declines in phytoplankton production and biomass are associated with nutrient limitation ^[4]. Owing to climate change, the sea surface temperatures in most oceans have increased over the past century ^[2], which has led to the strengthening of the density gradient and enhancement of thermal stratification. Consequently, the vertical mixing rate between the upper layer and the deep ocean declines, and the nutrient supply becomes limited ^[7]. Nutrients are gradually depleted in the mixed layer, which exerts an impact on the growth and survival of phytoplankton.

In addition to limiting the phytoplankton growth rate, the variability of nutrient can also be a determinate factor in the cell size evolution of phytoplankton. As observed in most oceans, small phytoplankton tends to dominate in nutrient-poor waters, while large phytoplankton is more abundant in nutrient-rich areas ^[8]. In fact, cell size is a key trait of phytoplankton, and it makes a significant difference in physiological and ecological functions. Specifically, cell size affects the growth rate ^[9], sinking rate ^[10], nutrient quota ^[11], and competitive ability ^[12]. Also, the variation of phytoplankton cell size has cascading consequences on higher trophic levels and even the entire marine ecosystem. Therefore, it is necessary to explore the evolutionary adaptation of phytoplankton cell size under varying nutrient availability.

Extensive research has focused on effects of nutrient supply on long-term behavior of phytoplankton ^[13,14]. Recent studies have shown that the phytoplankton species may respond to the change of nutrient supply by altering its physiological characteristics, such as the cell size or the elemental composition in the cells ^[15,16]. Especially, there is an obvious trend that the mean cell size of phytoplankton is decreasing when the nutrient becomes limiting in the mixed layer ^[17,18]. Marine ecologists have done some experimental and field studies to verify the relationship between the nutrient supply and phytoplankton cell size. Peter and Sommer ^[14] performed a laboratory experiment with three levels of nutrient limitation, and the nutrient stress is controlled by semicontinuous dilution. The results show that the nutrient supply has a crucial effect on the cell size of phytoplankton, and the size reduces significantly at the highest level of nutrient limitation.

However, the experimental studies are difficult to measure and evaluate the long-term impacts quantitatively. The trait-based model is an analytical method to understand the relationships between environmental factors and traits of phytoplankton ^[19,20]. Moreover, the adaptive dynamics method can be applied to the trait-based model to explore the evolutionary dynamics of the traits, which has been studied by many researchers ^{[21,22,23,24]}. The ecological model and theoretical analysis provide a good opportunity to reveal the dynamical mechanism of the population system ^[25]. In addition, the existing studies of the trait-based nutrient-phytoplankton model show that the cell size can only be an evolutionary stable strategy (ESS) in a nutrient-phytoplankton model ^[26]. They failed to find the evolutionary branching point, which implies the trait diversity of population. However, phytoplankton is extremely diverse with respect to its cell size in the marine ecosystem. Therefore, it is of considerable importance to explore whether the adaptive model admits an evolutionary branching point under varying nutrient levels.

The enhanced stratification structure of the ocean due to climate changes reduces the vertical mixing rate between the deep and the surface layers of the ocean. Hence the nutrient supply from the deep layer to the surface layer is limited, which significantly affects the growth and reproduction of phytoplankton. It is reasonable and important to examine the effects of the variation of nutrient supply on the evolution of biomass and cell size of phytoplankton. To achieve this objective, we combine the trait-based model and adaptive dynamics method to elucidate the responses of phytoplankton species to the variation of nutrient supply. The paper is organized as follows. In the next section, a trait-based nutrient-phytoplankton model is proposed, and some fundamental results about the ecological model are obtained. In Section 3, we develop an evolutionary model of the phytoplankton cell size and analyze the evolutionary dynamics of the singular strategy. Some numerical simulations are carried out to demonstrate the impacts of input nutrient concentration and vertical mixing rate on the persistence of phytoplankton species and the evolution of cell size in Section 4. In Section 5, we conclude our results and explain the potential biological meanings. Finally, the paper ends with the detailed proofs of our analytical results in Section 6.

2.   Ecological model and dynamics

2.1.   Model formulation

To explore the effect of the nutrient supply on the population and evolutionary dynamics of marine phytoplankton, we formulate a nutrient-phytoplankton model consisting of two state variables, including nitrogen concentration $N$ ($\mathrm{fmol\; N\cdot L^{-1}}$) and phytoplankton density $P$ ($\mathrm{cells\cdot L^{-1}}$). Assume that the ocean is composed of the mixed layer and the deep ocean (see Figure 1). The phytoplankton grows in the mixed layer, and the nitrogen is mainly supplied from the deep ocean. Moreover, the upper layer is assumed to be well mixed, thus the nutrient and phytoplankton are homogeneous there.

We develop a size-dependent nutrient-phytoplankton model under the following assumptions.

● If the phytoplankton cells are spherical, then the mean cell diameter $x$ is applied to describe the cell size; if the phytoplankton cells are not spherical, then the equivalent spherical diameter (ESD) can be used to characterize the cell size.

● The exchange rate between the mixed layer and the deep ocean is proportional to the vertical mixing rate $v$, but is inversely proportional to the mixed layer depth $z_m$. The amount of input nitrogen from the deep ocean is the product of the exchange rate $v/z_m$ and the input nitrogen concentration $N_{in}$. The water exchange also leads to some nitrogen being removed from the mixed layer.

● The proportion of recycled nitrogen from the phytoplankton debris is denoted by $r$ with $0 < r < 1$.

● The nitrogen is absorbed and utilized by phytoplankton, which depends on the phytoplankton nutrient quota $Q(x)$.

● The maximum growth rate of phytoplankton follows an allometric pattern, which means that the growth rate represented by $\mu(x)$ is a function of phytoplankton cell size.

● The nitrogen limitation on growth rate is characterized by Michaelis-Menten function $g(N)$ with half-saturation constant $K$, which is $g(N) = N/(N+K)$.

● The losses of phytoplankton are determined by four processes, including a constant size-independent mortality rate $m$, a size-dependent sinking rate $s(x)/z_m$, a dilution rate due to the water mixing $v/z_m$, and intraspecific competition $\alpha(0)$.

Based on the above assumptions, the interactions between the nitrogen and phytoplankton in the mixed layer can be described by

where $i = x, y$. Without loss of generality, let $i = x$.

From (2.1), we find that some ecological parameters of phytoplankton are size-dependent, such as the growth rate $\mu(x)$, sinking rate $s(x)$, and nutrient quota $Q(x)$. In addition, the competitive effect $\alpha(x-y)$ is also size-dependent. Note that all size-dependent functions are positive and bounded. Next, we describe other potential properties for these size-dependent parameters.

Growth rate $\mu(x)$. It is accepted that the maximum growth rate of phytoplankton is related to its cell size. Recent research shows that there exists a unimodal relationship between phytoplankton growth rate and its cell size, and the growth rate reaches its peak at an intermediate cell size ^[9]. It means that if the cell size is small, then the growth rate increases with the size; if the cell size is large, then the growth rate decreases with the size ^[20]. Hence, $\mu(x)$ satisfies the following inequalities

Nutrient quota $Q(x)$. The nutrient quota of phytoplankton is mainly determined by the protoplasm volume, which is approximately proportional to phytoplankton cell volume ^[19]. Hence, it is trivial to find that the nutrient quota increases with phytoplankton cell size, which is $Q'(x) > 0$.

Sinking rate $s(x)$. Phytoplankton cells sink at a certain rate since they are denser than the fluid around them. It is generally accepted that the larger cells sink faster than the smaller cells, which leads to $s'(x) > 0$, $s{''}(x) > 0$.

Competitive effect $\alpha(x-y)$. It is common ground that large phytoplankton is more competitive than small ones in terms of light capturing, nutrient storage, and so on. We assume that there is asymmetric competition among the phytoplankton species. That is to say, the large species suffers little competitive effects from small species, while the small species stands more competitive effect from large species. Hence, we adopt $\alpha(x-y)$ to describe the competitive relationship, which represents the effect of cell size $x$ on cell size $y$. When $x = y$, define $\alpha(0)$ as the intraspecific competitive effect. Moreover, $\alpha(x-y)$ is a decreasing function of the difference in the cell sizes of competitive phytoplankton species, which is

Assume that the resident population can be invaded by the mutant population, and the mutant population competes for nutrient with the resident population. If the mutant phytoplankton population with trait $y$ appears in the resident phytoplankton population with trait $x$, then we present the interactions between the resident population and the mutant population as follows

where $P_{y}$ is the density of mutant phytoplankton with trait $y$, $\alpha(x-y)$ is the competitive effect of the mutant population $P_y$ on the resident population $P_x$, $\alpha(y-x)$ is the competitive effect of the resident population $P_x$ on the mutant population $P_y$.

2.2.   Ecological dynamics

In this subsection, we shall study some fundamental dynamics of (2.4) including the non-negativity, boundedness, persistence, and the existence of the coexistence steady state. For simplicity, define $\mathbb{R}^3_+: = \{(N, P_x, P_y) \lvert N \geq 0, P_x \geq 0, P_y \geq 0\}$. We begin with the positive invariance and boundedness of system (2.4).

Theorem 2.1. The system (2.4) is positive invariant and dissipative in $\mathbb{R}^3_+$. More specifically, the system has the following property

Notes. Theorem 2.1 indicates that the full model (2.4) is biologically well defined, that is, any solutions with non-negative initial values are non-negative and bounded for all the future time.

Define

According to Theorem 2.1, we find the compact set $\Omega$ is a positively invariant set and global attracting set for (2.4).

Before exploring the dynamics of the full model (2.4), we investigate the global dynamics of the nutrient-phytoplankton model (2.1), containing one phytoplankton species with fixed cell size $x$. A threshold value for the trait-based nutrient-phytoplankton model (2.1) is introduced, which is the ecological reproductive index

$R^x_0$ represents the average amount of newborn phytoplankton produced by one unit of population $P_x$ during the life expectancy ^[27]. It is a key indicator to evaluate whether phytoplankton species can survive in a given environment or not. Next, we apply the ecological reproductive index of phytoplankton $R^x_0$ to characterize the global dynamics of (2.1).

It is not difficult to find that (2.1) may have two equilibria. Note that the phytoplankton-free equilibrium $E_0(N_{in}, 0)$ always exists, while the positive equilibrium $E^*(N_x^*, P_x^*)$ exists if and only if there exist positive constants $N_x^*$ and $P_x^*$ satisfying the following algebra equations

Theorem 2.2. If $R^x_0 < 1$, then system (2.1) has a unique equilibrium $E_0$ and it is globally asymptotically stable; if $R^x_0 > 1$, then $E_0$ is unstable and there exists a unique positive equilibrium $E^*$ being globally asymptotically stable.

Notes. By the dynamical systems theory, the equilibrium being globally asymptotically stable means an orbit with an arbitrary initial point in the domain will tend toward the equilibrium point. From Theorem 2.2, we find that phytoplankton-free equilibrium $E_0(N_{in}, 0)$ is globally asymptotically stable if the growth rate of the phytoplankton is smaller than its loss rate due to mortality, sinking, and dilution, i.e. $\mu(x)g(N_{in}) < m+(s(x)+v)/z_m$, which means the phytoplankton cannot survive in the environment. The positive equilibrium $E^*(N_x^*, P_x^*)$ to be globally asymptotically stable requires that the growth rate of phytoplankton is larger than its loss rate, i.e. $\mu(x)g(N_{in}) > m+(s(x)+v)/z_m$, then the phytoplankton species can reproduce and exist in the water column.

Next, we explore the persistence and coexistence steady state of the full model (2.4), which reduces to the single-species model (2.1) when $P_x = 0$ or $P_y = 0$. It is not hard to find that system (2.4) has three boundary equilibria

where $(N^*_i, P^*_i), \; i = x, y$ is the interior equilibrium of the single-species model (2.1).

Notice that $R^x_0 = \mu(x)g(N_{in})/\left(m+(s(x)+v)/z_m\right)$. For convenience, we define the following ecological reproductive indexes

We claim that $R^i_0 > R^i_1$ for $i = x, y$. It can be obtained from (2.6). Multiplying the second equation of (2.6) by $P_x^*Q(x)$ and adding the first equation of (2.6), we have

It is easy to find that the right-hand side of the above equation is greater than 0, thus we have $N_{in} > N_x^*$. Since $g(N)$ is a increasing function of $N$, then $g(N_{in}) > g(N_x^*)$. It follows that $R^y_0 > R^y_1$. Similarly, we can get $R^x_0 > R^x_1$.

We have the following theorem about the persistence and permanence of the full model (2.4).

Theorem 2.3. The nutrient $N$ in the full model (2.4) is persistent. In addition, population $P_y$ is persistent if

● $R_0^x > 1$ and $R_1^y > 1$; or

● $R_0^x < 1$ and $R_0^y > 1$.

Similarly, we have following conditions when population $P_x$ is persistent

● $R_0^y > 1$ and $R_1^x > 1$; or

● $R_0^y < 1$ and $R_0^x > 1$.

Moreover, if $\min\{R^x_1, R^y_1\} > 1$, then system (2.4) is permanent and has at least one positive steady state solution $\widetilde{E}^*(\widetilde{N}^*, \widetilde{P}^*_x, \widetilde{P}^*_y)$.

Notes. A species is persistent which means the species is bounded and is larger than a positive constant for quite a long time. Theorem 2.3 describes the possible competitive outcomes for phytoplankton species $P_x$ and $P_y$. Note that the ecological reproductive indexes of phytoplankton play a decisive role in the extinction or persistence of species. $R_0^i$ describes the reproductive capacity of the single species $P_i$ without competitor, while $R_1^i$ indicates the reproductive capacity of population $P_i$ under the competition of population $P_j$, where $i, j = x, y$. The population $P_y$ persists in the following two scenarios.

$(1)$ If the population $P_x$ does not exist in the environment, i.e. $R_0^x < 1$, then it requires that the growth rate of phytoplankton $P_y$ is larger than its loss rate due to mortality, sinking, and dilution, which is $\mu(y)g(N_{in}) > m+(s(y)+v)/z_m$.

$(2)$ If the population $P_x$ exists in the environment, i.e. $R_0^x > 1$, then it requires that the growth rate of phytoplankton $P_y$ is larger than its total loss rate due to mortality, sinking, dilution, and interspecific competition with $P_x$, which is $\mu(y)g(N_x^*) > m+(s(y)+v)/z_m+\alpha(y-x)P_x^*$.

The persistence for population $P_x$ is in the same case. Moreover, if we have $\min\{R^x_1, R^y_1\} > 1$, then both phytoplankton species $P_x$ and $P_y$ can survive and coexist. That is to say, the system (2.4) is permanent.

3.   Evolutionary model and dynamics

3.1.   Evolutionary model

Assume that the phytoplankton cell size is an evolving trait, and the resident population can be invaded by the mutant population. We apply the adaptive dynamics theory ^[28,29] to track the evolution of cell size. In this framework, evolutionary change is assumed to be much slower than ecological dynamics, thus evolutionary dynamics are considered when the resident population has reached a steady state.

Before introducing the evolutionary model, we first established the fitness function of the mutant population, which plays an important role in determining whether the mutant population can invade or not. The invasion of the mutant population can be reflected by the stability of equilibrium $E_1(N_x^*, P_x^*, 0)$ of the full model (2.4). If $E_1$ is stable, then the mutant population cannot invade the resident population, and they would be extinct eventually. Otherwise, if $E_1$ is unstable, then it implies that the mutant population can invade the resident population, and the mutant population and the resident population can coexist in the environment. The stability of $E_1$ can be determined by the following Jacobian matrix

Note that the submatrix of $J_{E_1}$ consisting of the first two rows and the first two columns has the same form as $J_{E^*}$, which is given in Section 6. From Theorem 2.2, we know that the eigenvalues of this submatrix have negative real parts provided by $R_0^x > 1$. Hence, the stability of $E_1$ is determined by the third eigenvalue of $J_{E_1}$, which is

If $\lambda$ is negative, then $E_1$ is asymptotically stable, which implies that the mutant population is extinct. If $\lambda$ is positive, then the mutant population can persist in the environment.

Based on the facts mentioned above, we obtain that the sign of $\lambda$ plays a critical role in determining the fate of mutant phytoplankton species. Therefore, we define the fitness function for the mutant population $P_y$ as follows

where $f_{x}(y)$ is the long-term exponential growth rate of $P_y$. It reveals the destiny of the mutant population. If $f_{x}(y) > 0$, then the mutant population can invade and appear in the resident population. Otherwise, the mutant population cannot invade the resident population and become extinct.

In the adaptive dynamics framework, another important function is the local fitness gradient, which decides the direction of evolution. The local fitness gradient $D(x)$ is defined by

where $\alpha'_y(0) = \left.\frac{\partial \alpha(y-x)}{\partial y}\right|_{y = x}$. We assume that the variation of the cell size is random and small such that the trait $y$ is very close to $x$. Hence, by linear approximation of the fitness function, we yield

Substituting $y = x$ into (3.1), we have $f_x(x) = 0$ by the second equation of (2.6). Hence, the sign of $f_x(y)$ is determined by the sign of $D(x)$ and the difference between $x$ and $y$. If $D(x) > 0$, then the mutant with trait larger than $x$ can replace the resident phytoplankton; if $D(x) < 0$, then the mutant with trait smaller than $x$ can take over the resident population.

Next, we formulate the evolutionary model for the phytoplankton species and explore the adaptive dynamics of the cell size $x$. According to the adaptive dynamics theory, the mutation rate of a phenotypic trait is proportional to the probability that a newborn individual in the resident is a mutant, the variance of mutation distribution, the resident population density, and the local fitness gradient. Hence, the evolutionary model for the cell size $x$ is described as

where $\theta$ is the probability that a newborn individual in the resident is a mutant, $\sigma^2$ is the variance of mutation distribution, $P_x^*(x)$ is the density of resident phytoplankton at steady state, and $D(x)$ is the local fitness gradient.

It is interesting to find that the adaptive model (3.3), a one-dimensional differential equation, can describe competitive effects among phytoplankton species of all cell sizes and reveal the evolutionary dynamics of multi-species systems.

3.2.   Evolutionary dynamics

The evolutionarily singular strategy $x^*$ is one of the most important studies of evolutionary dynamics, which is defined as $D(x^*) = 0$. In our case, the condition for evolutionarily singular strategy $x^*$ is given by

By calculating the equation (3.4), we can obtain the value of the evolutionarily singular strategy $x^*$.

Next, we focus on what the phytoplankton species go through after reaching the evolutionarily singular strategy $x^*$. The strategy may be evolutionarily stable, which means that the new mutant species can no longer invade the resident species; or undergo diversification, which implies that the mutant species can coexist with the resident species in the environment. This depends on the evolutionary stability or convergence stability of the evolutionarily singular strategy, which is defined as follows.

Evolutionarily stable strategy. The evolutionarily singular strategy is called an evolutionarily stable strategy if there are no mutants nearby $x^*$ can invade, which is determined by the following inequality

Convergence stable strategy. The evolutionarily singular strategy is known as convergence stable strategy if mutants nearby $x^*$ can be invaded by trait being even more close to $x^*$, which can be derived from the following inequality

According to the above definitions, we can further classify the evolutionarily singular strategy.

Continuously stable strategy. If the singular strategy is both evolutionarily stable and convergence stable, then it is said to be a continuously stable strategy (CSS). The fitness function has a maximum at CSS. It implies the endpoint of evolution, and the specific strategy prevents invasion by others.

Evolutionary branching point. If the singular strategy is convergence stable but lacks evolutionary stability, then it is an evolutionary branching point. The evolutionary branching point can lead to speciation, and the mutant population can coexist with the resident population.

Repellor. If the singular strategy lacks convergence stability, then we refer to the strategy as a repellor.

In our case, we have

where $\alpha''_{yy}(0) = \left.\frac{\partial^2 \alpha (y-x)}{\partial y^2}\right|_{y = x = x^*}$, $g'(N_x^*(x^*)) = \left.\frac{ { {\mathrm{d}} g(N_x^*)}}{ { {\mathrm{d}} N_x^*}}\frac{ { {\mathrm{d}} N_x^*}}{ { {\mathrm{d}} x}}\right|_{x = x^*}$, $P_x^{*'}(x^*) = \left.\frac{ { {\mathrm{d}} P_x^*}}{ { {\mathrm{d}} x}}\right|_{x = x^*}$.

Theorem 3.1. Assume that $R^x_0 > 1$. For the evolutionarily singular strategy $x^*$, one has

(i) if both (3.7) and (3.8) are negative, then $x^*$ is a CSS;

(ii) if (3.7) is positive and (3.8) is negative, then $x^*$ is an evolutionary branching point;

(iii) if (3.8) is positive, then $x^*$ is a repellor.

Notes. Theorem 3.1 provides the sufficient conditions for the evolutionarily singular strategy to be a CSS or an evolutionary branching point or a repellor, which is determined by the growth rate $\mu(x)$, sinking rate $s(x)$, and the competitive effect $\alpha(x-y)$. The biological meaning of CSS is the endpoint of the evolution process, which implies the phytoplankton with this strategy cannot be invaded by other phytoplankton species. CSS represents the most competitive cell size of phytoplankton species in the environment. At the evolutionary branching point, the cell size has the potential for adaptive diversification, which causes the phytoplankton to split into two species with different cell sizes. A repellor lacks convergence stability and thus is unattainable in the evolution process.

Many theoretical ecologists have applied the eco-evolutionary model, including both ecological and evolutionary changes, to study the ecological and evolutionary dynamics of population. In this paper, an eco-evolutionary model can be established based on the ecological model (2.1) and the evolutionary model (3.3), which is given by

Note that the eco-evolutionary model (3.9) exists a positive equilibrium $\hat{E}(N_x^*, P_x^*, x^*)$ when $R_0^x > 1$. The stability of the equilibrium $\hat{E}(N_x^*, P_x^*, x^*)$ can be ensured by the dynamics of evolutionary models. We have the following proposition.

Proposition 3.1. The following two expressions are equivalent.

(i) The equilibrium $\hat{E}(N_x^*, P_x^*, x^*)$ of the eco-evolutionary model (3.9) is locally asymptotically stable.

(ii) The evolutionarily singular strategy $x^*$ of (3.3) is convergence stable.

Notes. Proposition 3.1 suggests that the equilibrium of (3.3) and the equilibrium of the eco-evolutionary model (3.9) have the same stability. Therefore, to a certain extent, these two systems can be referred to be equivalent. In this paper, we adopt a more concise approach to studying the evolution of phytoplankton phenotypic traits by separating the ecological time scale from the evolutionary time scale. This is because the turnover of phytoplankton is very fast, a few hours or days, while evolution is a long-time behavior. Thus, we assume that the ecological and evolutionary time scales can be separable.

4.   Applications

In this section, we illustrate numerically the effects of environmental drivers on the ecological and evolutionary dynamics of phytoplankton. First, we need to determine the mathematical formula for size-dependent ecological parameters including the growth rate $\mu(x)$, nutrient quota $Q(x)$, sinking rate $s(x)$, and competitive effect $\alpha(x-y)$. Many researchers have studied the correlations between phytoplankton cell size and the ecological processes and apply functions of cell size $x$ to describe these relationships.

Considering the growth rate $\mu(x)$, Pu et al. ^[20] adopted the following formula to describe the size-dependent relationship, which reads

where growth-dependent coefficients $c_1, c_2, c_3$ are selected to satisfy the condition (2.2).

Regarding $Q(x)$, Litchman et al. ^[11] employed an exponential function to explain the relationship between the nutrient quota and phytoplankton cell size, which is

where $\beta$ is the nutrient quota coefficient and $\gamma$ is the nutrient quota exponent.

For spherical phytoplankton cells, the sinking rate $s(x)$ can be given by the Stokes equation ^[10]. It assumes that the sinking rate is proportional to the square of the spherical diameter, which reads

where $s_0$ is the sinking rate coefficient, and it is related to the density and viscosity coefficient of the water.

The asymmetric competition between phytoplankton species is characterized by the following concave-convex function

where $\alpha_0$ and $k$ are intraspecific and interspecific competition coefficients, $w$ is the strength of competition.

The default values for the parameters are listed in Table 1.

4.1.   Numerical investigation of the ecological dynamics

Based on the analytical results, we know that the persistence and the existence of coexistence steady state of (2.4) is highly dependent on nutrient-related parameters. Thus, some numerical results are carried out to study their effects on phytoplankton dynamics.

In the following numerical simulations, we apply equations (4.1)-(4.4) to describe the size-dependent relationships. We next explore the effects of input nitrogen concentration $N_{in}$ and the vertical mixing rate $v$ on the distribution of cell size and persistence of phytoplankton species.

The cell size distribution of viable phytoplankton varies with environmental conditions. From Theorem 2.2, the survival of a single phytoplankton population is determined by the ecological reproductive index $R^x_0$. By calculating $R^x_0$, we can obtain the cell size range of viable phytoplankton. As shown in Figure 2, the cell size range of viable phytoplankton becomes larger as the input nitrogen concentration $N_{in}$ increases, while the cell size range of viable phytoplankton becomes smaller as the vertical mixing rate $v$ rises. This is due to the fact that the survival of phytoplankton is determined by its growth rate and loss rate. As $N_{in}$ increases, it provides enough nutrients to increase the growth rate of phytoplankton, while an increase in the vertical mixing rate $v$ leads to a large loss rate of phytoplankton. Thus, phytoplankton species with a wider range of cell size can survive as the input nitrogen concentration $N_{in}$ increases, whereas only phytoplankton with medium cell sizes can persist with an increase in the vertical mixing rate $v$.

We explore the effects of nutrient input on the persistence of two competitive phytoplankton species. From Theorem 2.3, the persistence of phytoplankton species $P_x$ and $P_y$ is determined by the values of ecological reproductive indexes $R^x_0$, $R^x_1$, $R^y_0$ and $R^y_1$, which depend on the parameters related to nutrient input. Figure 3 shows the effects of the input nitrogen concentration $N_{in}$ and the vertical mixing rate $v$ on the ecological reproductive indexes. Here we set the cell sizes of the two phytoplankton species to $x = 2$ and $y = 3$, respectively.

The effect of nitrogen input concentration $N_{in}$ is shown in Figure 3 (a) when the vertical mixing rate $v = 3$. $R^x_0$, $R^x_1$, and $R^y_0$ increase rapidly to 1, and then $R^y_1$ slowly reaches 1. After that, $R^x_1$ drops slightly less than 1. The corresponding dynamics change from the extinction of two species to the persistence of $P_x$ to the coexistence of $P_x$ and $P_y$ to the persistence of $P_y$ (see Figure 4 (a)). For the nitrogen input concentration $N_{in}$, we have the following results.

● When the input nutrient concentration is very low, it is not sufficient to support the survival of phytoplankton, and thus both species become extinct.

● When the nutrient increases to a certain level such that $R^x_1 > 1$, the small phytoplankton species $P_x$ can survive in the environment.

● When the input nutrient concentration becomes adequate, both phytoplankton species $P_x$ and $P_y$ can coexist in the water column.

● As the nutrient concentration continues to increase, the large phytoplankton species $P_y$ dominates the competition and persists in the environment.

In general, small phytoplankton dominates the water column when the nutrient is limited due to its advantages in nutrient absorption and uptake, while large phytoplankton begins to dominate the competition when the nutrient is adequate.

We then investigate the effect of the vertical mixing rate $v$ on the persistence of phytoplankton when input nitrogen concentration is fixed $N_{in} = 9$. Figure 3 (b) shows that the vertical mixing rate $v$ changes $R_1^x$ and $R_1^y$ a lot, as does the persistence of phytoplankton species $P_x$ and $P_y$ (see Figure 4 (b)). For the vertical mixing rate $v$, we have the following results.

● When $v$ is relatively low ($v\in(0.17, 0.83)$) or relatively high ($v\in(4.47, 18.66)$), all the ecological reproductive indexes are all greater than 1, which means two phytoplankton species can coexist.

● When $v$ is extremely low ($v < 0.83$) or high ($v > 18.66$), the small phytoplankton $P_x$ dominates the environment and the large phytoplankton $P_y$ cannot survive.

● When $v$ takes an intermediate value ($v\in(0.83, 4.47)$), the large phytoplankton $P_y$ dominates the environment and the small phytoplankton $P_x$ cannot survive.

It implies that a too low or too fast mixing rate is not conducive to the survival of large phytoplankton. This may be due to the fact that when the mixing rate is small, the nutrient input is insufficient to support the survival of large phytoplankton. When the mixing rate is fast, large phytoplankton would be rapidly lost from the mixed layer. When the nutrients are sufficient and the mixing rate is not too high, large individual phytoplankton begins to dominate in the water column.

The joint effects of input nitrogen concentration $N_{in}$ and the vertical mixing rate $v$ on phytoplankton persistence are illustrated in Figure 5. There are four competitive outcomes of two phytoplankton species, including extinction, persistence of $P_x$, persistence of $P_x$ and coexistence of $P_x$ and $P_y$. Figure 5 shows the extinction or persistence of phytoplankton species with $x = 2$ and $y = 3$ by varying both $N_{in}$ and $v$ from 0 to 20, where the black region marked with $\bar{E_0}$ represents extinction; the green region marked with $E_1$ represents persistence of $P_x$; the blue region marked with $E_2$ represents persistence of $P_y$; and the red region marked with $\widetilde{E}^*$ represents coexistence of $P_x$ and $P_y$. First, we look at the effect of $N_{in}$.

● When $N_{in}$ is small, the nutrient is insufficient to support phytoplankton growth.

● As $N_{in}$ increases, small individuals gradually dominate in the water column.

● When $N_{in}$ is large enough, large individuals can coexist with small individuals.

● When $N_{in}$ increases further, large individuals begin to dominate in the environment.

This suggests that small individuals dominate in the environment when nutrient is limited, while large individuals gradually begin to dominate the competition under adequate nutrient supply.

The vertical mixing rate $v$ has an interesting and complex effect on the dynamics of phytoplankton, and it depends on the input nitrogen concentration $N_{in}$.

● When $N_{in}$ is small, a moderate $v$ allows small phytoplankton to survive, while too high or too low mixing rate leads to phytoplankton extinction.

● When $N_{in}$ takes intermediate value, a moderate $v$ favors the coexistence of both populations, while only small individuals $P_x$ can survive under too high or too low mixing rate.

● When $N_{in}$ is very abundant, a moderate $v$ is beneficial to the survival of large phytoplankton $P_y$, and too fast or too low $v$ enables the coexistence of large phytoplankton and small phytoplankton.

In brief, the input nitrogen concentration $N_{in}$ and the vertical mixing rate $v$ have significant effects on the phytoplankton dynamics.

In addition to environmental influences on the persistence of phytoplankton, cell sizes can also affect the survival of phytoplankton. Figure 6 demonstrates the effect of cell sizes on the persistence of phytoplankton with $N_{in} = 10$, $v = 3$. The numerical result shows symmetry. When one of the phytoplankton cell sizes is very small, the larger individual prevails; when one of the phytoplankton cell sizes is very large, the smaller phytoplankton individual prevails. In addition, phytoplankton coexistence is mostly distributed in areas with large differences in cell sizes between the two species. This validates the ecological niche theory, which represents the ecological niches of species incline to differ in some respects, as interspecific competition minimizes their overlap.

4.2.   Numerical investigation of the evolutionary dynamics

In this subsection, we numerically demonstrate the effects of nutrient input on the evolutionary dynamics of marine phytoplankton, such as the variation of cell size and trait diversity. We concentrate on how the optimal cell size $x^*$ of phytoplankton varies with the input nitrogen concentration $N_{in}$ and the vertical mixing rate $v$.

From Figure 7 (a), we find that the evolutionarily singular strategy $x^*$ increases monotonically with the input nitrogen concentration $N_{in}$. It means that the cell size of phytoplankton increases with the nutrient concentration of deep ocean brought about by the upwelling and water mixing. In other words, the cell sizes of phytoplankton species in eutrophic waters are generally larger than that in oligotrophic waters. It is consistent with the ecological phenomena and observations in the ocean ^[8]. It is observed that the marine ecosystem in nutrient-rich areas is predominated by large phytoplankton species, while small phytoplankton species tend to dominate in nutrient-depleted waters.

The input nutrient concentration not only changes the optimal cell size of phytoplankton but also alters the stabilities of the evolutionarily singular strategy. In other words, it can affect the trait diversity of phytoplankton species. As shown in Figure 7 (a), the black line indicates that the evolutionarily singular strategy $x^*$ is a CSS and the red line indicates that the evolutionarily singular strategy $x^*$ is an evolutionary branching point. It implies that there exists a threshold $N_{in}^*$, when the input nitrogen concentration $N_{in}$ is less than $N_{in}^*$, $x^*$ is a CSS; when the input nitrogen concentration $N_{in}$ is large than $N_{in}^*$, $x^*$ is an evolutionary branching point. The pairwise invasibility plot is a useful tool to numerically illustrate the stability of evolutionarily singular strategy, as shown in Figure 8. Fix the vertical mixing rate $v = 3$, we find that $x_1^*$ is a CSS when $N_{in} = 8$ ($N_{in} < N_{in}^*$) and $x_2^*$ is an evolutionary branching point when $N_{in} = 15$ ($N_{in} > N_{in}^*$). Hence, abundant nutrient supply is in favor of maintaining the high trait diversity of phytoplankton. This fact is in accord with the previous study, where it is found that there is a positive relationship between the nutrient supply and phytoplankton community structure ^[30]. The increasing supply allows the survival of phytoplankton having a large cell size and does not decrease the phytoplankton with a small cell size. Thus, it leads to the coexistence of small phytoplankton and large phytoplankton and promotes the trait diversity of phytoplankton species.

There is a nonmonotonic relationship between the evolutionarily singular strategy $x^*$ and the vertical mixing rate $v$ as illustrated in Figure 7 (b). As the vertical mixing rate increases, the cell size presents a unimodal trend with respect to $v$. When the vertical mixing rate is at a low level, the phytoplankton can evolve into large individuals with the increase in the mixing rate. This may be due to the fact that an increase in the mixing rate of the water brings more nutrients to support large phytoplankton growth in the upper layer. On the contrary, a drop in the vertical mixing rate can result in a decline in the cell size of phytoplankton. This explains the experimental results, that is, the cell size of phytoplankton decreases when the nutrient supply is becoming limited ^[14]. However, when the mixing rate is at a high level, the size of the phytoplankton starts to decrease as the mixing rate continues to increase. Smaller individuals have a greater survival advantage than larger individuals when the mixing rate is high. Thus, it means that a moderate mixing rate promotes a competitive advantage for large cells, while a low or high mixing rate reduces the optimal cell size of phytoplankton.

Comparing Figure 9 (a), (b) and (c), one observes that $x_1^*$ is a CSS when $v = 0.5$, $x_2^*$ is an evolutionary branching point when $v = 2.5$, and $x_3^*$ is a CSS when $v = 10$. It means that a moderate vertical mixing rate can promote the trait diversity of phytoplankton species, while a low or high vertical mixing rate is beneficial to maintain the evolutionary stability of the evolutionarily singular strategy. When the mixing rate is low, large phytoplankton cells cannot survive due to too poor nutrient conditions; when the mixing rate is high, large cells are removed from the mixing layer by the vertical mixing. Thus, small phytoplankton is more competitive at a low mixing rate or high mixing rate. This finding is in accordance with the ecological observation that moderate water mixing can enhance the diversity of phytoplankton ^[31]. Hence, the diversity of phytoplankton species is high in areas where the water mixing is intermediate, but only a few phytoplankton species dominate in the oceans where the water exchange is slow or intense.

5.   Discussion and conclusion

In this paper, we develop a dynamic model based on nutrient competition to explore the effects of nutrient supply on phytoplankton cell size evolution. The analysis shows that the nutrient supply, including the input nitrogen concentration and the vertical mixing rate, play major roles in the population and evolutionary dynamics of phytoplankton. We obtain that increasing input nitrogen concentration leads to larger cell size and a high level of diversity, while there is a unimodal relationship between cell size and vertical mixing rate, with a moderate vertical mixing rate promoting large cell size and trait diversity. It is noteworthy that when the vertical mixing rate is at a low level, the optimal cell size decreases sharply with decreasing water mixing. Thus, a decrease in the input nutrient concentration or a slowdown in the vertical mixing rate can result in the evolution of phytoplankton toward smaller cell size and less trait diversity.

Our results are consistent with ecological phenomena and observations ^[17]. It is well known that, in nutrient-enriched waters such as Peruvian upwelling, the phytoplankton population has the following features, such as high biomass and large cell size. While in nutrient-poor waters, for example, the ocean areas near the subarctic North Pacific, there are low phytoplankton diversities, and the small phytoplankton cells tend to dominate there. The reason for this phenomenon is that, in the case of poor-nutrient status, small cells have an advantage in growth rate, while large individuals have high nutrient requirements, and insufficient nutrients cannot support the survival of large cells. Therefore, natural selection tends to select small cells. As the nutrient supply increases, the nutrient in the mixed layer becomes abundant, and it no longer restricts the growth of large cells. Thus, large cells can coexist with small cells, so the trait diversity in nutrient-rich areas is high.

We find that the evolutionary model based on the single-species model admits evolutionary branching points when some parameter conditions are satisfied. The biological meaning of the evolutionary branching point is trait diversity, which means that more than one species can survive in the environment. In this paper, the size-dependent competition between phytoplankton species is taken into consideration in our model. Since the phytoplankton species compete for light and nutrients in the mixed layer, the competitive effect can be size-dependent. The results show that the competition can promote cell size diversity and the coexistence of marine phytoplankton with multiple cell sizes. Our current findings extend prior work on the evolutionary dynamics of phytoplankton.

The framework proposed in this paper can be used to predict the cell size evolution of marine phytoplankton under climate change. As we all know, climate change has significant impacts on the marine environment, especially limiting the nutrient supply from the deep ocean. The variation in the marine environment can be reflected by related parameters in our model, for example, a reduction in nutrient availability can be described by a decrease in input nutrient concentration or a slowdown in the vertical mixing rate. From the theoretical and numerical analysis, we obtain that decreasing input nitrogen concentration and slower water mixing may lead to phytoplankton evolving to smaller cell sizes and less species diversity. Therefore, we predict that marine phytoplankton will evolve into smaller cells and have less trait diversity in the future.

However, some limitations are worth noting. We only consider the interactions between nitrogen and phytoplankton and ignore the influence of higher trophic levels. For simplicity, the impact of predation by higher trophic levels (e.g., zooplankton) is included in the mortality rate of phytoplankton. However, due to the cascade effect to or from high trophic levels, the change in biomass and cell size of phytoplankton can have a significant impact on the zooplankton. Therefore, future research should introduce the zooplankton population into the nutrient-phytoplankton model and explore the coevolution of phytoplankton and zooplankton. In addition, we only consider the effect of the variation of nutrient supply in our model. In fact, temperature also has an important effect on phytoplankton evolution, and rising temperature can indirectly or directly affect the evolutionary behavior of phytoplankton cell size. It has been shown that increasing temperature will indirectly change the nutrient supply, while the rising temperature will directly change the growth rate of phytoplankton, making phytoplankton with faster generation times. Therefore, it will be of great significance to study the combined effects between altered nutrient supply and changed temperature on the size evolution of phytoplankton.

6.   Some mathematical proofs

Proof of Theorem 2.1.

Proof. For any given initial value $\left(N(t_0), P_x(t_0), P_y(t_0)\right)\in \mathbb{R}^3_+$, we obtain

for all $t \in [t_0, +\infty)$.

Similarly, we get

for all $t \in [t_0, +\infty)$.

Next, we prove $N(t)$ is non-negative when $t \geq t_0$. If $N(t_0) = 0$, then $N(t) > 0$ based on the change rate of nitrogen ${\mathrm{d}} N/ {\mathrm{d}} t > 0$. Next, we consider $N(t_0) > 0$. If the statement is not true, then there would exist a $t_1 > t_0$ such that $N(t_1) = 0$ and $N(t) \geq 0$ for $t \in(t_0, t_1)$. Given that $N(t_0) > 0$, it implies ${\mathrm{d}} N / {\mathrm{d}} t \leq 0$ when $t = t_1$. From the first equation of (2.4), we derive

which is a contradiction. Therefore, $N(t) \geq 0$ for $t\geq t_0$.

We are now in the position to prove (2.4) is dissipative. Let

From (2.4), we have

It is easy to find that $\lim_{t\to\infty}V(t)\leq N_{in}$. Hence, the system (2.4) is dissipative. The proof is complete.

Proof of Theorem 2.2.

Proof. First we consider the case when $R^x_0 < 1$. The Jacobian matrix of (2.1) evaluated at $E_0$ is

The two eigenvalues of $J_{E_0}$ are

It is easy to find that $\lambda_{11} < 0$, and $\lambda_{12} < 0$ holds true if and only if $R^x_0 < 1$. Thus, when $R^x_0 < 1$, all the eigenvalues of $J_{E_0}$ have negative real parts, so $E_0$ is locally asymptotically stable. In addition, $E_0$ is the unique equilibrium of (2.1) when $R^x_0 < 1$. From Theorem 2.1, we know that (2.1) is dissipative. By Poincar$\mathrm{\acute{e}}$-Bendixson criterion, (2.1) has no closed orbit. Thus, $E_0$ is globally asymptotically stable. When $R^x_0 > 1$, we have $\lambda_{12} > 0$. It implies that one eigenvalue of $J_{E_0}$ has positive real part, so $E_0$ is unstable.

Next, we explore the existence and stability of $E^*$. Rearranging the algebra equations of (2.6), we can obtain that $E^*$ exists if the two curves defined by $N_1(P)$ and $N_2(P)$ have intersection in $\mathbb{R}^2_+$, where the two curves are described by

From (6.1), we find $N_1(P)$ is a downward-opening parabola and its axis of symmetry is $P = -((1-r)m z_m+s(x)+v)/(2z_m \alpha(0)) < 0$. Moreover, $N_1(0) = N_{in} > 0$. Hence, $N_1(P)$ decreases with $P$ in $\mathbb{R}^2_+$. In addition, we know that $N_2(P)$ is an increasing hyperbola, and its vertical asymptote is $P = (\mu(x)z_m-mz_m-s(x)-v)/(\alpha(0)z_m) > 0$. Also, $\lim _{P\to\infty}N_2(P) = -K < 0$. It follows that $N_2(P)$ increases with $P$ in $\mathbb{R}^2_+$. Thus, $N_1(P)$ and $N_2(P)$ have a unique interior intersection in $\mathbb{R}^2_+$ if and only if $N_1(0) > N_2(0) > 0$, which is equivalent to $R^x_0 > 1$. Therefore, there exists a unique positive equilibrium $E^*$ of (2.1) provided that $R^x_0 > 1$.

The stability of $E^*$ can be derived from the Jacobian matrix of (2.1) evaluated at $E^*$

The trace of $J_{E^*}$ is

The determinant of $J_{E^*}$ is given by

It follows that all the eigenvalues of $J_{E^*}$ have negative real parts, so $E^*$ is locally asymptotically stable. The global stability of $E^*$ can be derived from the Bendixson-Dulac criterion. Let

Choosing a Dulac function $B(N, P_x) = 1/P_x$, we yield

By the Bendixson-Dulac criterion, (2.1) has no closed orbit in the plane. Therefore, $E^*$ is globally asymptotically stable. The proof is complete.

Proof of Theorem 2.3.

Proof. According to the proof of Theorem 2.1, we have for any $\epsilon > 0$

for time $t$ large enough. Therefore, for time $t$ large time, we have

Choose $\epsilon = v N_{in}/(2z_m)$, then we have

which implies that

Therefore, the nutrient $N$ is persistent.

Notice that the $\omega$-limit set of the $N-P_x$ subsystem of the full model (2.4) when $P_y = 0$ is $E_0(N_{in}, 0)$ if $R^x_0 < 1$, while the $\omega$-limit set is $E^*(N^*_x, P^*_x)$ if $R^x_0 > 1$. Thus according to Theorem 2.5 ^[32], we can conclude species $P_y$ is persistent in $\mathbb{R}^3_+$ if one of the following cases hold.

$(1)$ In the case that $R^x_0 < 1$, we have the subsystem $N$-$P_x$ converges to $E_0$ globally from Theorem 2.2. Thus,

$(2)$ In the case that $R^x_0 > 1$, we have the subsystem $N$-$P_x$ converges to $E^*$ globally from Theorem 2.2. Thus,

Similarly, we are able to show that species $P_x^*$ is persistent if $R_0^y > 1$ and $R_1^x > 1$; or $R_0^y < 1$ and $R_0^x > 1$.

Note that $R_0^i > R_1^i$ for $i = x, y$, therefore, if both species $P_x$ and $P_y$ are persistent, then it requires that $\min\{R^x_1, R^y_1\} > 1$.

Because $N$ is always persistent. Therefore, the full model (2.4) is permanent if $\min\{R^x_1, R^y_1\} > 1$. According to the fix point theorem, we can conclude that the full system has at least one interior equilibrium point under the condition that $\min\{R^x_1, R^y_1\} > 1$.

Proof of Proposition 3.1.

Proof. The stability of $\hat{E}$ of (3.9) can be derived from the Jacobian matrix

The equilibrium $\hat{E}(N_x^*, P_x^*, x^*)$ is locally asymptotically stable provided that all the eigenvalues of $J_{\hat{E}}$ have negative real parts. It is equivalent to the eigenvalues of the submatrix of $J_{\hat{E}}$, which consists of the first two rows and the first two columns, have negative real parts and the third eigenvalue $\frac{1}{2} \theta \sigma^2P_x^*D'(x^*)$ is negative. From Theorem 2.2, we know that the eigenvalues of the submatrix have negative real parts. Thus, the equilibrium $\hat{E}(N_x^*, P_x^*, x^*)$ is locally asymptotically stable if and only if $D'(x^*) < 0$. This completes the proof.

Acknowledgments

This work is supported by National Natural Science Foundation of China (No. 12071068) and China Scholarship Council (No. 201906620063). The work of YK is particially supported by NSF-DMS (Award Number 1716802 & 2052820); and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472).

Conflict of interest

The authors declare there is no conflict of interest.