Evolution in fluctuating environments is predicted to disfavor genetic canalization and instead select for alternative strategies, such as phenotypic plasticity or possibly bet-hedging, depending on the accuracy of environmental cues and type of fluctuations. While these two alternatives are often contrasted in theoretical studies, their evolution are seldom studied together in empirical work.

We used experimental evolution in the nematode worm Caenorhabditis remanei to simultaneously study the evolution of plasticity and bet-hedging in environments differing only in their temperature variability, where one regime is exposed to faster temperature cycles between 20°C and 25°C, with little autocorrelation between parent and offspring environment, while the other regime had slowly increasing temperature with high autocorrelation in temperature between parent and offspring. We exposed worms for 30 generations to either fluctuating or slowly increasing temperature, these two environments had the same average temperature over evolutionary time. After experimental evolution, we scored size at sexual maturity and fitness in full siblings reared in two different temperatures, optimal 20°C and mildly stressful 25°C.

The datasets contain data on adult size, total reproduction and individual fitness of full-sibs from the two evolutionary selection regimes, where each regime was replicated in 6 replicate populations.

Experimental evolution in the fluctuating environment resulted in the evolution of increased body size plasticity but not increased bet-hedging, compared to evolution in the slowly changing environment. Plasticity followed the temperature size rule as size decreased with increasing temperature and this plastic response was adaptive. In addition, we documented substantial standing genetic variation in body size, which represents a potential for further evolutionary change.

Experimental evolution

We used C. remanei nematode worms, strain SP8 which has been lab adapted at 20°C and subsequently exposed to 30 generations of experimental evolution in two regimes (Increased warming and Fast temperature cycles). The experimental evolution has been previously described in detail in Lind et al., (2020). Briefly, in the Increased warming experimental evolution regime, the temperature gradually raised from 20°C to 25°C, which is a novel and mildly stressful temperature. This gradual change over 30 generations represents an increase of 0.1°C every 2.13 days and results in a correlated parental and offspring environment. In the second regime, Fast temperature cycles, the temperature fluctuated every second generation between 20°C and 25°C, resulting in 14 temperature shifts but no exposure to the intermediate temperatures. The environmental change is deterministic (every second generation), but since parents and offspring would end up in either the same or in different temperature, it represents uncorrelated parental and offspring environment. The generation time in C. remanei is temperature dependent; 4 days long in 20°C and 3.4 days long in 25°C. Despite these differences, the average temperature and the total chronological time of experimental evolution were identical for both regimes, at 22.5°C and 110 days, respectively.

Each evolutionary regime consisted of six replicate populations. The populations were maintained on 92 mm Petri dishes poured with NGM agar in climate chambers set to 60% relative humidity. In order to prevent bacterial and fungal contamination, the agar and bacterial LB contained the antibiotics streptomycin and kanamycin, and the fungicide nystatin. The plates were seeded with 2 ml of an antibiotic-resistant OP50-1 (pUC4K) strain of E. coli (Stiernagle, 2006) that served as a source of food. Every 1-2 days, a piece of agar containing approximately 150 worms of mixed ages was cut and transferred to a new plate containing fresh bacteria. This resulted in populations with overlapping generations that were maintained in a constant exponential growth phase. After the experimental evolution, populations were expanded for two generations and frozen in -80°C for later revival and subsequent phenotypic assays.

Experimental set-up

Each replicate population of each of the two selection regimes was run in a separate block resulting in 12 experimental blocks in total. For logistical reasons, we focus on females, since they are responsible for population growth rate and their fitness is straightforward to measure.

Briefly, populations were revived from freezing and exposed to 25°C for 3 generations, to avoid any maternal effects associated with freezing. The third generation was split into eight families, each family consisting of one male and one female worm. From each family, we randomly picked eight offspring females (full siblings) and placed four females in 20°C and four in 25°C. Since our focus was evolution in females, their fitness was assessed by mating them with standardized males from the ancestral line. For the detailed description of the experimental set up, see supplementary figure 2.

Phenotypic assays

Daily reproduction

Female and male worms were transferred to a new plate every 24 hours, and viable offspring were counted two days later. The female worm was discarded after dying, or when reproduction ended. 

Body size

Worms in 20°C reach their peak size at day 4 of adulthood (Lind et al., 2016). The peak size in 25°C is on day 2 of adulthood, which was determined during pilot assays (supplementary figure 3). Photographs of worms were taken during their peak size using a Lumenera Infinity2-5C digital microscope camera mounted on a Leica M165C stereomicroscope. Body size was measured from the photographs using ImageJ 1.46r (https://imagej.nih.gov/ij/) as total cross-section area.

Individual fitness

We used the age-specific reproduction data to calculate rate-sensitive individual fitness λind for each individual (McGraw & Caswell, 1996), which is analogous to the intrinsic rate of population growth. Individual fitness was calculated by constructing a population projection matrix for each individual, and then calculating the dominant eigenvalue of this matrix, following McGraw & Caswell, (1996). Since we kept the population size and age structure constant during experimental evolution, individual fitness is the most appropriate fitness measure for this study (Mylius & Diekmann, 1995).

Thermal reaction norms

To test whether the degree of phenotypic plasticity has evolved, we used linear mixed-effect models to separately estimate the thermal reaction norms of body size and individual fitness, using the package lme4 (Bates et al., 2015) in R. The models included either body size (area) or individual fitness (λind) as response variables. The full model included three fixed effects: mean-standardized temperature as a covariate, the experimental evolution regime as a categorical factor, and their interaction. We expect this interaction to be significant if the degree of plasticity has evolved. Experimental line and dam identity were included as random effects. Significance of the fixed effects was evaluated using Wald χ2 tests. Pseudo-R2 values were calculated as the squared correlation coefficient between fitted values from the model and observed values.

Selection

To test if temperature responses in size are adaptive, we estimated the selection on body size and compared it to the observed temperature response. Selection on body size (area) was estimated using mixed-effect models in R with the package lme4 and individual fitness (λind) as the response variable. The full model included the following fixed effects: area, area2, temperature, experimental evolution regime and all interactions except for interactions involving both area and area2 together. Experimental line was included as random effect. Significance of fixed effects was evaluated using Wald χ2 tests.  From the full model, we obtained temperature-specific estimates of the slope and the squared term between λind and body size. For each temperature, the optimal size (i.e. the area that maximizes fitness) was calculated as: -b/(2×c), where b = the temperature-specific slope from the full model (i.e. the linear selection gradient, estimating the relationship between body area and λind) and c = the temperature-specific squared term from the full model (i.e. the quadratic selection gradient, estimating the relationship between area2 and λind). Confidence intervals of the temperature-specific optimal sizes were generated by bootstrapping, implemented in the boot package using 10 000 bootstrap replicates.

Within family coefficient of variance

To test whether the degree of diversifying bet-hedging has evolved, we tested whether the experimental evolution regimes differed in the mean within family variance within temperatures. For each family, we used the trait values of the offspring (within a temperature) to calculate the within family variance. To account for differences in trait means, we used within family means to mean-standardize the variance by calculating the within family coefficient of variance (CV): CV = σ/μ, where σ = within family standard deviation, and μ = within family mean. An ANOVA was used within each temperature to test if the evolution regimes differed in their mean within family CV.

Since it is more difficult to detect differences in variances than differences in means, we also performed power calculations on our ability to detect whether within-family CVs differ between the selective regimes. Balanced one-way ANOVA power calculations were performed to estimate the effect sizes possible to detect with power ranging from 0.70 to 0.95. Effect sizes, η2, were obtained for our sample size of N = 48 per selection regime and a significance level of 0.05. η2 is calculated as the sum of squared explained by the treatment (here: selection regime) divided by the total sum of squares and has an equivalent interpretation as an R2.

Genetic variance and correlations

For body size, genetic variance and genetic correlations across temperature were estimated using animal models in the package MCMCglmm. Univariate models were used to estimate genetic variance, whereas bivariate models were used to estimate genetic correlations, both models using Gaussian family for trait distribution. An inverse Wishart prior with parameters V = 1 and nu = 0.02 were used in both univariate and bivariate models. Pedigree data linking offspring to parents, based on full-sib relationships, was included in the models. Convergence of the models were ensured by evaluating diagnostic plots of posterior distributions, using the convergence diagnostic half-width test by Heidelberger and Welch (1983), and by ensuring that the autocorrelation between MCMC samples was close to zero.

For univariate models, body area was used as response variable. Temperature (categorical), experimental evolution regime, and their interaction, were included as fixed effects. Genetic variance (VG), variance due to differences between experimental lines, and residual variances were estimated separately as random effects in the full model for each temperature-by-evolution regime combination. The full model ran for 4.2×106 MCMC iterations, 0.2×106 samples were discarded as burnin, and the thinning interval was 4000, resulting in a sample size of 1000 MCMC-samples. Reduced models, subset by temperature-by-evolution regime combination, were used to assess statistical significance of VG, by comparing deviance information criterion (DIC) of models with versus without genetic variance included.

Broad sense heritability (H2 = VG/VP, where VP = total phenotypic variance after accounting for variance due to experimental line effects) and broad sense evolvability (I2 = VG/mean2, Hansen et al., 2003, 2011) were used to estimate the population’s evolutionary potential of body size. This was estimated for each temperature-by-evolution regime combination. Evolvability measures the expected percentage change in a trait per generation under unit strength of selection. Compared to heritability, evolvability is independent from the environmental variance and represents a measure of the evolutionary potential that is comparable across traits, populations, and species when applied to traits with a natural zero and which are strictly positive (Hansen et al., 2011).

Genetic correlations of body size were estimated using bivariate animal models in MCMCglmm. Genetic correlations was estimated separately for the two regimes. Body size was the response variable and was treated as two traits (size at 20°C and at 25°C). Random effects included genetic covariance between the temperatures, whereas VG, variance due to differences between experimental lines, and residual variances were estimated separately for each temperature. The full models ran for 2.05×106 MCMC iterations (burnin: 0.05×106 samples, thinning interval: 2000 samples), resulting in a sample size of 1000 MCMC-samples. Reduced models without genetic covariance were used to access the statistical significance of the genetic covariance, by comparing the DIC of models with versus without genetic covariance included. The genetic correlation of body size across temperatures was calculated by dividing the genetic covariance by the product of the genetic standard deviation of the two temperatures. This was done on the posterior distributions, in order to carry the error forwards in the analyses.

To compare posterior distributions of H2, I2, and genetic correlations across temperatures and selection regimes, we calculated, within each MCMC sample, the pairwise differences in these measures and checked if the posterior distributions of these differences had a 95% credibility interval that included zero. Pairwise comparisons of distributions were performed between evolution regimes within temperature, and between temperatures within evolution regimes.