EFFECT OF TEMPERATURE ON ADAPTIVE EVOLUTION OF PHYTOPLANKTON CELL SIZE∗

We present a simple nutrient–phytoplankton model that incorporates adaptive evolution and allometric relations. This model allows us to examine the patterns and consequences of adaptive changes in the cell size of phytoplankton under the effect of changes in water temperature. A theoretical study reveals that the ecological reproductive index can be used to characterize the evolutionary dynamics of the nutrient–phytoplankton model. Numerical analysis suggests that phytoplankton should evolve toward the small sizes typical of picophytoplankton as the water temperature increases. This study provides a framework for studying the adaptive evolution of phytoplankton cell size in water ecosystem.


Introduction
Phytoplankton is a polyphyletic group of single-celled primary producers that are ubiquitous in aquatic ecosystems. A diverse phytoplankton assemblage composed of species of different sizes, has existed for hundreds of millions of years, the result of long-term evolution [16]. The cell size of phytoplankton ranges over several different orders of magnitude. Body size is one of the most fundamental traits of organisms, affecting almost all aspects of their physiology and ecology. Phytoplankton cell size affects not only physiological rates but also community structure and ecological function [11,23]. Thus, it is thought to be a promising ecophysiological trait for modeling and tracking population dynamics in response to biotic and abiotic factors [16].
Water temperature is a key environmental factor for the aquatic ecosystem, and has essential impacts on the nutrient uptake, growth, and death of phytoplankton. Therefore, water temperature can alter the productivity and composition of phytoplankton communities, thereby affecting global biogeochemical cycles [21,28,29].
Recent studies [22] show that temperature alone is able to explain 73% of the variance in the relative contribution of small cells to total phytoplankton biomass, regardless of differences in trophic status or inorganic nutrient loading. There also exists an interesting phenomenon whereby small phytoplankton species dominate the equatorial and subtropical oceans while larger species are more abundant in subpolar regions [30]. This phenomenon is the result of long-term evolution related to water temperature. This evidence demonstrates that water temperature plays an important role in phytoplankton evolution. Hence, it is important and reasonable to incorporate all possible direct temperature effects on phytoplankton evolution into any mathematical modeling approach. Nevertheless, most mathematical modeling studies on the evolution of phytoplankton have ignored the impact of water temperature [15,16]. Continued and increasing global warming raises two interesting questions: how do the phytoplankton communities respond to global warming, and how does the cell size of phytoplankton evolve under the impact of global warming over the long term?
There have been many conceptual, experimental, and theoretical attempts to predict the evolution of phytoplankton cell size [2,14,15,19]. Adaptive dynamics provide a useful framework for modeling the evolution of quantitative traits. This approach assumes that evolutionary and ecological dynamics occur on different time-scales, separates the two dynamical processes analytically, and draws on the feedback between ecological and evolutionary processes [12,15]. In the adaptive dynamics approach to modeling evolutionary changes in dynamical models of ecological communities, the evolutionary processes are usually governed by two important indexes, i.e., invasion fitness and fitness gradient. Species always evolve to maximize their fecundity in a specific environment. In [5], an ecological reproductive index is defined to measure the fecundity and characterize the growth of phytoplankton, which also allows a comprehensive analysis of the role of temperature on phytoplankton's growth and reproductive characteristics. This suggests another question to be addressed: what role does the ecological reproductive index play in the adaptive dynamics of phytoplankton? Motivated by the above considerations, the principal aims of this study are twofold. The first is to develop a new reproductive index that characterizes the evolutionary dynamics of phytoplankton cell size. The second is to explore the patterns and consequences of the evolutionary trait (i.e., cell size) and investigate the impacts of water temperature on the evolutionary dynamics. In Section 2, we establish a nutrient-phytoplankton model that consider the effect of water temperature on the adaptive evolution of phytoplankton cell size, and then explore the global dynamics of the model with the help of a newly defined ecological reproductive index, R 0 . Section 3 formulates an evolutionary model and theoretically studies the evolutionary behavior of phytoplankton cell size with the help of R 0 and pairwise invisibility plots. Section 4 is devoted to investigating the effects of water temperature on the evolution dynamics of phytoplankton. Finally, Section 5 ends this paper with a discussion of the main results.

Ecological model
In this section, we focus on the dynamical ecological interaction between the phytoplankton and the nutrient (say phosphorus). Phytoplankton cells take a wide variety of different forms, including filamentous shapes and the presence of spines. Without any loss of generality, the volume of differently shaped cells (V cell ) can be converted into the estimated spherical diameter (ESD) via the equation ESD=2(3V cell /4π) 1/3 [24]. To facilitate the following discussion, we use the ESD defined in [24] and applied in [16] to represent the phytoplankton cell size.

Model formulation
Let N (t) be the density of soluble mineral nutrient in units of milligrams of nutrient per stere and let P (t) be the density of phytoplankton at time t in units of cell number per stere. The dynamical interaction between the limiting nutrient and the phytoplankton is described by where D(N in − N ) models the nutrient exchange with an external source; µ opt (x) is the maximum specific growth rate (day −1 ) of phytoplankton, s(x) is the loss rate (day −1 ), Q(x) is the amount of nutrient in an individual cell (mg · cell −1 ), which are all size dependent (i.e., function of ESD x); m is the specific mortality rate (day −1 ) of phytoplankton, which is size-independent; f (N ) represents the nutrient limitation experienced by phytoplankton; h(T ) and k(T ) represent the temperature limitations on the growth rate and loss rate of phytoplankton, respectively. In (2.1), the comprehensive impact of the cell size and water temperature on the growth and loss of phytoplankton is assumed to be of multiplicative type. The maximum specific growth rate of phytoplankton, µ opt (x), does not bear a simple relationship with size over the possible range of phytoplankton cells [16]. Many factors related to cell size impact the loss rate s(x) of phytoplankton, as do the natural mortality, sinking, intraspecific competition, and predation by higher trophic levels (such as zooplankton and fish) [31,32]. Therefore, in different cases, µ opt (x) and s(x) can take different forms.
To characterize the effect of temperature on the adaptive evolution of phytoplankton, we first define an ecological reproductive index of phytoplankton as which is a function of cell size x. R 0 is a key indicator of phytoplankton fecundity, measuring the average newborn phytoplankton produced by one unit of phytoplankton during the phytoplankton life span [5].
Proof. For the positive equilibrium E * (N * , P * ), N * and P * can be determined from the following system of algebraic equations Then, by (2.3), it is easy to show that, if R 0 > 1, where the ecological equilibrium E * (N * (x), P * (x)) is determined by the cell size x and temperature T . Because f (N ) > 0, such an N * (x) is unique. Thus, there is a unique positive equilibrium E * (N * (x), P * (x)) of (2.1).
Consider the Lyapunov function defined by Differentiating V (N, P ) along the solutions of (2.1) and using (2.3), we obtain Obviously, dV /dt = 0 if and only if (N, P ) = (N * , P * ). This implies that, when

Evolutionary model
In this section, we establish a resident-mutant model to explore the effect of temperature on the adaptive evolution of phytoplankton. Consider the body or cell size as the evolutionary trait of phytoplankton, and assume that the rate of change of this trait is slower than the ecological dynamics.

Model formulation
When a mutant phytoplankton population with a slightly different cell size y appears, it is reasonable and realistic to ask which population will become dominant in the system. To derive the invasion fitness for mutant phytoplankton and model the dynamics of the resident-mutant population, we modify the resident phytoplankton model (2.1) to where P m denotes the population density of mutant phytoplankton. Before the occurrence of the mutation, the resident population density is stable at (N * (x), P * (x)) under system (2.1). Just after the occurrence of a rare mutation, however, the resident and mutant populations are close to the ecological equilibrium (N * (x), P * (x), 0) of (3.1). That is, the stability of (N * (x), P * (x), 0) is crucial in determining whether the mutant phytoplankton can invade. The stability of this ecological equilibrium is now analyzed. The Jacobian matrix of (3.1) at (N * (x), P * (x), 0) takes the form Note that J sub is the Jacobian matrix of model (2.1) at E * (N * (x), P * (x)). All the eigenvalues of J sub have negative real parts when R 0 > 1. Note that the third eigenvalue of J(N * (x), P * (x), 0) is which denotes the per capita growth rate of a very scarce mutant phytoplankton with cell size y. To facilitate the following discussion, we define this quantity as F (y, x). If F (y, x) < 0, then (N * (x), P * (x), 0) is asymptotically stable, which means that the mutant phytoplankton population fails to invade and tends to become extinct; if F (y, x) > 0, then (N * (x), P * (x), 0) is unstable, which means that the abundance of the mutant phytoplankton will initially increase, that is, the mutant phytoplankton population can invade. Because of its important role in determining the fate of the mutant phytoplankton, F (y, x) is called the invasion eigenvalue or invasion fitness. This measures the initial exponential growth rate of the mutant phytoplankton population in an environmental set dominated by the resident population [33]. If the mutants are rare and have a relatively small effect, then the cell size of the phytoplankton will evolve through both successive invasions and replacements along the direction predicted by the sign of the local selection gradient g(x), which can be obtained by taking the derivative of F (y, x) with respect to the mutant cell size y and evaluating it at the resident phytoplankton cell size x [33]. Then, The continuity of the selection gradient g(x) with respect to x and the small difference between cell sizes x and y in the early stage of mutation mean that the mutant's fitness can be linearly approximated by According to (3.4), a positive selection gradient g(x) allows the invasion of mutant phytoplankton with bigger cell sizes, whereas a negative g(x) is conducive to the successful invasion of mutant phytoplankton with smaller cell sizes.
When mutants with higher fitness values appear, evolution occurs. Following the results of [1,9], if the mutation processes are homogeneous and the mutations are rare and have a relatively small effect, it can be assumed that the rate of change in the mean cell size of the phytoplankton will be proportional to the fitness gradient where r is the birth probability of mutant phytoplankton, σ 2 is the variance of the mutation distribution, and P * (x) is the equilibrium population density of phytoplankton satisfying (2.3). The factor (1/2)rσ 2 is called the mutation rate. Together with the population density P * (x), this determines the rate of evolutionary change [8]. Equation (3.5) governs the variation of the cell size.

Evolutionary analysis
One of the central issues of evolutionary dynamics is to characterize the evolutionary process at community levels, such as whether it is a continuously stable strategy, an evolutionary repeller, or some other process. According to Definition 3.1, the evolutionarily singular strategy x * satisfies g(x * ) = 0. In addition to the invasion fitness (3.2) and fitness gradient (3.3), we now study the evolutionarily singular strategy in terms of the basic reproductive index defined in (2.2).
Proof. The evolutionarily singular strategy x * occurs when g(x * ) = 0. From (2.3), it follows that A direct calculation then leads to This completes the proof. Theorem 3.1 tells us that the evolutionarily singular strategy x * appears at the extreme point of the function R 0 (x).

Definition 3.2 ( [33]
). An evolutionarily singular strategy x * is said to be convergence stable if the repeated invasion of nearby mutant strategies into nearby resident strategies will lead to the convergence of resident strategies x * .
The condition for convergence stability is given by [10,12,13] where g(x) is the selection gradient defined by (3.3).

Definition 3.3 ( [33]
). An evolutionarily singular strategy x * is said to be evolutionary stable if x * is a maximum fitness point with respect to the mutant trait value y.
The evolutionarily singular strategy x * is stable against the invasion of neighboring strategies. We consider the second derivative of the invasion fitness with respect to the mutant trait value y, and evaluate it at the singular point x * . The condition for evolutionary stability is then given by [7,12,27]

Definition 3.4 ( [33]
). An evolutionarily singular strategy x * is said to be a continuously stable strategy if x * is both evolutionary stable and convergence stable.

Definition 3.5 ( [12]
). An evolutionarily singular strategy x * is said to be an evolutionary repeller if x * is neither convergence stable nor evolutionary stable.
Convergence stability and evolutionary stability identify the dynamical behavior of evolutionarily singular strategy x * . We next analyze the convergence stability and evolutionary stability of the evolutionarily singular strategy x * with the help of the ecological reproductive index. Theorem 3.2. Let x * be an evolutionarily singular strategy of (3.5).
Proof. For system (3.5), we have From (2.4) and (3.7), it follows that By (3.6), (3.9) and (3.10), if ∂ 2 R 0 /∂x 2 | x=x * < 0, then and so x * is a continuously stable strategy. If ∂ 2 R 0 /∂x 2 | x=x * > 0, then and x * is an evolutionary repeller. In general, for the given trade-off functions, on the curve defined by the ecological reproductive index with respect to the cell size, the evolutionarily singular points x * of system (3.5) at the local "peaks" of the curve denote the continuously stable strategies. The evolutionarily singular points x * of system (3.5) at the local "valleys" of the curve denote the evolutionary repellers. Theorems 3.1 and 3.2 reveal that the phytoplankton cell size will always evolve to maximize its fecundity for the current environment.
The evolutionary dynamics of system (3.5), characterized by Theorems 3.1 and 3.2, can be intuitively represented with the help of the basic reproductive index and pairwise invisibility plots. Figure 1 illustrates such a scenario, where the parameter values and trade-off functions are deliberately specified such that the curve defined by R 0 (x) admits both "peaks" and "valleys". In Figure 1-(a), Theorem 3.1 implies that the two solid points and one hollow point satisfy ∂R 0 /∂x| x=x * = 0, and thus denote evolutionarily singular strategies. By Theorem 3.2, the two solid points at the local "peaks" of the curve satisfying ∂ 2 R 0 /∂x 2 | x=x * < 0 are continuously stable strategies, whereas the hollow point at the local "valley" of the curve satisfying ∂ 2 R 0 /∂x 2 | x=x * > 0 is an evolutionary repeller. The pairwise invisibility plot ( Fig.  1-(b)) shows that, in the dark (+) region, the mutants have positive fitness values, which implies that the mutant population can invade; in the white (−) region, the mutant have negative fitness, which means the mutant population can not survive. In particular, in the dark (+) region on the left side of the vertical line through x * 2 , x * 1 is convergence stable and the local fitness gradient points toward the singular strategy x * 1 ; in the dark (+) region on the right side of the vertical line through x * 2 , x * 3 is convergence stable and the local fitness gradient points toward the singular strategy x * 3 .

Impact of temperature on evolution of cell size
In this section, we analytically investigate how the evolutionary model responds to variations in water temperature.
Then, the evolutionarily singular point x * of (3.5) is decreasing with respect to water temperature T .

Proof. From Theorem 3.1 and (3.3), it follows that
By (2.4), we have Then, which implies that the effect of temperature on the evolutionarily singular point x * is related to the temperature-dependent loss rate. A direct calculation produces Hence, (4.1) implies that k(T ) and H(x * ) have opposite monotonic trends. Equation (4.3) indicates that x * is decreasing with respect to T when (4.1) holds.
We next discuss the biological significance of (4.1). It is assumed that phytoplankton cells are spherical. Therefore, the loss rate of a cell with a ESD of x increases in proportion to the square of its diameter [16,17]. This fact implies that s (x * ) > 0 and s (x * ) > 0. For species with a relatively large cell size, as the cell size increases, the maximum specific growth rate tends to decline with a decelerating rate [16], presumably because it is more efficient for smaller cells to acquire resources [25]. An examination of picophytoplankton shows that the maximum specific growth rate tends to increase with size, mainly for the limited supply of cellular catalysts in extremely small cells [25]. These patterns require µ opt (x) to satisfy µ opt (x) > 0 and µ opt (x) < 0 for small x, and µ opt (x) < 0 and µ opt (x) > 0 for relatively large x. Equation (4.2) together with positive s (x * ) results in which implies that the evolutionarily singular point x * should fall in the range where µ opt (x * ) > 0 and µ opt (x * ) < 0. The monotonicity of s(x * ) and µ opt (x * ) leads to The loss rate is temperature-dependent and various candidates have been suggested to account for this dependence [17]. The temperature-dependent loss rate k(T ) is chosen to be a monotonic increasing function [17]. Then we have Equations (4.4) and (4.5), deduced by empirical patterns, ensure that (4.1) holds.
To clarify the effect of temperature on the evolution of phytoplankton cell size, (2.1) should be specified using some empirical models. Table 1 illustrates some typical examples of the functions in (2.1) satisfying (4.1) of Theorem 4.1. In Table  1, the effect of temperature on the growth of phytoplankton is described by the so-called cardinal temperature model with inflexion (CTMI) [26], which has been validated for a large range of phytoplankton species [6]. In addition, the nutrient quota of a cell is assumed to be directly proportional to its volume, i.e., the cube of its ESD [16]. The parameter values are listed in Table 2.
µopt(x) cell size dependent growth rate µopt(x) = x a 1 x 2 + a 2 x + a 3 [16] s(x) cell size dependent loss rate s(x) = αx 2 [16] h(T ) temperature dependent growth rate [3] k(T ) temperature dependent loss rate  Figure 2 shows the results of numerical simulations of the ecological reproductive index with respect to cell size along the gradient of water temperature T , from low temperature (e.g., 4 • C) to high temperature (e.g., 28 • C). From Theorems 3.1 and 3.2, we know that the red solid points at the peaks of the curves defined by R 0 (x) represent the continuously stable strategies x * at different temperatures.
LetR 0 be the basic reproductive value at the given temperature T and the corresponding continuously stable strategy x * , i.e.,R 0 = R 0 (T, x * ), which measures the evolutionarily stable fecundity of phytoplankton at different temperatures. When the water temperature is below the optimal level (e.g., 20 • C), the evolutionarily stable fecundity of phytoplanktonR 0 increases, but the cell size x * decreases with an increase in water temperature (Fig. 2-(a)). When the water temperature is above the optimal level, bothR 0 and x * decrease as the temperature increases ( Fig.  2-(b)). Moving the temperature away from the optimal level has a negative effect on the fecundity of phytoplankton. . Basic reproductive index with respect to cell size x with varying temperature. The solid point at the local "peak" of the curve represents the continuously stable strategy. (a) Water temperature below the optimal level. (b) Water temperature above the optimal level. The horizontal arrows represent the directions as x * evolves with increasing temperature. The vertical arrows represent the directions as the evolutionarily stable fecundity of phytoplanktonR0 varies with increasing temperature. Figure 3 gives a full picture of the effect of water temperature on the evolutionarily stable fecundity of phytoplanktonR 0 and evolutionarily stable phytoplankton cell size x * . Figure 3-(a) shows that the evolutionarily stable fecundity of phytoplankton first increases and then decreases as the water temperature rises. The evolutionarily stable fecundity of phytoplankton achieves a maximum value at some optimal temperature. Figure 3-(b) elaborates the relationship between cell size and water temperature. It shows that the phytoplankton cell size decreases as the temperature rises, and is always above a positive minimum level.

Conclusions
Adaptive dynamics are an effective approach for studying the evolutionary phenotypic changes in evolving populations when the fitness is density-or frequency-dependent [12]. This study has investigated the evolutionary dynamics of a single population of phytoplankton with varying water temperature.
A nutrient-phytoplankton model incorporating the effect of temperature and cell size was first formulated and then fully analyzed. Based on the globally asymptotically stable internal equilibrium, we proposed an evolutionary model of cell size and examined its evolutionary dynamics by investigating both the convergence stability and evolutionary stability. Equation (3.5) is much more general and flexible in terms of admitting different concrete functional forms of cell size and temperature. The criteria established to characterize the evolutionary dynamics of (3.5) (i.e.,Theorem 3.1 and 3.2) are independent of the functional forms of the cell size and temperature. It is the limiting factors, not their concrete functional forms, which assert a significant impact on the phytoplankton dynamics.
The ecological reproductive index R 0 measures the fecundity of phytoplankton and plays an important role in characterizing their evolutionary dynamics. Theoretical analysis has revealed that the evolutionarily singular strategy x * of system (3.5) occurs at the points where the derivative of R 0 is equal to zero. Moreover, the evolutionarily singular strategy x * of (3.5) at the local "peaks" of the curve defined by R 0 (x) is a continuously stable strategy, and the evolutionarily singular strategy at the local "valleys" is an evolutionary repeller. Numerical simulations have been conducted to confirm the theoretical analysis.
The important insights obtained from accurate computations indicate that the continuously stable strategy is affected by the temperature dependent loss rate rather than the growth rate. To adapt to increasing temperatures, the cell size of phytoplankton evolves to decrease and maximize its fecundity, but never falls below a minimum value. The continuous stable fecundity depicted by the ecological reproductive index first increases and then decreases as the water temperature rises.
In summary, our analytical and numerical analysis suggests the following conclusions: • In mathematical terms, unlike the traditional analysis of evolutionary dynamics using the invasion fitness and fitness gradient, this paper has proved that the ecological reproductive index is important in determining the type of singular points attained by the phytoplankton evolutionary model and characterizing the evolutionary dynamics.
• In ecological terms, it has been revealed that the evolution of the phytoplankton cell size is affected by the temperature-dependent loss rate. When the temperature increases within a reasonable range, the cell size of phytoplankton evolves to decrease. Thus, global warming will cause phytoplankton to evolve into picophytoplankton.
As phytoplankton is an assembly of organisms that arose from different evolutionary processes and present different characteristics (e.g., physiological, morphological), the extension of our evolutionary framework to include several other traits and abiotic factors should enable a more comprehensive understanding of the dynamics of aquatic ecosystems. A key challenge for future research is to investigate not only the impact of temperature on the evolution of primary producers (i.e., phytoplankton), but also the influence of comprehensive climate changes on the networks of evolutionary changes throughout the food web.