The coccolithophore E. huxleyi is widely distributed over the global ocean (e.g., Hagino et al., 2011), viable over a temperature range from 4 to 28°C (e.g., Watabe and Wilbur, 1966; Fielding, 2013), with maximum growth rates in the temperature range 18–25°C (Watabe and Wilbur, 1966; Zhang et al., 2014). In general, growth rates increase with temperature, with clear differences between Arctic and Atlantic strains (Daniels et al., 2014; Zhang et al., 2014). Zhang et al. developed thermal reaction norms (TRNs) for six E. huxleyi isolates originating from the central Atlantic near the Azores, Portugal (38°34′N; monthly SST range 16–22°C), and five isolates originating from coastal waters near Bergen, Norway (60°18'N; monthly SST range 6–16°C), all kept at 15°C in culture. The fitted growth rates for the Bergen isolates were higher in the range 7–22°C; they were higher for the Azores isolates in the range 26–28°C, with a crossover point near 24°C.

While isolates of E. huxleyi from various regions around the globe have a core set of common genes, there is considerable genetic variability across its global distribution (Hagino et al., 2011; Read et al., 2013). According to Read et al., “Genome variability within this species complex seems to underpin its capacity both to thrive in habitats ranging from the equator to the subarctic and to form large-scale episodic blooms under a wide variety of environmental conditions.” Thus, the cultures in Schlüter et al. (2014), originating from a single cell take from waters (~10°C) near Bergen, Norway, would not necessarily be expected to grow at the maximum rate observed for the species at 15°C, even though they were apparently growing in an exponential manner.

The model is formulated to simulate, as closely as possible, the experimental protocol followed during the laboratory studies (Schlüter et al., 2014). The main trait is the growth rate (d⁻¹), which is a function of a single environmental variable/stressor—temperature. The cultures were grown at three different pCO₂ levels: 400, 1100, and 2200 μatm, but in this model only the experiments at the “ambient” level, 400 μatm, are simulated. The maximum growth rate as a function of temperature has long been considered to be an important predictor of the rate of primary production, along with incoming irradiance and nutrient availability (Eppley, 1972; Bissinger et al., 2008). A recent analysis of observations of growth rate as a function of temperature specifically for E. huxleyi has been published (Fielding, 2013): non-zero growth rates have been observed over the range of temperatures from 2 to 27°C, with a very sharp decline at ~27°C as also observed by Schlüter et al. (2014). The most dense range of observations are for standard temperatures 10 and 15°C, where growth rates from near zero to the maximum at that temperature have been observed (Figure 2A from Fielding, 2013). Several fits to the 99th quantile of the data (i.e., 1% of the data points exceeded the fitted function) were carried out. We adopted the power law fit, which is the best fit according to the criteria used by Fielding. Figure 1 shows the dependence of growth rate on temperature for the power law fit (Fielding, 2013). The red line shows the power law fit that passes through the growth rate at 15°C observed by Schlüter et al. (2014) (1.15 d⁻¹, black diamond), and the vertical dashed red line shows its rapid drop off near 27°C. The dashed black line shows the dependence on temperature of the maximum growth rate observed for E. huxleyi, according to Fielding (2013).

In the model, genotypes are formulated in equal intervals along the trait axis, the growth rate μ(T) (d⁻¹), where T is the temperature. Thus, there are potential genotypes along the trait axis from μ_max = 0 to μ_max ≈ 2 d⁻¹ [corresponding to a temperature ~27°C, where the growth rate drops precipitously to zero (Fielding, 2013; Schlüter et al., 2014, Supplementary Information)].

The model is a simple exponential growth equation for each genotype i: (1)dNidt=μi(T)Ni where: N_i is the number of cells in genotype i, and μⁱ(T) is the growth rate of genotype i. At 15°C the maximum possible growth rate is μ_max = 1.29 d⁻¹ (Figure 1, Fielding, 2013) and, as shown by the solid diamond in Figure 1, the realized growth rate was 1.15 d⁻¹ (Schlüter et al., 2014).

During the laboratory experiments, 10⁵ cells were transferred every 5 days from the existing batch cultures into fresh culture medium to initiate the next batch culture. In the model, the total number of cells across all genotypes N_T (= ∑_iN_i) is “normalized” to 10⁵ every timestep (0.2 d) with the same normalization factor applied to each genotype. In the standard simulations, random mutation was allowed once each day: at 15°C the growth rate was 1.15 d⁻¹, so that a mutation occurred slightly more often than 1 per generation. Each mutation produced a single cell (in 10⁵ cells). Mutations were allowed at genotypes with growth rates less than and equal to 1.15 d⁻¹. Those mutant genotypes with growth rates less than 1.15 d⁻¹ grew less quickly than the original genotype—so were not “fixed.” Because of the normalization each timestep, those mutants consisted of less than 1 cell, so were set to zero (i.e., not fixed) when they dropped to 0.1 cells. When multiple non-zero genotypes exist, the mean realized growth rate across all non-zero genotypes is given as μ_mean = (∑_iN_iμ_i)/N_T.

According to Huertas et al. (2011), “Experimental measures of mutation rates in phytoplankton range from 10⁻⁵ to 10⁻⁷ mutations per cell per generation.” So the rate of one mutation per 1.15 generations in a culture of 10⁵ cells is at the high end of the published rates. We carried out sensitivity studies with mutations occurring both more or less frequently than 1 per day: generally the simulations at lower rates resulted in a slowing down of the increase in biomass of favorable mutations but without any material change in the final results.

The magnitude of each mutation (distance of mutant genotype along the growth rate trait axis from its parent genotype) was determined from a random number generator. Two cases were simulated. First, with a “flat” pdf where the new genotype was equally probable anywhere from the lowest to highest allowed growth rate for that temperature, and, second with a Gaussian normal pdf, where the width of the distribution across genotypes was changed for different simulations. The original genotype and genotypes resulting from mutations were tracked separately. For a flat pdf, only mutations from the original genotype were allowed, since the origin was not important. For the Gaussian pdf, mutations were allowed from the original genotype or from the mutant genotypes (with a probability proportional to their relative biomasses). In both cases, only one mutation was allowed per day; sensitivity simulations with more or fewer mutations per day did not materially affect the results.

To illustrate the description above without explanation until the Results section, Figure 2A shows an initial genotype at the start of a simulation at 15°C, and Figure 2B shows a histogram of the distribution (on a logarithmic axis) of total cell biomass in each genotype—after 3 years of simulation at 15°C. The blue bars represent the relative biomass in all non-zero genotypes; “extinct” or not “fixed” genotypes were set to 10⁻⁷ as explained above.

This model considers only one trait, the maximum growth rate at a given temperature. However, two other traits were considered. First, the mortality rate μ_mort is considered to follow a power law scaling increasing as a function of temperature (Brown et al., 2004; McCoy and Gillooly, 2008; and, specifically for phytoplankton, Regaudie-de-Gioux and Duarte, 2013). However, the observations do not show any sudden drop in realized growth rate when the cultures were warmed from 15 to 26.3°C, suggesting that this scaling is inappropriate for short term changes and more appropriate for asymptotic “steady state” conditions. Second, cell size is also a function of temperature, but calculations based on the experimental results suggest that it is a minor effect. Hence, neither mortality nor cell size are considered further.

The other important environmental variable in the study of Schlüter et al. (2014) is pCO₂. While only the 400 μatm case is considered here, it is likely that higher concentrations of CO₂ result in a reduction in the height of the “fitness window,” in analogy with the “thermal window” concept for animals (Pörtner and Farrell, 2008; Denman et al., 2011).

For the first 3 years, five replicates of the coccolithophore Emiliania huxleyi were grown in culture at 15°C for each of three pCO₂ levels: 400, 1100, and 2200 μatm. At the end of 3 years, the temperature of the cultures was raised in intervals of 1°C d⁻¹ to a final temperature of 26.3°C at which they were grown for an additional 1 year. Model simulations of the laboratory experiments with pCO₂ concentration of 400 μatm were performed as follows:

Based on arguments (Orr, 1998, 2005) critical of the concept that adaptive evolution proceeds according “micromutationism” or infinitesimal mutations as first postulated by Fisher (1930), the first set of simulations (for 3 years at 15°C) allowed for the random mutations to obey a flat pdf over the growth rate from 0 to 1.16 d⁻¹, the approximate observed growth rate μ_mean at 15°C (Figure 1). There were 58 equally wide genotypes over this space (numbered from 0 to 57). Initially they were all zero except for genotype #57, representing the measured mean growth rate 1.15 d⁻¹, shown as the solid diamond in Figure 1. The initial biomass of genotype #57 in each simulation is shown as the height of the pale yellow bar in Figure 2A.

In the simulations at 15°C, genotype #57 outcompeted all other mutant genotypes, all of which had lower growth rates (and hence lower fitness) than genotype #57. [We take the fitness of a mutant genotype m relative to the parent genotype p to be the ratio of their realized growth rates: W_mp = μ_m/μ_p, and of the relative fitness of mutant genotype i relative to mutant genotype j to be W_ij = μ_i/μ_j (Lenski et al., 1991; Schlüter et al., 2014). If W_mp (W_ij) is greater than 1, then genotype m (i) has a higher fitness than genotype p (j).] Therefore, after 3 years, in all five replicate simulations the most abundant genotype was #57 (Figure 2B), with its realized growth rate of 1.15 d⁻¹. The horizontal dashed line represents 10⁻⁵ of the total biomass, i.e., 1 cell. So by the end of the simulation, almost all genotypes except that with the highest realized growth rate (i.e., genotype #57) had biomass less than 1 cell and were effectively extinct. In the particular simulation shown in Figure 2B (one of 5 replicates) genotype #55 also had barely more than 1 cell, but it was a very recent mutation and its biomass would have quickly dropped below 1 cell, as can be seen in Figure 3. Figure 3A shows the magnitude of all the mutations over the 3 years: they are randomly distributed evenly over all genotypes between #0 and #57.

Figure 3B shows the time paths of the logarithm of the biomass of genotypes 20, 55, 56, and 57. Whenever there is a mutation to genotype # 20 (blue dashed line), it quickly dies out because its fitness relative to the parent genotype #57 is small, W_{20 57} ~ 0.36. Genotype #55 (green) dies out more slowly, W_{55 57} ~ 0.97, and genotype #56 dies out even more slowly, W_{56 57} ~ 0.98. The solid red line shows the total biomass of genotype #57, consisting mostly of the original genotype (pale yellow portion in Figure 2B) plus mutations to that genotype (lower blue portion), which is also shown by the solid black line in Figure 3B. Note that the vertical scale in Figures 2B, 3B is the logarithm of biomass, so actually only 0.016% of the biomass in genotype #57 consisted of mutant cells after 3 years (all 5 replicate simulations had a similar fraction of biomass in genotype #57 ~0.02%. These cells have growth rates equal to the original clone, but they may have other genes that differ from the original clone, as pointed out by Schlüter et al. (2014). However, the simulations in the next section started with 10⁵ cells of what is assumed to be the original genotype #57, which were then warmed instantaneously (in the model) to a temperature of 26.3°C.

In the simulations described here, there are 92 genotypes (numbered 0 to 91) between a growth rate of 0 d⁻¹ and the maximum growth rate μ_max at 26.3°C (solid red line in Figure 1), with each genotype interval being 0.02 d⁻¹ wide as before. Thus, the highest genotype #91 is centered at 1.83 d⁻¹ with its upper limit being μ_max(26.3°C) = 1.84 d⁻¹. Similarly, μ_max(15°C) = 1.15 d⁻¹ occurs at the center of genotype #57.

Five replicate simulations at 26.3°C with a flat pdf of the magnitudes of random mutations all quickly ended up with mutant genotypes with the highest growth rate or relative fitness eventually dominating the culture. Figure 4 tracks the time history of the biomass of the four highest genotypes: 88, 89, 90, and 91 for one of the 5 simulations. A mutation to #89 occurred first (on day 31), followed by one to #88 (on day 59), then one to #91 (on day 72), and lastly one to #90 (on day 219), with the four jumps afterwards indicating subsequent mutations to this genotype. These are not the only mutations to these genotypes then, just the first ones during this simulation. Note that #88 grew more slowly than #89 (because the fitness of #89 was relatively higher, but #91 outpaced them both, for the same reason. They all outcompeted the original clone #57 because of their significantly higher relative fitness. This coexistence of four mutant clones is an example of clonal interference observed in asexual populations (e.g., Muller, 1932; Gerrish and Lenski, 1998; Imhof and Schlötterer, 2001).

Given enough generations, the clone with the highest fitness will dominate the culture. Reducing the mutations from once a day (slightly greater than one generation) to every other day did not affect the end result, only marginally the rate of getting there. Figure 5 shows time series of the growth rate for the five replicate simulations, each with a new “seed” for the random mutations. The timing and nature of the increase in growth rate also depends on how soon mutations at the higher genotypes occur. Regardless, in each of these simulations, there is no significant increase in growth rate in the first ~20 days, then it rapidly increases eventually to the growth rate of the highest mutant genotype. This behavior does not follow the linear increase in μ_mean from 1.15 d⁻¹ at the start of the year to 1.33 d⁻¹ (with no discernable initial lag) that was observed in the laboratory (Schlüter et al., 2014). Moreover, the asymptotic mean growth rate after 1 year in these simulations is much higher than the eventual observed mean growth rate in the laboratory experiments (shown by the blue dashed line).

In the simulations with a flat pdf for random mutations, the invariable domination by mutations near the maximum possible growth rate μ_max suggests that in the laboratory experiments mutations with infinitesimal magnitudes might have been more likely (as suggested by Fisher, 1930), rather than mutations with “large” magnitudes being as likely as infinitesimal mutations (as argued by Orr, 1998, 2005, and others).

In these simulations, the magnitudes of mutations follow a Gaussian normal random pdf about the parent genotype. The width of the Gaussian normal pdf of mutation magnitude about the genotype undergoing a mutation is given by s (measured in genotype intervals of 0.02 d⁻¹). Hence, for s = 1 the distance from the center of the parent genotype interval to the outsides of the two adjacent genotypes is ±1.5s, i.e., the probability of a mutation occurring in the parent genotype or the adjoining genotypes is 86.6%. For s = 2 the distance from the center of the parent genotype interval to the outsides of the two adjacent genotypes is ±0.75s, i.e., the probability of a mutation occurring in the parent genotype or the adjoining genotypes is reduced to 54.7%.

Three sets of five replicate simulations were performed, with s = 1, 2, and 3. For s = 1, the growth rates after 1 year were well below those observed and these results are not discussed further. In Figure 6A, for s = 2, we show the distribution of biomass among genotypes after 1 year on a logarithmic scale for the replicate with the median mean growth rate. Although there appear to be about 20 genotypes in play, a linear plot of these genotypes (Figure 6B) shows that only about 4 or 5 genotypes account for almost all (at least 99%) of the biomass. We have therefore allowed only the 10 highest biomass mutant genotypes to undergo mutations themselves, with a probability inversely proportional to their relative biomass. Thus, from Figure 6B, on the next day, a mutation would be most likely be from #63, next most likely from #61, next most likely from #64 and so on, with the magnitude of each mutation determined randomly from a Gaussian normal pdf with s = 2. Figure 7A shows, for the same replicate simulation as in Figure 6, the time evolution of the biomasses of the original clone (#57, red) as well as the five mutant genotypes with the most biomass after 1 year: 61, 63, 64, 65, and 66, again demonstrating clonal interference, especially with #63 outcompeting #61, and with the higher fitness genotypes 64, 65, and 66 increasing each at a greater rate. Figure 7B shows the realized growth rate for the same simulation (heavy solid line), the mutations each day (jagged light solid line), and ± the standard deviation of the distribution of the biomasses of the active genotypes (see Figure 6A), along the trait or genotype axes (light dashed lines)

In this set of five replicate simulations with s = 2, the mean growth rate increased from 1.15 to 1.27 d⁻¹ after 1 year. In addition, the growth rate did not begin to increase until after ~110 days, roughly 125 generations (Figure 8). In the laboratory study (Schlüter et al., 2014), the realized growth rate increased, more or less linearly without discernible lag, to a value (from their fitted line) of 1.33 d⁻¹ at the end of 1 year. Thus, the increase in realized growth rate in these simulations was too small and the initial lag of ~110 days was unrealistic.

Another set of five simulations was performed with s = 3. Figure 9 shows the time path of the mean growth rate for each of the five replicate simulations (light lines), the overall mean growth rate (brown solid line) and the fitted linear increase in growth rate from the experiments (blue dashed line, Schlüter et al., 2014). After 1 year, the mean growth rate was 1.47 d⁻¹ (range 1.39–1.69 d⁻¹), exceeding that in the laboratory experiment. These simulations tend to exhibit lags of ~70–130 days followed by plateaus, especially the simulation with the highest growth rate after 1 year. In that simulation, after ~70 days the growth rate increased steeply to a plateau ~1.33 d⁻¹ (genotype #66), then at ~200 days it increased again to a plateau ~1.55 d⁻¹ (genotype #77) followed by another steep rise starting at ~320 days. Figure 10A shows the distribution of genotypes after 1 year for this simulation: here there were 12 genotypes with biomasses greater than 1 cell. Figure 10B tracks the time history of the simulation for genotypes 66, 77, 81, 82, and 85. The first large magnitude mutation from the parent genotype #57 to genotype #66 (on day 15) is rare (3s) but possible for a Gaussian normal pdf with s = 3; the second from #66 to #77 (on day 149) is also rare but possible. Figure 10B again shows clonal interference, as mutants with higher relative fitness eventually out-compete less fit mutants (cf. Figures 8, 9 in Gerrish and Lenski, 1998).

To summarize the results of this set of simulations, the overall mean growth rate increased roughly linearly after about 70 days, but again among the replicate simulations there was a lag of ~70–130 days before the growth rate(s) start to increase, contrary to the experimental data (Schlüter et al., 2014).

Initially, the model was set up with 58 genotypes (each 0.02 d⁻¹ wide) spanning the range of possible growth rates at 15°C from 0 to ~1.16 d⁻¹ (Figure 2 and Fielding, 2013). Mutations were equally probable (a “flat” pdf) to all 58 genotypes (#0 to #57) under the assumption that the original genotype #57 was the genotype with the highest possible growth rate and hence the maximum relative fitness. So only mutations to that same genotype survived or were “fixed,” while mutations to other genotypes (#0 to #56), all with lower relative fitness, became extinct (or failed to be “fixed”), as shown in Figures 2, 3. At the end of the 3 years, in all five replicate simulations, mutations to genotype #57 contributed ~0.02% to the total biomass shown in that genotype. Although these cells had the same growth rate as the original clones, they presumably possessed other genes that apparently did not affect their growth rate at 15°C.

Simulations at 26.3°C had 92 possible genotypes, each of width 0.02 d⁻¹, spanning the growth rates between 0 and 1.84 d⁻¹, the maximum possible growth rate for this clone (the value of the red line in Figure 1 at 26.3°C). The initial genotype was still #57, assumed to be the single genotype existing after 3 years of growing at 15°C, ignoring the ~0.02% of the cells that were mutants in the simulations but with the same growth rate at 15°C. All simulations exhibited a time lag after the temperature shift from 15°C before there was any significant increase in the realized growth rate, which was not observed in the laboratory experiments of Schlüter et al. (2014).

For a flat pdf of random mutations, the genotype created from the largest magnitude favorable mutation eventually dominated and replaced the original genotype. Usually, after 1 year the dominant genotype is #90 or #91 (Figures 4, 5). Figure 5 shows that there is typically a lag of ~15–30 days (~17–34 generations) before the realized growth rate starts to increase significantly. Then very quickly (over ~20 days) the growth rate climbs rapidly to that for the highest mutant genotype, where it remains for the rest of the year unless there is a subsequent mutation to a higher genotype. This behavior, two plateaus separated by an abrupt increase from the first to the second, is completely inconsistent with, and with a much larger final growth rate than, the final fitted growth rate from the laboratory experiments.

For the flat pdf, each mutation has a 62% chance of having a relative fitness less than that of the original clone (57 of 92 possible genotypes). If the rate of mutations were to decrease by a factor of 10, then the lag time would be 10 times longer, but the time for a given favorable mutation to increase would be the same because it is a function of the generation time or growth rate.

It was concluded from these simulations that an initial lag followed by an abrupt increase in growth rate results from allowing mutations of the largest magnitude to have the same probability as mutations of the smallest magnitude. If small magnitude mutations were to have a much higher probability than large magnitude mutations, then possibly the transition to higher grow rates after increasing the temperature to 26.3°C would be more gradual and continuous. To allow for mutations of very small magnitude, in subsequent simulations the magnitude of each mutation along the trait axis was chosen randomly from a Gaussian normal pdf, N(μ_max, s) in statistical notation, centered on the growth rate of the parent genotype, with a width along the trait axis scaled by the standard deviation s (in genotype intervals).

The first set of five simulations shown here, with s = 2 genotypes wide, generated a much too small increase in mean growth rate μ_mean, reaching only 1.27 d⁻¹ compared with the fitted value of 1.33 d⁻¹ from the laboratory experiments. Again, the simulations show a long lag of ~110 days before any significant increase occurred (Figure 8). A second set of five simulations with s = 3 was then performed (Figure 9). A larger number of genotypes were generated from mutations, with many of them still viable (1 cell or more) at the end of 1 year (Figure 10). The mean growth rate of the five replicate simulations increased more or less linearly for the latter two thirds of the year, with a slope roughly double that observed. But there still remained a lag of ~70 days before the mean growth rate started to increase. Overall, for a larger s, i.e., larger magnitude mutations, the lag time was shorter. For s = 1, 2, and 3, the lag time was respectively ~150, ~120, and ~70 days. Clearly, adaptive evolution by genetic mutations, modeled in the manner described here, cannot alone explain the laboratory results (Schlüter et al., 2014) because all simulations were characterized by initial lags upon warming to 26.3°C.

A possible explanation for the immediate continuous increase in observed growth rate after increasing the temperature from 15 to 26.3°C is that it was a “plastic” response of the cells to their changing environment. The idea of plasticity “buying time” for genetic adaptation to take place is central to the concept of “plastic rescue” avoiding extinction (e.g., Lande, 2009; Chevin et al., 2010; Kopp and Matuszewski, 2014). There are many definitions of plasticity: Whitman and Agrawal (2009) list 11, but perhaps their simplest is “Phenotypic plasticity” is “the capacity of a single genotype to exhibit variable phenotypes in different environments ….” According to Reusch (2014), reviewing evidence of plasticity in marine animals and plants: “Phenotypic plasticity broadly defines the adjustment of phenotypic values of genotypes depending on the environment, without genetic changes.”

Most studies of plasticity compare different phenotypes of the same genotype in different environments. More relevant in our case would be studies that consider the temporal change in phenotypic properties of a given genotype in response to a change (continuous or abrupt) in its environment, but such studies are rare. And different species generally have different plastic responses to the same change in environment (different “reaction norms,” (e.g., Pigliucci, 2005; Whitman and Agrawal, 2009; Reusch, 2014). There is no clear information on how plasticity might be modeled in this case, especially what ultimately limits the rate and magnitude of plastic response. Currently, both the energetic costs and limits of plasticity are research questions of considerable interest (e.g., DeWitt et al., 1998; Pigliucci, 2005).

Given that the generation time of E. huxleyi in the laboratory experiments was less than 1 day, a plastic response would have to be “transgenerational” or heritable, for which there is mounting evidence in animals (Munday, 2014; Walsh et al., 2014), in clonal plants (e.g., Latzel and Klimešová, 2010), and in phytoplankton (Schaum et al., 2013; Schaum and Collins, 2014). For asexual reproduction the hypothesized mechanism is “epigenetic inheritance” whereby an environmental change causes genes to be expressed, which continue to be expressed in succeeding generations if the environmental change continues (Latzel and Klimešová, 2010; Schaum et al., 2013; Schaum and Collins, 2014; van Oppen et al., 2015). In the laboratory experiments of Schlüter et al. (2014), it is assumed that the change in temperature from 15 to 26.3°C caused the expression of an existing gene in essentially all cells in the culture in the first and succeeding generations which mediated a slow increase in growth rate, without the lag that would result from the favorable mutation of a single cell dividing sufficiently often to start to compete with the original population.

Here we assume that the plasticity is heritable, and model it as a first order restoring function without lag: (2)dμi(t)dt=(μmax(i) − μi(t))/T where μ_i(t) is the realized growth rate and μ_max(i) is the maximum growth rate, both for the initial genotype i, according to the power law fit (Fielding, 2013) through the original growth rate at 15°C (solid red line in Figure 1). T (or τ) is the “e-folding” time for the response. The assumed mechanism is that because of energy costs the plastic response of a given genotype i cannot exceed its μ_max(i) value for 15°C, even though it is now at a higher temperature. The rate at which μ_i(t) approaches μ_max(i) decreases as it gets closer, the logic being that the larger the plastic response, the more energy that is required.

The solution to the ordinary differential Equation (2) is standard: (3)μi(t)=(μi(0)−μmax(i))e−t/T+μmax(i)

This function has the form shown Figure 11, where the dashed curve shows the full function to its asymptote. The solid curve has a T value set to 281 d, so that the initial slope of Equation (3) at small t is equal to the slope of the linear fit to the observations over 1 year (Figure 1A in Schlüter et al., 2014). For the initial genotype, μ_max = 1.29 d⁻¹, so that the observed (fitted) growth rate after 1 year of 1.33 d⁻¹ could not be reached by means of plasticity alone. In fact for the value chosen for T, the plastic response at the end of 1 year would result in a mean growth rate of only 1.25 d⁻¹ (solid curve, Figure 11).

Three sets of five replicate simulations were then performed at 26.3°C with s = 2, 2.5, and 3, but with a plastic response of the original genotype. Now the growth rate of the initial genotype (initially μ_i(0) = 1.15 d⁻¹) increased with time according to Equation (3). As its growth rate increased, the biomass of the original genotype was added to any in the appropriate growth rate interval that had accumulated from mutations. Due to the plastic response, the original genotype had moved from #57 to #62 on the growth rate axis after 1 year. Figure 12A shows the results of 1 replicate simulation with s = 2.5; in all five replicate simulations with s = 2.5, the original genotype represented most of the biomass in that interval. Figure 12B shows the time evolution of the five other genotypes with the most biomass after 1 year: clonal interference between genotypes 63, 64, and 66. Genotype #66, because it had the highest relative fitness, outcompeted genotypes 63, and 64, even though they had mutated earlier. Figure 13A shows the time history of the growth rates for all five replicate simulations, the ensemble mean growth rate (solid brown line), and the (dashed blue) line fitted to the laboratory results. The plastic response eliminated the initial lag in the increase of the growth rate, consistent with the observations. The ensemble growth rate after 1 year was 1.34 d⁻¹ (slightly greater than that of the observations, 1.33 d⁻¹). In Figure 13B, simulated random measurement error has been added to the simulated overall mean growth rate, for comparison with the observations (Figure 1A in Schlüter et al., 2014). Thus, the plastic response, as formulated, removed the lag in response present in all previous simulations, giving time for mutations to new genotypes eventually to dominate the culture toward the end of the year. In addition to the maximum growth rate, the only other trait reported on by Schlüter et al. (2014) was cell diameter. The cell diameter was significantly smaller (~10%) at 26.3°C, but only for the cultures grown at the “ambient” level of pCO₂: 400 μatm.

Schlüter et al. (2014) maintained cultures of Emiliania huxleyi, all originating from the same single cell, at three different pCO₂ levels for 4 years, increasing the temperature from 15 to 26.3°C at the end of the third year. Although the growth rate at 15°C decreased with increasing pCO₂, the growth rate increase over the last year (at 26.3°C) was greater at successive higher pCO₂ levels, such that the growth rates at the end of the experiment were closer together than during the first 3 years, suggesting that at the higher temperature, the cells were affected less by CO₂ concentration. Without some information about energy costs of adaptation, it is not clear how to model either the effects of mutation (or of plasticity) in response to two simultaneous stressors.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.