JSTOR

Primary Producer Dynamics

Our model assumes that the specific nutrient uptake rate of a primary producer, f_i[N, Q_i], is a function of the environmental nutrient concentration N and of its intracellular nutrient content Q_i (Morel 1987; Ducobu et al. 1998): where f_max, i is the maximum nutrient uptake rate of species i, K_N, i is its half‐saturation constant, Q_min, i is its minimum intracellular nutrient content, and Q_max, i is its maximum intracellular nutrient content. In this equation, the nutrient uptake rate of the primary producer increases with environmental nutrient availability according to Michaelis‐Menten kinetics, is highest when primary producers are starved (Q_i = Q_min, i), and reduces to 0 when primary producers are satiated with nutrients (Q_i = Q_max, i). Experimental support for these assumptions is provided in several studies of the nutrient uptake kinetics of phytoplankton species (e.g., Morel 1987; Ducobu et al. 1998; Passarge et al. 2006).

The intracellular nutrient content (also known as nutrient quota) of primary producers increases due to nutrient uptake and declines due to dilution by growth (Droop 1973): where μ_i[Q_i] is the specific growth rate of primary producer species i. The specific growth rate, in turn, is an increasing function of the intracellular nutrient content according to the Droop equation (Droop 1973; Grover 1991): where μ_max, i is the maximum specific growth rate of species i. This equation states that the specific growth rate of the primary producer is positive if its intracellular nutrient content exceeds the minimum value Q_min, i. That is, Q_min, i corresponds to the intracellular nutrient content of the primary producer that just secures its basal metabolism and survival, while higher nutrient contents enable population growth (Sterner and Elser 2002).

The population dynamics of the primary producers are driven by their growth rates, mortality rates, and the grazing rate by herbivores: where ; P_i is the population abundance of primary producer species i; d_i is its mortality rate; Z is the herbivore abundance; and n is the number of phytoplankton species. The grazing rate on species i, g_i[P, Q], is a function of the population abundances and nutritional quality of all primary producers in the system, where the bold P and Q represent vector notation.

Herbivore Dynamics

Grazing by herbivores depends on the nutritional quality of the primary producers. We consider an evolutionary scenario in which all primary producer species are identical except for their nutrient uptake rate and, hence, their intracellular nutrient content. We therefore assume that the herbivore has the same search rate, a, and handling time per prey item, h, for all primary producer species. However, the herbivore has a preference for the most nutritious species. The grazing rate of the herbivore on primary producer species i can then be modeled as a multispecies type II functional response (Holling 1959; Křivan 1996): where α_i[P, Q] is the probability that the herbivore will consume encountered individuals of primary producer species i, given the population abundances and nutritional quality of all primary producers in the ecosystem. The underlying assumption that herbivores can detect the nutritional value of their food is supported by several experimental studies (Cowles et al. 1988; DeMott 1989; Buskey 1997; John and Davidson 2001; Schatz and McCauley 2007; Martel 2009). For instance, various zooplankton species can sense the nutrient content of their prey through chemical cues such as the chemical composition of their outer cell surface (Martel 2009) and leakage of nutrient‐rich molecules (Buskey 1997; Kiørboe 2008).

Optimal foraging theory predicts that herbivores will display a stepwise switch to less nutritious prey when the more nutritious prey has been reduced below a critical threshold abundance (e.g., Charnov 1976). However, this stepwise switch emerges from optimal foraging theory only if the forager is a perfect predator with complete and immediate knowledge of all prey types (Stephens and Krebs 1986; McNamara and Houston 1987). Deviations from these idealized assumptions lead to gradual switches in diet, for instance, if predators need some recognition time to identify their prey (McNamara and Houston 1987) or detect their prey only within a limited range (Berec and Křivan 2000). Indeed, in reality, predators often show gradual rather than stepwise changes in diet composition (Davies 1977; Mittelbach 1981; Stephens and Krebs 1986). Following previous studies, we therefore model gradual switching behavior of the predator by a sigmoid function (Egas et al. 2004; Rueffler et al. 2007). We assume that the probability that the herbivore will consume primary producer species i, α_i[P, Q], depends on the nutritional quality of focal species i in comparison to the average nutritional quality of all primary producers in the community: where S measures the selectivity of the herbivore for more nutritious primary producers and Q_mean[P, Q] is the average nutritional quality of the primary producers in the community: The shape of equation (6) is illustrated in figure 1. We note that if the herbivore is nonselective ( ), then each primary producer species is consumed with the same probability irrespective of its nutritional quality. Conversely, if the herbivore is highly selective ( ), then it switches abruptly to prey species exceeding average nutritional quality.

(29 KB)

Figure 1: Probability that the herbivore consumes primary producer species 1 (α₁) as a function of the nutritional quality of primary producer species 1. The curves are derived from equation (6), and each curve corresponds to a different selectivity of the herbivore ( , 0.1, 0.5, 1, and 5 cells fmol⁻¹). For the purpose of illustration, we assume that the average nutrient content of all primary producers in the community is fmol cell⁻¹.

Open New Window

We assume that the herbivore assimilates the ingested primary producers with efficiency e_i (Loladze et al. 2000): which is defined as the ratio of the nutrient content of primary producer species i to the nutrient content of the herbivore, q. Accordingly, high nutritional quality of primary producers or low nutritional demands of the herbivore will result in high assimilation efficiencies. Herbivores usually have a higher nutrient content per unit carbon than primary producers (Elser et al. 2000). When expressed on a per‐unit carbon basis, this precludes assimilation efficiencies greater than 1.

The population dynamics of the herbivore may then be written as where m is the specific mortality rate of the herbivore.

Nutrient Dynamics

We assume that the nutrient available in the environment is consumed by primary producers and recycled by both primary producers and herbivores (Grover 1997): where the first term on the right‐hand side describes nutrient uptake by the primary producers and the other two terms describe nutrient recycling due to the mortality of primary producers and herbivores, respectively.

The total amount of nutrient in the ecosystem, T, includes the freely available nutrient in the environment as well as the nutrient contained in the primary producers and herbivore: Evaluation of the time derivative of this equation shows that the total amount of nutrient remains constant (i.e., ). In other words, our model ecosystem is a closed system with respect to nutrients.

Adaptive Dynamics of Nutrient Uptake

To model evolutionary changes in nutrient uptake rate, we consider a resident community consisting of one limiting nutrient, one primary producer species, and one herbivore species. In the resident population of the primary producer, a novel mutant phenotype may appear. The mutant resembles the resident primary producer in every respect except its maximum nutrient uptake rate, f_max. This evolving trait indicates the active uptake rate of nutrients across cell membranes and generally correlates with the number of transport proteins assembled for nutrient uptake (Cornish‐Bowden 1995).

We will indicate the mutant and resident by the subscripts m and r, respectively. Initially, the mutant phenotype is rare relative to the resident phenotype. The success of a mutant will therefore depend on its invasion fitness, w, defined as the net specific growth rate of the mutant in the resident population: If the invasion fitness is negative ( ), the mutant will become extinct. Conversely, if the invasion fitness is positive ( ), the mutant will thrive and establish a new population. This new population can then again be invaded by another mutant phenotype, and so on. Following the rationale of adaptive dynamics (Metz et al. 1992; Geritz et al. 1998), we assume that mutations are rare events and that the phenotype of the mutant differs only slightly from that of the resident phenotype. Thus, trait evolution proceeds gradually, as a series of successful invasions by mutant phenotypes.

Two standard assumptions are made to derive invasion fitness. First, we assume that the population dynamics of the mutant phenotype obeys the same rules as the resident phenotype. Substituting equation (4) into equation (12), invasion fitness can be written as Second, we assume that the mutant phenotype is sufficiently rare not to affect the population dynamics of the resident community (i.e., ). Hence, the average nutritional quality of the community is determined by the resident (i.e., , which implies ). Substituting equations (3) and (5)–(7) into equation (13), the invasion fitness then reads This equation shows that an increase of the intracellular nutrient content of the mutant will increase its specific growth rate (the first term on the right‐hand side) but will also increase the mutant’s susceptibility to grazing compared to the resident (the last term on the right‐hand side). In addition, the last term shows that invasion fitness depends on herbivore selectivity. If the mutant has a higher intracellular nutrient content than the resident, an increase of herbivore selectivity will increase the mutant’s susceptibility to grazing. The intracellular nutrient contents of the mutant and resident, in turn, depend on their maximum nutrient uptake rate (via eqq. [1], [2]), which is the evolving trait in our analysis.

Trait evolution resulting from successful invasions can be analyzed graphically using pairwise invasibility plots and mutual invasibility plots (Geritz et al. 1998). For this purpose, we calculated invasion fitness numerically for various combinations of mutant and resident phenotypes. The resident community may be at equilibrium or may display nonequilibrium dynamics. If the resident community reached equilibrium, then we calculated the intracellular nutrient content of the mutant in equilibrium with the nutrient availability set by the resident community. According to equation (2), this yields We note that . That is, an increase in the maximum nutrient uptake rate of the mutant will result in a higher intracellular nutrient content and, therefore, a higher susceptibility to grazing. Subsequently, invasion fitness was calculated by inserting the equilibrium values of the resident community and the above intracellular nutrient content of the mutant into equation (14). If the resident community shows persistent nonequilibrium dynamics, then the net specific growth rate of a mutant may vary between positive values when nutrients are plentiful and herbivores are scarce and negative values when the environmental conditions worsen. In this case, we ran a full simulation of the resident community over a long time span, from to days. Subsequently, we calculated the long‐term average of the net specific growth rate of a mutant invading the fluctuating resident community to assess its invasion fitness in a variable environment (Metz et al. 1992). Numerical simulations were based on the NDSolve routine in Mathematica 6.0 (Wolfram Research, Champaign, IL).

Parameterization of the Model

Our model is parameterized for plankton communities. This allows realistic choices of parameter values, because many of the process descriptions in the model have been verified and measured for phytoplankton and zooplankton species (table 1). However, many of our results may probably be generalized to other plant‐herbivore communities as well.

Table 1: Model variables and parameter values

Open New Window

Phosphorus is an essential element in phytoplankton and zooplankton metabolism and an important limiting nutrient in many natural waters (Schindler 1977; Wu et al. 2000; Elser et al. 2007). Parameter values were therefore obtained from studies on phosphorus‐limited phytoplankton (Chlorella vulgaris; Passarge et al. 2006). We choose a copepod species (Acartia tonsa; Kiørboe et al. 1985; Walve and Larsson 1999; Hirst and Kiørboe 2002) as the herbivore, since copepods tend to be more selective grazers than generalist filter feeders such as Daphnia. Model variables and parameters are listed in table 1.

Ecological Dynamics

To set the stage, we first consider a community consisting of only one primary producer species and one herbivore and investigate the ecological dynamics of this simple community along a productivity gradient. In line with expectation, the model predicts low productivity of the primary producer at low nutrient levels, supporting an equilibrium community of plants and herbivores. Nutrient enrichment results in higher growth rates and higher nutritional quality of the primary producer, which can drive the community from a stable equilibrium toward plant‐herbivore oscillations. This pattern matches Rosenzweig’s (1971) classic paradox of enrichment and has also been found in other stoichiometrically explicit plant‐herbivore models (e.g., Diehl 2007; Jäger et al. 2008). A similar pattern can be observed for variation in maximum nutrient uptake rate of the primary producer. When maximum nutrient uptake rate is low, the primary producer grows poorly and supports low herbivore densities (fig. 2A). An increase in maximum nutrient uptake rate enhances the growth rate and nutrient content of the primary producer, which may shift the community dynamics from a stable equilibrium to plant‐herbivore oscillations (fig. 2B). These ecological dynamics are summarized in figure 3A for different combinations of the total nutrient concentration, T, and maximum nutrient uptake rate of the primary producers, f_max.

(31 KB)

Figure 2: Community dynamics resulting from different maximum nutrient uptake rates. A, A low maximum nutrient uptake rate of the primary producer ( pmol cell⁻¹ day⁻¹) leads to an equilibrium with low‐nutritional plants and low herbivore abundance. Note that the nutrient line is partly masked by the primary producer. B, In contrast, a high maximum nutrient uptake rate ( pmol cell⁻¹ day⁻¹) leads to high nutritional quality of the plants and plant‐herbivore oscillations. The total nutrient concentration was μM. Other parameters as in table 1.

Open New Window

(34 KB)

Figure 3: Ecological and evolutionary dynamics along a productivity gradient. A, Nutrient enrichment and high maximum nutrient uptake rates of the primary producer shift the ecological dynamics from a stable equilibrium of plants and herbivores to plant‐herbivore oscillations (limit cycles). Solid lines indicate the boundaries between regions with different ecological dynamics, as determined by numerical bifurcation analysis. Between the empty system (ES) and stable plant‐herbivore equilibrium, a very narrow region (too narrow to be visible in the graph) consists of primary producers only. The dashed line represents the maximum nutrient uptake rate of the CSS phenotype, . This shows an evolutionary reduction of the maximum nutrient uptake rate with an increasing total amount of nutrient in the ecosystem. The CSS phenotype produces stable plant‐herbivore communities (solid symbols) in nutrient‐poor ecosystems but plant‐herbivore oscillations (open symbols) in nutrient‐rich ecosystems. B–D, The amplitude of the plant‐herbivore oscillations generated by the continuously stable strategy phenotype increases with the total nutrient concentration. Open symbols in B–D represent the minima and maxima of the plant‐herbivore oscillations. The simulations assume mildly selective grazing ( cells fmol⁻¹). Other parameters as in table 1.

Open New Window

Evolutionary Dynamics

As a next step, we consider invasion of novel mutant phenotypes of the primary producer into the resident plant‐herbivore community. Pairwise invasibility plots (PIPs) show which mutant phenotypes can invade which resident phenotypes (fig. 4). If invasion fitness is positive for mutants with a higher maximum nutrient uptake rate than the resident phenotype, then the maximum nutrient uptake rate will increase during evolution. Conversely, if invasion fitness is positive for mutants with a lower maximum nutrient uptake rate than the resident phenotype, then the maximum nutrient uptake rate will decrease. This results in trait evolution, which is conceived as a gradual process that occurs in small steps. Evolution therefore proceeds along the diagonal of the PIP.

(31 KB)

Figure 4: Selective grazing and the evolution of nutrient acquisition. A–C, Pairwise invasibility plots. Black and gray shading indicate regions with positive and negative invasion fitness, respectively; white dots represent continuously stable strategy (CSS) phenotypes. A, Under nonselective grazing ( cells fmol⁻¹), evolution leads to runaway selection toward high maximum nutrient uptake rates. B, Under mildly selective grazing ( cells fmol⁻¹), evolution converges to a CSS phenotype with an intermediate maximum nutrient uptake rate. C, Under strongly selective grazing ( cells fmol⁻¹), evolution converges to a CSS phenotype with a low maximum nutrient uptake rate. D, The maximum nutrient uptake rate of the CSS phenotype as function of herbivore selectivity. The gray horizontal band indicates the range of f_max values that yield a stable equilibrium of plants and herbivores. Values of f_max above this gray band lead to plant‐herbivore oscillations, while the herbivore cannot persist if f_max falls below this gray band. Mildly selective herbivores favor CSS phenotypes that produce plant‐herbivore oscillations (open symbols ), whereas highly selective herbivores stabilize the community dynamics (solid symbols). The total nutrient concentration was μM. Other parameters as in table 1.

Open New Window

Our results show that evolution of the maximum nutrient uptake rate of the primary producer depends on the selectivity of the herbivore (fig. 4). If the herbivore is nonselective, mutants have positive invasion fitness whenever they have higher nutrient uptake rates than residents (fig. 4A). This favors runaway evolution toward high maximum nutrient uptake rates. However, if the herbivore is selective, its preference for more nutritious primary producers favors evolution toward intermediate or low maximum nutrient uptake rates (fig. 4B, 4C). Thus, selection for growth at low maximum nutrient uptake rates is countered by selection against grazing at high maximum nutrient uptake rates, and the evolutionary dynamics converge to an intermediate phenotype. In PIPs, this phenotype is found at the intersection of the two zero isoclines separating regions of positive and negative invasion fitness. In our model, this singular phenotype, hereafter written as , cannot be invaded by any other phenotype. Hence, is an evolutionarily stable strategy (ESS; sensu Maynard Smith and Price 1973). We note that is not only an ESS but also a continuously stable strategy (CSS; sensu Eshel and Motro 1981), because it constitutes an evolutionary attractor. The exact value of will depend on the traits of the primary producers and herbivores as well as the environmental conditions. We did not find examples of evolutionary branching in our simulations. This implies that once evolution has converged to the CSS phenotype, the ecological dynamics is captured by a simple community consisting of a single primary producer phenotype and a herbivore population.

We investigated the ecological dynamics of primary producers and herbivores and the evolutionary changes in maximum nutrient uptake rate of primary producers along a productivity gradient, assuming a mildly selective herbivore (fig. 3). On ecological timescales, nutrient enrichment would shift the population dynamics from stable equilibria to plant‐herbivore oscillations and would improve the nutritional quality of primary producers. Evolution counters this effect by selection for a lower maximum nutrient uptake rate in ecosystems with a higher total amount of nutrient (fig. 3A). Nevertheless, despite this reduction in nutrient uptake rate, nutrient enrichment still yields sufficiently high nutritional quality and growth rates of primary producers to drive the plant‐herbivore community toward oscillations (fig. 3). Apparently, if herbivores are only mildly selective, evolutionary reduction of the maximum nutrient uptake rate is not strong enough to prevent destabilization of the ecological dynamics through nutrient enrichment.

However, this result depends on herbivore selectivity. As the selectivity of the herbivore increases, the fitness cost of a high nutrient content of the primary producer also increases, and is therefore reduced to lower values (fig. 4). Hence, highly selective herbivores suppress the maximum nutrient uptake rate of the CSS phenotype more strongly. Yet, selective herbivores cannot suppress the maximum nutrient uptake rate beyond a lower bound. The gray horizontal band in figure 4D indicates the narrow range of values that yield stable equilibria (corresponding to the narrow range of stable equilibria at μM in fig. 3A). If would be suppressed below this gray band, the herbivore population would become extinct (fig. 3A) and evolution would increase the maximum nutrient uptake rate of primary producers once the herbivores are gone. In other words, highly selective herbivores drive the CSS phenotype into this narrow gray band and thereby stabilize the community dynamics in nutrient‐rich environments (fig. 4D).

Coexistence or Competitive Exclusion?

The results above assumed that trait evolution is a gradual process taking small steps only. However, one can relax this premise to consider invasion of novel phenotypes radically different from the resident phenotype. This could offer opportunities for coexistence that cannot be reached by gradual evolution (Egas et al. 2004). As a rule of thumb, two phenotypes can coexist if they are each able to invade an established resident population of the other phenotype, a condition termed mutual invasibility (Geritz et al. 1998). Opportunities for coexistence can be deduced from mutual invasibility plots (MIPs; fig. 5). A MIP is obtained by superposition of two PIPs, one indicating whether phenotype 1 can invade a resident population of phenotype 2 and the other indicating whether phenotype 2 can invade a resident population of phenotype 1 (Geritz et al. 1998).

(35 KB)

Figure 5: Pairwise invasibility and mutual invasibility. A, D, Invasion of phenotype 2 (P₂) into the resident population of phenotype 1 (P₁). B, E, Invasion of phenotype 1 into the resident population of phenotype 2. Black and gray shading in these pairwise invasibility plots indicate regions with positive and negative invasion fitness, respectively; white dots represent the continuously stable strategy phenotype. C, F, The resulting mutual invasibility plots (MIPs). A–C, If grazing is nonselective ( cells fmol⁻¹), each phenotype can be invaded by phenotypes with a higher maximum nutrient uptake rate. The MIP in C shows that this leads to competitive exclusion, where the phenotype with the highest maximum nutrient uptake rate always wins. D–F, If grazing is strongly selective ( cells fmol⁻¹), phenotypes with high maximum nutrient uptake rates can be invaded by phenotypes with lower maximum nutrient uptake rates. The MIP in F shows that in a large part of the trait space phenotypes with lower maximum nutrient uptake rates win against phenotypes with higher maximum nutrient uptake rates. In addition, the MIP in F shows a narrow region in trait space, indicated in black, where the two phenotypes coexist. The total nutrient concentration was μM. Other parameters as in table 1.

Open New Window

In our model, selective grazing is decisive in determining the outcome of mutual invasibility. If the herbivore is nonselective, the phenotype with the highest maximum nutrient uptake rate always wins (fig. 5A–5C). However, if the herbivore is selective, phenotypes with high maximum nutrient uptake rates can be displaced by phenotypes with lower maximum nutrient uptake rates, which are eaten less due to their lower nutritional quality (fig. 5D–5F). In case of selective grazing, a narrow region in the MIP displays species coexistence, where both phenotypes can invade each other (fig. 5F). Species coexistence is readily confirmed with a food‐web model comprising two primary producers (fig. 6). One primary producer has a high maximum nutrient uptake rate and will be called the nutrient exploiter. The other primary producer has a low maximum nutrient uptake rate and will be called the grazing avoider. We note that the community dynamics produces plant‐herbivore oscillations. Thus, in this example, the food‐web model does not predict stable coexistence but nonequilibrium coexistence of the primary producers (fig. 6).

(31 KB)

Figure 6: Nonequilibrium coexistence of the grazing avoider and nutrient exploiter. The grazing avoider and nutrient exploiter are identical except for their maximum nutrient uptake rate (grazing avoider: pmol cell⁻¹ day⁻¹; nutrient exploiter: pmol cell⁻¹ day⁻¹). The simulations assume mildly selective grazing ( cells fmol⁻¹). The total nutrient concentration was μM. Other parameters as in table 1.

Open New Window

A Stoichiometric Trade‐Off

In our model, enhanced nutrient uptake rates improve the competitive ability of primary producers at the expense of an increased susceptibility to grazing (see also de Mazancourt et al. 2001; Loeuille and Loreau 2004), so that phenotypes with an intermediate nutrient uptake rate are selected. We did not specify an explicit trade‐off between two or more species traits. Instead, the trade‐off between competitive ability and grazing susceptibility resulted from a change in a single trait, the maximum nutrient uptake rate. Trade‐offs driven by changes in a single trait are described as “emergent trade‐offs” in the epidemiological literature (André et al. 2003; van Ballegooijen and Boerlijst 2004; Alizon and van Baalen 2005). For instance, the replication rate of a parasite is a single trait that favors both its transmission and its virulence. This has contrasting effects on parasite fitness, since a fast replication rate ensures high parasite transmission but reduces host survival (Alizon and van Baalen 2005). Consequently, a trade‐off emerges, and intermediate replication rates of the parasite are selected. In this respect, our model bears a remarkable resemblance to host‐parasite models and, to our knowledge, is among the first to address an emergent trade‐off in the ecological stoichiometry of primary producers (see also Menge and Weitz 2009).

Interestingly, recent work provides experimental support for the stoichiometric trade‐off predicted by our study. Passarge et al. (2006) investigated the competitive ability for phosphorus and light of five freshwater phytoplankton species and showed that the strongest competitors had the highest intracellular phosphorus contents. Although the grazing susceptibility of their species is not known, these results suggest that strong competitors would constitute more nutritious food for zooplankton and are likely to suffer more from grazing than weaker competitors. Sunda and Hardison (2010) reported considerable variation in the ammonium uptake rates of marine phytoplankton species. Four of their algal isolates are known to be poorly grazed by zooplankton, and only these four isolates had unusually low ammonium uptake rates and associated growth rates for their size. These studies present two intriguing examples of a trade‐off between competitive ability and grazing susceptibility, driven by interspecific differences in nutrient acquisition.

Hence, low nutritional quality of primary producers can be an advantageous defense against selective herbivores (Moran and Hamilton 1980). However, primary producers may also adopt alternative defense mechanisms to cope with selective grazing, such as the production of secondary compounds (Stamp 2003; Van de Waal et al. 2009), development of spines and thorns (Milewski et al. 1991; Tollrian and Harvell 1999), and reinforcement of cell walls (Van Donk et al. 1997; Hamm et al. 2003). These defenses may enable primary producers to sustain high nutritional quality without suffering intense grazing by selective herbivores. Primary producers that invest in alternative defense mechanisms may thus overcome a stoichiometric trade‐off between competitive ability and grazing susceptibility. The interplay between the stoichiometric trade‐off described here and investments in alternative defense mechanisms offers an interesting avenue for further research.

Evolutionary Convergence to Nonequilibrium Dynamics

Our results show that evolution of the nutrient uptake rate allows sustained oscillations of the limiting nutrient, primary producers, and herbivores. This result adds to a long‐lasting debate on the evolution of the stability of ecological communities. Theory predicts oscillations in predator‐prey communities (Lotka 1925; Volterra 1926; Rosenzweig 1971), and predator‐prey oscillations have been reported in many experimental studies (Gause 1934; Fussmann et al. 2000; Yoshida et al. 2003; Benincà et al. 2009). Most predator‐prey models, however, also predict stable equilibria for at least part of their parameter space. Hence, the question is whether evolution will drive predator‐prey interactions toward ecological stability or nonequilibrium dynamics.

Several theoretical studies have suggested a general tendency of natural selection to stabilize population dynamics, with nonequilibrium dynamics evolving only under stringent trade‐off postulates (Doebeli and Koella 1995; Ebenman et al. 1996; Zeineddine and Jansen 2005; but see Ferrière and Gatto 1993). However, other theoretical studies have shown that evolution may destabilize predator‐prey interactions (Abrams and Matsuda 1997; Abrams 2000; Jones and Ellner 2007). Indeed, long‐term experimental studies with ample time for evolutionary changes displayed sustained predator‐prey fluctuations, indicating that evolution did not prevent nonequilibrium dynamics (Yoshida et al. 2003; Benincà et al. 2008).

We systematically explored our model and found that nutrient enrichment triggers evolutionary convergence toward lower maximum nutrient uptake rates (fig. 3A). For mildly selective herbivores, nutrient enrichment overrides the evolutionary reduction of the maximum nutrient uptake rate, which has a destabilizing effect on plant‐herbivore interactions. Accordingly, nutrient enrichment leads to plant‐herbivore oscillations despite evolutionary reduction of the maximum nutrient uptake rate (fig. 3A). This illustrates that Rosenzweig’s (1971) classic paradox of enrichment is rather robust to evolutionary changes in the nutrient uptake rates of primary producers, at least for mildly selective herbivores. Our model predictions therefore lend theoretical support to the persistence of nonequilibrium dynamics over evolutionary timescales.

However, highly selective grazing by herbivores can stabilize plant‐herbivore interactions (fig. 4D). This occurs because a higher selectivity by herbivores favors a stronger evolutionary reduction of the nutrient uptake rate of primary producers. A strong reduction in maximum nutrient uptake suppresses the nutritional quality and growth rate of primary producers, which both have a stabilizing effect on plant‐herbivore interactions. This result echoes many earlier findings that selective predation stabilizes predator‐prey interactions (e.g., Murdoch and Oaten 1975; Fryxell and Lundberg 1994). However, the interesting point in our work is that evolution plays a key role. The parameter range that allows ecological stability at high nutrient levels is very narrow (fig. 3A). Yet, if the herbivore is sufficiently selective, evolution drives the nutritional quality of plants to low values and parks the plant‐herbivore interactions in this narrow lane of ecological stability (fig. 4D).

Coexistence of Nutrient Exploiters and Grazing Avoiders

The stoichiometric trade‐off between competitive ability and grazing susceptibility enabled coexistence of nutrient exploiters and grazing avoiders (figs. 5F, 6). This result is in line with several plant‐herbivore models predicting coexistence of primary producers through a competition‐predation trade‐off (Armstrong 1979; Holt et al. 1994; Leibold 1996; Křivan 2003; Yoshida et al. 2003). These models usually assume trade‐offs between two or more traits. For instance, Yoshida et al. (2003) considered a trade‐off between the half‐saturation constant and food value of primary producers, such that strong competitors had a higher food value. Indeed, their model predicted coexistence of primary producers with different food values and competitive abilities. In contrast, our model considers the evolution of only a single trait (i.e., the maximum nutrient uptake rate) while all other species traits are held constant. Apparently, the trade‐off between competitive ability and grazing susceptibility driven by this single trait creates sufficient niche differentiation for species coexistence.

However, our model predicts that gradual evolution will converge to a single primary producer with an optimal nutrient uptake rate that cannot be invaded by any competitor. Furthermore, we did not observe evolutionary branching, a process that describes the divergence of phenotypes and could have permitted coexistence on evolutionary timescales (Geritz et al. 1998). These results are in line with those of Shoresh et al. (2008), who investigated the evolutionary dynamics of a standard resource competition model and found far fewer coexisting species than would be expected from the number of ecological niches. Thus, in the long run, evolution may restrict the coexistence of primary producers.

Prospects for Further Investigation

Our model is evidently a simplification of natural communities, which host numerous species competing for multiple resources and engaged in many trophic interactions. For instance, our model assumes that primary producers and herbivores are exclusively nutrient limited. In reality, nutrient enrichment may alleviate primary producers from nutrient limitation, shifting their growth toward light‐limited conditions (Huisman and Weissing 1995; Passarge et al. 2006). Likewise, herbivores may become limited by the energy content rather than the nutrient content of their food when feeding on nutritious plants (Sterner and Elser 2002; Hall 2004; Diehl 2007). It is not immediately obvious how an environmental gradient from nutrient to energy limitation will affect the evolution of nutrient acquisition, and it would be interesting to study these aspects in further detail.

Another simplification is that we studied variation in a single trait only. In reality, phytoplankton species face trade‐offs among multiple traits. This affects opportunities for species coexistence but may also have implications for the evolutionary dynamics. For instance, Litchman et al. (2007) showed that the maximum nutrient uptake rate, the half‐saturation constant of nutrient‐limited growth, and the minimum nutrient content of phytoplankton species are all positively correlated with each other. The evolutionary implications of such correlated traits are unclear and certainly merit further study.

Finally, we note that herbivores may adapt to evolutionary changes in the nutritional quality of primary producers. Hence, further study of the coevolution of their ecological stoichiometry could be worth pursuing (Abrams 2000; Loeuille and Loreau 2004).

Conclusions

Our model analysis shows that high nutrient uptake rates of primary producers favor their competitive ability but also enhance their food quality for herbivores. In the presence of mildly selective herbivores, theory predicts that these two selection pressures converge to intermediate nutrient uptake rates that are evolutionarily stable but induce plant‐herbivore oscillations at ecological timescales. Highly selective herbivores favor a stronger evolutionary reduction in the nutrient uptake rate of primary producers, which may stabilize the plant‐herbivore interactions. The model is firmly based on a mechanistic understanding of resource competition and the ecological stoichiometry of plant‐herbivore interactions, which may facilitate experimental tests of the model predictions. Empirical investigation of the trade‐off described in this article is likely to provide further insights into the evolution of the ecological stoichiometry of primary producers.