Eco-evolutionary dynamics in response to selection on life-history

Understanding the consequences of environmental change on ecological and evolutionary dynamics is inherently problematic because of the complex interplay between them. Using invertebrates in microcosms, we characterise phenotypic, population and evolutionary dynamics before, during and after exposure to a novel environment and harvesting over 20 generations. We demonstrate an evolved change in life-history traits (the age- and size-at-maturity, and survival to maturity) in response to selection caused by environmental change (wild to laboratory) and to harvesting (juvenile or adult). Life-history evolution, which drives changes in population growth rate and thus population dynamics, includes an increase in age-to-maturity of 76% (from 12.5 to 22 days) in the unharvested populations as they adapt to the new environment. Evolutionary responses to harvesting are outweighed by the response to environmental change (∼ 1.4 vs. 4% change in age-at-maturity per generation). The adaptive response to environmental change converts a negative population growth trajectory into a positive one: an example of evolutionary rescue.


Introduction to Supplementary Information
This supplementary information (SI) has been provided to give more detail about our study system, methods and analysis. Further analyses include plots of juvenile/adult ratios to see changing demography, plots to show how repeatable replicate time series were and a full statistical model representing the direct and delayed density dependent effects on survival and recruitment. A full explanation of the multivariate approach we took to jointly model age and size at maturity and how this model was used to predict co-dependent age and size at maturity for each assay time point is provided.
We present here in the Supplementary Information two plots to supplement our analysis of molecular variance that describes how a number of our AFLP markers were probably under directional selection and that we found significant levels of population differentiation in quantitative traits. Finally, we present a full description of the methods used to estimate rate of change of population size from life-history that allows us to consider how well information from the life-history assays can predict changes in population growth.

Model System
The soil mite, Sancassania berlesei (Michael), is used for the experimentally harvested model system. This species has a plastic life-history in that growth rates, survival and fecundity vary with resource allocation (Beckerman et al. 2003). Furthermore, juveniles compete with adults, and are morphologically similar, and the juvenile phase is typically short relative to the adult phase (Plaistow et al. 2004). Juveniles can be identified as distinct from adults due to the relative size and shape of different body parts. Adults are competitively dominant over juveniles as they can move to new food patches faster and pull juveniles off existing food patches. But juveniles can feed amongst detritus and underneath adults, and it should be remembered that there are many more juveniles than adults. If small juveniles are culled, maturation rates into the adult class increase suggesting the juveniles are food limited and adults are better competitors (Cameron & Benton 2004). Under equilibrium conditions the population has noisy equilibrium like dynamics. Previously we have shown that the lifehistory of this organism can respond to changes in current and previous environments and there are high levels of genetic variation in heritable life-history traits (Beckerman et al. 2002;Beckerman et al. 2006;Plaistow et al. 2006). And we have shown that this system lends itself well to investigating the population dynamics consequences of harvesting (Cameron & Benton 2004) and with a short generation time (10-50 days) it is suitable for asking evolutionary questions over contemporary timescales. General laboratory and rearing techniques can be found in work we have published previously (Benton et al. 2001).

Further Analysis and Supporting Graphs (a) Population dynamics SI
The dynamics can be summarised as having two main features, an initial decline in mean abundance of all stages that takes approximately 40 weeks (total population, the minimum is different for each stage) and then a second phase where the mean and variance in density increased ( Figure S1). In the unharvested control the mean density levels off at the end of the experiment around week 70. The increase in variance and the levelling off is particularly evident in the juvenile and juvenile:adult ratio plots (Figure S1 b).
Using an a priori list of candidate models that include recruitment to and survival within each life-history stage, harvesting, time and their interactions models that captured the logged population dynamics of the mites best were chosen based on their weighted AIC values such that we did not assume there was one best model (Table S1). Log transformations were undertaken as a precaution to ensure normality and reduce patterning in the residuals. The best set of models with ΔAIC <6 is presented to estimate the approximate effect of each harvesting treatment on recruitment, survival and population growth of each stage; not as a means to predict population dynamics (Table S2).
The best statistical model of the population dynamics of each stage, including direct and delayed density dependence (where stage densities are influenced by current and past mite stage densities), shows significant changes in the dynamics of all stages over time (Table S1 and S2). However, the sign and magnitude of these changes differ between stages (stage*week interactions, for example: The average coefficients from the set of selected models show that juvenile survival increased over time (Survival:Week coef = 0.0019) whereas adult survival did not change over time (Survival:Week coef = 0.00). This is what leads to the changes in the relationship between juveniles and adults over time ( Figure S1b). The increasing amplitude of the cycles in in all stages under all harvesting treatments occurs around week 60 coinciding with the period where juvenile development periods approach the 21 day mark (see LH assay timepoint 3 results, week 63). 21 days represents when mite populations would receive their 3 rd harvesting event per generation.
Imposition of harvesting leads to immediate changes in age structure: harvesting adults significantly reduces adult densities via reducing adult survival (i.e. direct lagged effect of harvest, Adult dynamics model: Survival coefficient for adult harvested = -0.11 log scale = 11%) and reduces the density of potential recruits via indirectly reducing juvenile densities (Juvenile dynamics model: Survival coefficient for adult harvested = -0.17 log scale = 16%) (Table S2). So despite removing 40% of adults each week, this change in adult numbers is partially compensated for via a reduction in juvenile densities. This occurs due to a significant negative relationship between juvenile densities and adult recruitment (such that large juvenile densities result in reduces adult recruitment (model averaged coefficient -0.09 (log scale)) partially offsets these direct and indirect negative effects by reducing competition such that recruitment increases leading to rapid replacement of harvested adults.
Post harvesting in unharvested and adult harvested populations, adult numbers become less variable but remain high whereas when released from juvenile harvesting adult numbers increase but their amplitude remains relatively unchanged.
Inter-replicate variability is small ( Figure S1d) and where data has been missing due to staff absence at weeks 20, 56, 83 and 97; the data were spline interpolated over weeks 19-21 but due to the nature of interpolations poorly predicting peaks and troughs the average of the two peaks or troughs adjacent to the missing value were used on weeks 56, 83 and 97. Figure S1. Panel plot showing the average across 6 replicates of the pre-harvest weekly census for (a) juveniles,(b) eggs, (c) Mean total population size fitted using a prediction from a General Additive Model (4 degrees of freedom selected by AIC) and (d) Unharvested adult number data.

(b) Multivariate Analysis of Variance of Harvest treatments on Post-Harvesting Phenotype
To determine the role of harvesting on the final phenotypes in either high or low food juvenile growth environments we used MANOVA to model log age-at-maturity and log sizeat-maturity as a function of the factorial harvest treatment (har =Unharvested, Juvenile or Adult Harvested) and food availability (Low or High). In addition we included a priori defined measures of density to control for differences between replicates in competitionrelated food availability as numbers change. These include total survival (surv), the summed daily densities (i.e. weighted density =densWt) and the median density (densmed) as a consequence of evolved changes in growth rates. Each of these juvenile density covariates can influence adult phenotype in a unique way and could be a consequence of any evolved changes in investment in life-history as well as differences between genotypes from each family selected for the life-history assay (see life-history assay methods). To control, as would be done in the analysis of a randomised block design, for the random differences between populations (pop) within treatment (har) and random difference between families (fam) within populations these nested factorial covariates were included as fixed effects (whether significant or not). To generate Figure 2e the covariates for survival, weighted density and median density were standardised with the same mean values across all treatments at the last timepoint. The interaction between harvest treatments and food availability on the phenotype is significant (Table S1) and therefore our full model is our minimal model:    Figure S2. Panel plot summarising goodness-of-fit of the minimal manova of the post-harvesting phenotype at the end of the experiment where we compare model fits (open circles) with observed family means of age and size at maturity (closed circles) for (a) unharvested (mean difference between paired observed and fits is zero; age/size: t= -0.1177/0.0069, df=27, p>0.90), (b) juvenile harvested (age/size: t= -0.1153/-0.1198, df=27, p>0.90) and (c) adult harvested treatments (age/size: t= -0.1368/-0.0168, df=27, p>0.87) when each density covariate are kept at observed replicate means for each food availability. Model residuals for age and size are distributed evenly (d) and a plot of residuals vs fitted values suggests the distribution is relatively constant with a change in mean.

Age residuals
Size residuals Table S5. The table displays the age and size at maturity of female mites from each life-history assay time point either as measured family means or MANOVA predicts controlling for differences in survival between tubes within assays (see explanation at 2a above) or . A separate model, as described above, was built for each assay time point. For assay time point 1, where covariate information for the median and weighted density are unavailable, mid points between assay time point 0 and 2 are used to make best estimates. We have shown our MANOVA is a good predictor of the original data (see Figure S2). See Figure 2 to understand variation between families at timepoint 4. WT ´= wild type, U = unharvested, J = juvenile harvested, A = adult harvested.

Predictions using tube means
Predictions using treatment means

(d) Calculation of R 0 values from the life-history data
To understand the feedback from adaptation to the dominant local environmental conditions, strong resource competition/ starvation, on population dynamics we calculated the basic reproductive rate (R 0 ) based on the family mean life-history trait values for age-at-maturity, size-at-maturity and daily survival rates per harvest treatment from our low food phenotype full sib females (See life-history assay methods in SI and Box 1). Low food phenotypes are chosen as they better represent the food limited conditions experienced in the microcosm tubes. Estimates of fecundity were obtained from the regressions in the main text based on the results of an earlier experiment (Plaistow et al. 2006). We worked with peak fecundity produced between, and including, day 2 to 7 post maturity and not fecundity per lifetime. The effects of size on fecundity under low food conditions is not as pronounced as under high food conditions and is dependent on adult survival. The effects of size on fecundity measured up to 8 days post maturation is neutral, but peak daily fecundity in days 2-4 post maturation is strongly linked to increased size. Therefore adult harvested phenotypes are under strong selection to invest in increased body size where unharvested phenotypes are not. Fecundity estimates for unharvested and juvenile used the neutral size-fecundity relationship whereas for adult harvested phenotypes the relationship between size and fecundity 2-4 days post maturation was used.
Where R 0 = (l x *m x ) Where R 0 is the sum of l x , the chance of an individual surviving to age x, times m x , the number of offspring produced during the time from age x-1 to x. Due to the overlapping generations we correct R 0 as follows: R 0 =exp((ln R 0 )/T c .) R 0 is corrected by T c , the average cohort lifespan, which is the average time from the birth of an individual to the birth of one of its offspring. As the average time from maturity to peak fecundity is 3 days we took T c to be the average time to sexual maturity plus 3 days.
We created a function to calculate R 0 in R2.14.1 and using 1000 resamples from the family mean life-history data combined with resamples of fecundity coefficients (selected from a normal distribution with the same mean and sd as our data) we estimated bootstrapped mean and confidence intervals of R 0 per assay time point per harvest treatment.
This analysis was undertaken to show how a suboptimal phenotype with low R 0 has changed during the study and changed a declining population into an increasing one. But also to show how differences in the life-history over time and between treatments manifest themselves in differences in population dynamics.
We compared our estimates of R 0 from life history data with the average population growth rates across all replicates per treatment. We expect high correlation between estimates of R 0 and population growth rates as R 0 per generation, corrected for generation time or not, corresponds to other measures of population growth rate such as λ (Caswell 2000, p128-129).