A fitness trade-off between seasons causes multigenerational cycles in phenotype and population size

Although seasonality is widespread and can cause fluctuations in the intensity and direction of natural selection, we have little information about the consequences of seasonal fitness trade-offs for population dynamics. Here we exposed populations of Drosophila melanogaster to repeated seasonal changes in resources across 58 generations and used experimental and mathematical approaches to investigate how viability selection on body size in the non-breeding season could affect demography. We show that opposing seasonal episodes of natural selection on body size interacted with both direct and delayed density dependence to cause populations to undergo predictable multigenerational density cycles. Our results provide evidence that seasonality can set the conditions for life-history trade-offs and density dependence, which can, in turn, interact to cause multigenerational population cycles. DOI: http://dx.doi.org/10.7554/eLife.18770.001


Introduction
In many organisms, reproduction is confined to seasonal fluctuation in periods of high resource, in which both fecundity (reproduction) and viability (survival) selection can occur, and periods of low resources, when reproduction stops and natural selection occurs only through viability selection. Consequently, sequential episodes of reproduction and survival caused by seasonality could be a major source of fluctuations in the strength and direction of natural selection (Darwin, 1859;Lack, 1954;Fretwell, 1972;Schluter et al., 1991;Bell, 2010;Bergland et al., 2014), giving rise to classic life-history trade-offs (Lack, 1947;Roff, 1992;Stearns, 1992;Garland, 2014). More specifically, traits that confer a fecundity advantage, but which are associated with a survival cost, will experience natural selection in one season that is opposed by selection in the subsequent season (Levins, 1968;Michod, 2006;Bell, 2010;Bergland et al., 2014). The sequential rather than simultaneous nature of trade-offs driven by seasonality could have important consequence for the trait distribution within and across generations (Levins, 1968;Grafen, 1988;Michod, 2006) and population dynamics (Ozgul et al., 2010).
One way by which life history trade-offs might arise from seasonal variation in resources is via body size (Ozgul et al., 2010(Ozgul et al., , 2014. Large individuals usually have higher fecundity (Mueller and Joshi, 2000;Schulte-Hostedde and Millar, 2004), but they could also have lower survival during the non-breeding season when resources are scarce (Stockhoff, 1991;Reznick et al., 2000;Munch et al., 2003;Monaghan, 2008). In addition, large individuals might take more time to grow and require more resources for maintenance (Munch et al., 2003), which could negatively impact their survival probability (Kingsolver and Huey, 2008).
This association between fitness and body size in seasonal environments could have important consequences for population dynamics, particularly when selection on body size is density-dependent (Mueller, 1997;Sinervo et al., 2000;Travis et al., 2013). Differences in the selective advantage of body size across seasons could also shed light on the long-standing question about why population densities of many species fluctuate periodically over time (Elton and Nicholson, 1942;Kendall et al., 1999;McCauley et al., 2008;Yan et al., 2013). For example, when body size is positively related to fecundity, but small individuals survive better in the non-breeding season (Stockhoff, 1991;Munch et al., 2003;Monaghan, 2008;Betini et al., 2014), these opposing patterns of selection could cause population cycles if selection is density-dependent. Specifically, if smaller offspring have higher survival when density is high, then the population will be composed of smaller than average individuals with lower average fecundity in the following breeding season. This lower mean fecundity will reduce population growth rates even though larger individuals will be favoured through fecundity selection. As population size declines, the strength of density-dependent viability selection on body size will also decline, which could cause net selection to reverse and favour larger individuals due to the fecundity benefit of being larger. As large individuals increase in frequency, population size should also increase via an improvement of reproductive output, returning populations to high densities. Although changes in the intensity and direction of natural selection caused by environmental variation are widespread (Schluter et al., 1991;Bell, 2010;Thompson, 2013;Bergland et al., 2014), we have little information about whether opposing episodes of natural selection could arise from seasonality and indirectly affect population dynamics through the feedback loop between ecological (density dependence) and evolutionary (selection and evolution) processes (Chitty, 1960;Krebs, 1978;Hairston et al., 2005;Pelletier et al., 2009;Ozgul et al., 2010;Schoener, 2011).
Here, we investigated how a seasonal fitness trade-off related to body-size could affect changes in population size and body size over time using replicate populations of Drosophila melanogaster exposed to repeated changes in food resources. In addition to standard breeding conditions for Drosophila, we also created a 'non-breeding season' by manipulating the food medium to prevent eLife digest Many wild populations go through long cycles in abundance that span several generations. The traditional explanation for such "multigenerational" cycles is that they are driven by predator/prey relationships, the classic example being oscillations between the numbers of lynx and snowshoe hares.
Population cycles could also be driven by seasonal changes. For example, traits that help animals to produce large numbers of offspring during the breeding season may reduce the ability of the animal to survive the non-breeding season. Body size is one such trait. Large individuals tend to produce more offspring, but their larger body size means that they find it harder to survive when food is scarce. As a consequence, large individuals should have an advantage and be more common when the population size is low and there are enough resources for all individuals. However, small individuals should be more abundant when population size is high. This trade-off caused by seasonality could set the population in motion towards predictable, multigenerational cycles.
To test this idea, Betini et al. established populations of fruit flies that went through 'breeding' and 'non-breeding' seasons. This was achieved by periodically altering the flies' food to prevent the females from laying eggs (in the lab, fruit flies do not normally have non-breeding seasons). Over 58 generations, the number of flies in each population cycled between peaks of high and low numbers.
When the population contained relatively few flies, there was strong selection for large flies because they have high reproductive success. Hence, the population grew. When the population was large, meaning that the flies had to compete for a limited amount of food, there was strong selection for small flies because they are better able to survive on limited resources. However, small flies also produce fewer offspring on average, resulting in a decrease in population size. When the flies all had sufficient food during the non-breeding season, these regular cycles completely disappeared.
A major challenge will be to understand how common this phenomenon is in the wild. Virtually all organisms live in seasonal environments but whether they face strong trade-offs in the expression of traits is not well understood. This is primarily because of the difficulty in following individuals throughout the year. females from laying eggs during this period. Thus, in this system, breeding and survival were restricted to two sequential and distinct seasons (hereafter 'breeding' and 'non-breeding'). The number of days and amount of food in each season was determined so that both fecundity and nonbreeding survival were density-dependent (Betini et al., 2013a(Betini et al., , 2014(Betini et al., , 2015, which is an important feature of many populations. In Drosophila, as in many other species, the positive correlation between body size and fecundity is well known (Mueller and Joshi, 2000; Appendix 1) and we previously demonstrated that small individuals have higher survival during the non-breeding season when abundance is high (Betini et al. 2013a(Betini et al. , 2014. In addition, populations of D. melanogaster do not show evidence of multigenerational cycles (Mueller and Joshi, 2000), even when kept under the same conditions as our breeding season (i.e. 'aseasonal populations'; Appendix 2). We, therefore, hypothesized that density dependence and opposing episodes of fecundity and viability selection on body resulting from seasonality could cause predictable and repeatable fluctuations in both population and body size. Specifically, we predicted that seasonal fitness trade-offs would cause population size and body size to undergo multigenerational cycles between periods of high abundance, when small individuals predominate, and periods of low abundance, when large individuals are more frequent. Furthermore, we predicted that populations not exposed to viability selection in the nonbreeding season would lack periodic fluctuations in population size and body size.
We tested whether seasonality could result in multigenerational cycles in population size and body size using three experiments ( Figure 1). In the first experiment, we submitted 45 replicate populations to the seasonal treatment described above and tracked the total number of individuals and body size at the end of the breeding and non-breeding season for 58 generations (the 'longterm control' treatment; Figure 1). In the second experiment, we tested the role of viability selection during the non-breeding season by tracking 13 additional populations over 31 generations using a similar protocol to the 'long-term control', but in which we experimentally prevented viability selection in the 'non-breeding' season by providing high levels of food during this season (the 'stop-selection' treatment; Figure 1). This protocol also maintained direct density effects on fecundity and survival, similar to the ones observed in the 'long-term control'. In order to address potential environmental changes in the lab, we conducted a third experiment using the same protocol as in the 'long-term control', but under the same initial conditions and at the same time as the 'stop selection' treatment. This 'short-term control' experiment also had 13 replicate populations tracked over 31 generations ( Figure 1).
In addition to these experiments, we also developed a mathematical model to investigate the contributions of both viability selection and delayed density dependence to population dynamics. The 'stop-selection' experiment was designed to eliminate viability selection, but might have also reduced potential effects of past densities on fecundity and survival. Such delayed density-dependent effects can also cause populations to cycle (Stenseth et al., 2003;Yan et al., 2013), or lead to more complex dynamics, such as chaos (May 1973). One potential mechanism for delayed density dependence are carry-over effects, which we have previously identified in this seasonal system (Betini et al., 2013a). Thus, we first investigated if lag effects were present in all three experiments, and then used the mathematical model to understand whether they played a role in the dynamics of their populations.

Long-term control
Over 58 generations, the 45 replicated seasonal populations showed a predictable increase in abundance during the breeding season, where the food medium allowed females to lay eggs, and a decline in the subsequent non-breeding season (Figure 2a), as is typical of many natural seasonal systems. However, the autocorrelation functions (ACF) also revealed that these short, seasonal cycles were embedded within longer multigenerational cycles where average population size fluctuated 3fold (insert in Figure 2a). In these populations, the ACF function was characterized by stationary periodic dynamics, which resulted in an oscillatory decay to zero (Figure 2a inset).
To investigate the presence of viability selection for small body size and whether this selection was density-dependent, we measured female dry weight in 38 generations from 25 different populations. As expected, there is a negative correlation between population size at the end of the nonbreeding season and body size at the end of the non-breeding season (Pearson's product-moment correlation = À0.64; t = À4.48, p<0.001), suggesting that density negatively impact body size in the non-breeding season. Mean survival during the non-breeding season was 71% (±0.21 SD) and survival was density-dependent (b survival = À0.001, t = À27.73, p<0.001). Mean female dry weight was significantly lower after the non-breeding season (0.279 mg, n = 3620 females) than before the nonbreeding season (0.381 mg; n = 5258 females; standardized values = 0.577 before and À0.566 after the non-breeding season; Welch t-test: t = À35.90, df = 1,589.240.24, p<0.001; Figure 2b) and this viability selection was density-dependent ( Figure 2c; Table 1, Appendix 3). That is, when population size was high at the start of the non-breeding season, there was stronger selection for smaller flies and this was driven by changes in mean dry weight after the non-breeding season rather than changes in the mean dry weight before the non-breeding season (Appendix 3). Average dry weight measured after the non-breeding season also showed multigenerational cycles (Figure 2d), as indicated by the autocorrelation function ( Figure 2d, inset), varying between average peaks of 0.32 mg and lows of 0.23 mg.

'Stop selection' experiment
In order to test the role of viability selection in these multigenerational cycles, we experimentally eliminated viability selection during the non-breeding season in 13 additional populations. Unlike the 'long-term controls' (Figure 2a), there was no evidence of multigenerational population cycles in these 'stop selection' populations (Figure 3a inset). In addition, body size did not significantly decline after the non-breeding season (Figure 3b; average body size was 0.337 mg before and 0.331 after the non-breeding season; n = 1290 and n = 689 females, respectively; standardized values: 0.028 before and À0.052, respectively; Welch t-test: t = À1.733, df = 1,479.400.40, p=0.081). There was also no evidence of density-dependent selection ( Figure 3c; Table 1) and no evidence of cycles in body size (Figure 3d inset).

'Short-term control'
Over 31 generations, and similar to the 'long-term control' (Figure 2a), the 'short-term control' exhibited evidence for multigenerational cycles (Figure 4a). These 'short-term control' populations also experienced viability selection that was density-dependent. Overall, female dry weight significantly decreased from an average of 0.371 mg before the non-breeding season (n = 968 females) to 0.276 mg after the non-breeding season (n = 623 females; standardized values = 0.487 before and À0.761 after the non-breeding season; Welch t-test: t = À30.601, p<0.001; Figure 4b), but the magnitude of this viability selection was stronger when densities were higher at the start of the nonbreeding season (Table 1, Figure 4c). dry weight before (blue bars) and after (black bars) the non-breeding season. Vertical bars indicate the mean female dry weight before and after the non-breeding season; (c) Increased population size in the non-breeding season led to stronger directional selection for smaller flies. (d) Time series of female dry weight measured at the end of the non-breeding season over 38 generations. In (a and (d), the autocorrelograms (insets) showed evidence of cycles in both population size and body size. In (a solid blue line denotes mean population size for each generation from all replicates and dotted lines denote ±1 s.d. In (d, the horizontal line within each box represents the median value, the edges are 25th and 75th percentiles, the whiskers extend to the most extreme data points, and points are potential outliers. DOI: 10.7554/eLife.18770.004

Delayed density dependence
We statistically investigated whether fecundity and survival were influenced by density in past seasons (i.e. delayed density dependence) in all three experiments: 'long-term control', 'short-term control' and 'stop-selection'. We used linear mixed effect models with vial (population) as a random effect and densities going back up to two generations as fixed effects. In the 'long-term control' and 'short-term control', fecundity and survival were influenced by density in the current and past season ( Table 2 and Table 3). In contrast, the 'stop-selection' experiment showed evidence of the delayed effects on fecundity ( Table 2), but not on survival ( Table 3), meaning that the 'stop-selection' treatment eliminated both viability selection on body size as well as delayed effects of density on survival.

Mathematical model
A mathematical model including delayed density dependence as well as the effects of body size on survival (viability selection) and fecundity resulted in multigenerational cycles in population size (before red line in Figure 5a), similar to those observed in our 'long-term control' and 'short-term control' (Figure 2a and Figure 4a, respectively). The model without viability selection on body size and delayed density dependence (i.e. including only effects of current abundance on fecundity and survival), resulted in the elimination of multigenerational cycles (after red line in Figure 5a), as observed in our 'stop selection' populations ( Figure 3a). The exclusion of only viability selection eliminated the fitness trade-off in body size and allowed larger flies to survive to breed. This led to unstable population dynamics (i.e. the population crashed after 10 generations; Figure 5b) because larger flies have greater fecundity and there was a negative interaction between body size and abundance. The exclusion of delayed density dependence alone eliminated the cycles (Figure 5c), suggesting that viability selection and delayed density dependence are necessary for these persistent multigenerational cycles to occur. The model with low heritability (h 2 = 10 À5 ) also generated multigenerational population cycles ( Figure 5d).

Discussion
Our results provide empirical and mathematical evidence that the interplay between density-dependence and evolutionary trade-offs caused by seasonality can have important consequences for population dynamics. In our experimental system, seasonal variation in resources resulted in a fitness trade-off and multigenerational population cycles. These cycles were observed in the absence of predation, long-term fluctuations in resources, or any explicit negative-frequency dependence, all of which are known to cause regular fluctuations in population size (Lindströ m et al., 2001;Yan et al., 2013). Thus, the emergence of population cycles caused by three common characteristics of natural populations (seasonality, a life-history trade-off and density-dependence), suggests that these   Figure 4 continued on next page ecological and evolutionary processes could contribute to fluctuations in population size over a wider range of taxa and environments than has previously been considered. Our experimental elimination of viability selection and delayed density-dependence and our mathematical model confirmed the importance of both of these evolutionary and ecological processes in the persistence of multigenerational population cycles resulting from seasonal environments.
Recently, it has been proposed that eco-evolutionary dynamics are essential for understanding oscillations in population size (Hairston et al., 2005;Hiltunen et al., 2014), but evidence for evolutionary change (i.e. genetically based change in ecologically relevant traits) to feedback and influence demography and population dynamics is rare (Schoener, 2011;Hendry, 2013); but see (Cameron et al., 2013). Our results suggest that the evolutionary response to selection might not be required to observe feedback between ecological and evolutionary processes, as long as selection is strong. In flies, as in many other species, body size is highly heritable (Prout and Barker, 1989), and it is reasonable to suppose that, in our populations, offspring from smaller flies tend to be smaller. However, even when genetic variation is low, strong selection can still affect trait distributions within generations (Grafen, 1988), potentially affecting population size. In our experimental Figure 4 continued (insets), (b) female dry weight before the non-breeding season (blue bars) was higher than after the non-breeding season (red bars; vertical bars indicate the mean female dry weight before and after the non-breeding season). (c) increased population size in the non-breeding season led to stronger directional selection for smaller flies. In (a) solid black line denotes mean population size for each generation from all replicates and dotted lines denote ±1 s.d. In (d, the horizontal line within each box represents the median value, the edges are 25th and 75th percentiles, the whiskers extend to the most extreme data points, and points are potential outliers. DOI: 10.7554/eLife.18770.007 Table 2. Parameter estimates obtained from linear mixed effect models to investigate the effects of current and past density on fecundity in the 'long-term control', 'stop-selection' and 'short-term control' experiments. B refers to population size at the beginning of the breeding season in the current season and NB refers to population size at the beginning of the previous non-breeding season system, selection for body size during the non-breeding season was strong enough to drive observed cycles, potentially in the absence of any genetic change across generations. Simulations with our mathematical model with essentially zero heritability (h 2 = 10 À5 as opposed to 0.3) also generated multigenerational population cycles because viability selection changed the distribution of body size, which subsequently affected fecundity. Thus, changes in the trait distribution within a generation caused by seasonal fitness trade-offs could generate feedback loops between ecological and evolutionary process in traits across generations. We showed that variation in non-breeding abundance can cause density-dependent selection, but our mathematical model suggested that long-lasting effects of past density on survival were essential to generate the multigenerational cycles. These effects were strong in both the 'long-term control' and 'short-term control', but not in the 'stop-selection' treatment. One mechanism that can generate delayed effects is a density-mediated carry-over effect. We have shown how competition during the non-breeding season can negatively impact the physiological conditions of the survivors, reducing their breeding output in the following breeding season (Ratikainen et al., 2008;Betini et al., 2013a). This carry-over effect could also negatively impact survival if individuals during the non-breeding season were in lower body condition because of high abundance at the beginning of the season. Delayed density dependence can also arise from maternal effects, whereby offspring Table 3. Parameter estimates obtained from linear mixed effect models to investigate the effects of current and past density on survival in the 'long-term control', 'stop-selection' and 'short-term control' experiments. NB refers to population size at the beginning of the non-breeding season in the current generation. B1, NB1, B2 and NB2, refers the population size at the beginning of each season going back 1 or two generations, respectively. In the 'long-term control', R 2 LMM ( attributes are affected by parental abundance, which have previously been shown to generate population cycles (Plaistow and Benton, 2009). Both carry-over effects and maternal effects have the potential to interact with viability and fecundity selection because they can affect individual fecundity and survival via changes in body size and condition. Our mathematical model also suggested that the delayed density dependence alone would cause short cycles and population crashes after only a few generations. It is clear from the model that the addition of viability selection allowed for sustained multigenerational cycles, as was observed in our empirical results. These results are in general agreement with previous studies. For example, in a predator-prey system, evolutionary changes in the prey population extended density cycles compared to the typical quarter-period phase lags between prey and predator densities (Yoshida et al., 2003). This evolutionary effect is now believed to be widespread in lab systems (Becks et al., 2012;Hiltunen et al., 2014). Thus, it is possible that evolutionary processes affecting trait distributions within and between generations can buffer strong negative feedback caused by ecological processes that are typical in many natural systems (e.g. overcompensation caused by density dependence), representing a different form of evolutionary rescue (Bell, 2013;Carlson et al., 2014).
Opposing viability and fecundity selection on body size was necessary for the persistence of multigenerational population cycles. Such fitness trade-offs between seasons are likely more common than currently appreciated (Schluter et al., 1991;Schmidt et al., 2005;Conde, 2006, Schmidt andPaaby 2008;Kingsolver and Diamond, 2011;Behrman et al., 2015), especially given that many organisms use vastly different habitats over the annual cycle. Although body size caused fitness trade-offs between seasons in our experimental system, many other traits might also exhibit seasonal trade-offs. For example, sexually selected traits enhance mating success, but might also reduce survival during the non-breeding season (Emlen, 2001). Alternatively, increased development rates can accelerate age of first reproduction, but are often associated with reduced survival during the adult stage (Stearns, 1992;Kingsolver and Huey, 2008).
Seasonal fitness trade-offs also represent a form of balancing selection, which can maintain increased levels of genetic variation at evolutionary equilibrium (Barton and Keightley, 2002). Seasonal fluctuations in selection might allow for the maintenance of increased adaptive potential in populations (Shaw and Shaw, 2014), facilitating more rapid responses to directional environmental change (Bradshaw and Holzapfel, 2006;Bell, 2010;Huang et al., 2016). This temporal variation in selective pressures might also help to explain why some species are more successful at colonizing new areas, as suggested by studies of seasonal changes in genetic variation in wild populations of D. melanogaster (Schmidt and Conde, 2006;Bergland et al., 2014;Behrman et al., 2015). Thus, the consequences of seasonal fluctuations in resources are, therefore, not restricted to purely ecological pathways. Seasonality can also alter the strength and direction of natural selection and might maintain a population's adaptive potential over longer evolutionary timescales.
Examples of delayed density dependence causing cycles have been documented in natural populations subjected to strong seasonality (Merritt et al., 2001;Stenseth et al., 2003;Yan et al., 2013). We showed how opposing episodes of selection driven by a life-history trade-off arise in seasonal environments, which could interact with delayed effects to influence population size. Recently, it has been proposed that results from lab studies on consumer-resource dynamics were likely confounded by fast evolutionary changes (Hiltunen et al., 2014). The same could be true for lab and field studies on population cycles, if changes in trait distributions are linked to population size. To understand whether this is a common phenomenon in natural populations, one needs to tease apart the evolutionary from the ecological consequences of variation in density. As in our lab system, this is a constraint that needs to be overcome with either improved controlled experiments or mathematical approaches. Nevertheless, given that variation in density is common in many natural systems, measurements of selection and population density in both seasons will provide a better understanding of the dynamics of natural populations and how environmental change might alter such dynamics.

Materials and methods 'Long-term control' experiment
To simulate seasonality in populations with non-overlapping generations, we changed food composition to generate two distinct 'seasons' (n = 45 populations). During the breeding season, we placed adults in vials (28 Â 95 mm) with a dead yeast-sugar medium to lay for 24 hr (day 0), after which adults were discarded and larvae were allowed to mature to adults (16 days). Individuals were then transferred (day 17) to a non-breeding environment for 4 days. The non-breeding season consisted of an empty vial of the same size as the breeding vials, but food was provided by a pipette tip filled with 0.200 ml of 5% water-sugar solution per day from the top of the vial. This solution was sufficient for many flies to survive (~95% survival at low population size) but did not allow females to lay eggs (Bownes and Blair, 1986;Betini et al., 2013a). In both seasonal treatments, a consistent amount of food was provided regardless of population size, which caused reproduction and survival to be density-dependent (Betini et al., 2013a). Each of 45 replicate populations was repeatedly transferred between breeding and non-breeding environments for 58 generations and together they are referred to as 'long-term controls'. The quality and amount of the medium during the breeding season mimics the scenario well explored in other Drosophila systems (Mueller and Joshi, 2000), where the food is of lower quality and more limited for adults compared to larvae. Dynamics of populations experiencing only this medium (i.e. only the breeding season) do not show evidence of cycles (Mueller and Joshi, 2000); (Appendix 2-figure 1).

'Stop-selection' experiment
To investigate whether opposing selective pressures caused multigenerational cycles in seasonal populations, we experimentally stopped viability selection for a smaller body size in the non-breeding season while preserving density-dependent survival. To do this, we exposed 13 new replicate populations to the same seasonal change in food resources (described above) over 31 generations, but provided unlimited access to food during the non-breeding season (0.8 ml/day instead of 0.2 ml/day; initial population size was 5 males and 5 females). Average mortality was reduced from 28% (±0.4; mean ± s.e.) in the 'long-term control' populations to 0.5% (±1) in these new 'stop selection' treatments. After four days in the non-breeding season, we haphazardly selected the survivors that moved to the breeding season. To preserve the density-dependent survival process during the nonbreeding period, the number of survivors for each population was calculated based on a logistic non-breeding survival function parameterized from our 'long-term control' populations (Wilson, 1994; Figure 6) where Su is survival (number of survivors at the end of the non-breeding season divided by the total number of individuals at the beginning of the non-breeding season), N is the population size at the beginning of the non-breeding season, and v and w are constant to be estimated from the data. The function was parameterized with our 'long-term control' populations (45 replicates over 43 generations) using the non-linear function nls in R (R Core Team, 2015). Individuals that moved to the following breeding season were haphazardly selected.

'Short-term control' experiment
Our 'long-term control' and 'stop-selection' populations were initiated at different periods and with different population sizes. Thus, differences in the lab environment and initial conditions could have influenced the dynamics of these populations. For these reasons, we also initiated an additional 13 replicate populations at the same time and same initial population size as the 'stop-selection' populations. In these populations, we used the same protocol as in the 'long-term control', but initial population size (5 males and 5 females) and number of generations (n = 31) was the same as in the 'stop-selection' experiment.

Measuring viability selection on body size
We measured linear selection differentials (S) for female body weight (i.e. population mean dry weight after -population mean dry weight before the non-breeding season) in all three experiments ('long-term control', 'stop-selection' and 'short-term control'). We used female dry body weight as a proxy for body size. We experimentally confirmed that the mean body size declined along with mean dry weight in surviving females after the non-breeding season by placing individuals from the stock population in the non-breeding season in either high (n = 300) or low (n = 20) abundance. We measured the thorax size in 10 individual females sampled before the non-breeding season and 10 females sampled from the population after the non-breeding season. Mean thorax size after the non-breeding season was significantly lower than mean thorax size prior to the non-breeding season (Welch t-test: t = À2.87, df = 37.66, p=0.006), whereas there was no significant differences in average thorax in the low abundance treatment (Welch t-test: t = À0.55, df = 36.96, p=0.586; average size before = 1.318 mm; average size after the non-breeding season for the high abundance treatment was 1.225 and 1.297 mm for the low abundance treatment).
We obtained female dry weight from 5% of total population size before and after the non-breeding season. We then sexed, dried and weighed individual females to the nearest 0.001 mg using a microbalance. To calculate the linear selection differentials (S) for female body weight, we standardized dry weight measures to a mean of zero and unit variance prior to calculating selection (Brodie et al., 1995). The significance of our estimates of viability selection during the non-breeding season was assessed using a Welch t-test that tested for differences between mean female dry body weight before and after the non-breeding season. We then used a linear mixed effect model (LMM) to investigate whether S was density-dependent. In the LMM, the selection differential calculated for each population in each generation was the response variable, number of individuals at the beginning of the non-breeding season was fitted as a fixed effect and population was fitted as a random effect.
In the 'long-term control', we measured female dry weight in 38 generations from an arbitrary number of replicate populations (16 to 25 population per generation) in generations 9, 10,15,16,21,22,25,26,. Total number of individuals measured was 5258 before and 3620 after the non-breeding season. In the 'stop-selection' and 'short-term control', we sampled females in seven populations in generations 1 to 25 (n = 1591 females and n = 1979 females in the 'short-term control' and 'stop selection' treatments, respectively), with the exception of generation 13 in the 'short term control'.

Documenting multigenerational cycles
To investigate whether population size cycled in our experimental populations, we used the autocorrelation function method (Turchin and Taylor, 1992;Box et al., 2013). Cycles can be inferred by investigating the shape of the estimated ACF: population cycles are characterized by stationary periodic dynamics, which result in an oscillatory decay to zero of the ACF estimates. We calculated the ACF for each treatment using the mean population size across all replicates at the start of the breeding and non-breeding season in each generation. For all-time series, we excluded the first three generations to avoid any transient effects caused by the initial population size. We also detrended the time series to eliminate any overall temporal trends in abundance by subtracting the fitted values from a linear regression of average population size on generation ('long-term control': b = 0.468, s. e.=0.177, t = À2.65, p=0.011; 'short-term control': b = À1.855, s.e.=0.686, t = À2.71, p=0.01).
To investigate temporal changes in body size in the 'long-term control' and 'stop-selection' populations (measured as female dry weight), we used the same autocorrelation function method described above for the population size, including detrending the data.

Documenting delayed density dependence
We investigated the effects of past abundance on fecundity and survival with LMM with vial (population) as a random effect in all three experiments: 'long-term control' and 'short-term control' and 'stop-selection' treatment. In previous work, we have shown that fecundity in this system is influenced by an interaction between current abundance at the beginning of the breeding season and at the beginning of the previous non-breeding season (Betini et al., 2013a(Betini et al., , 2013b. Thus, in this study, we also investigated if this interaction significantly explained variation in fecundity in all three experiments. Fecundity was calculated as the number of individuals at the end of the breeding season (i.e. number of offspring) divided by the number of individuals at the beginning of the same breeding season (number of parents). Because we did not know how environmental effects (e.g. physiological condition) could influence survival in our system, we incorporated the breeding and non-breeding abundance going back two generations as explanatory variables in the model to explain survival. Lag effects have been documented in natural populations subjected to strong seasonality (Merritt et al., 2001(Merritt et al., , 2001Stenseth et al., 2003). Survival was calculated as the number of individuals at the end of the non-breeding season divided by the number of individuals at the beginning of the non-breeding season. Explanatory variables were standardized before analysis by subtracting the sample mean from each observation and dividing each value by the sample standard deviation and response variables were log transformed.

Mathematical model
We developed a mathematical model to investigate the contributions of both viability selection and delayed density dependence to population dynamics. The model linked quantitative trait evolution with demography and was similar to integral projection modelling approaches that are commonly used to study ecological and evolutionary processes concurrently (Slatkin, 1979(Slatkin, , 1980Ellner and Rees, 2006). In particular, the number of individuals with body size y at the next time step t + 1, n (y, t + 1), was a function of a projection function k (y) and n (y, t), such that n (y, t + 1) = k(y) * n (y, t). The projection function k (y) accounted for the additive genetic and environmental contributions to body size, as well as both density-independent and density-dependent differences in fecundity or survivourship among body sizes, depending on season.

Model description
Two population sizes were modeled in discrete time: population size at the beginning of the breeding season, and population size at the beginning of the non-breeding season within each generation. For each season, the mean breeding value for body size and the mean phenotypic value for body size were modelled, where the phenotypic value of an individual was its breeding value plus a random environmental effect. We assumed that breeding and phenotypic values were normally distributed and that the variance in breeding values and the variance in environmental effects remained constant across generations. The assumption of a constant variance in breeding values assumed that mutational variance and changes in genetic architecture (including linkage and epistatic interactions) could restore variance that was initially depleted by selection. Slatkin's model provides for changes in the variance of breeding values due to selection and considers the effect of genetic architecture, such as linkage (Slatkin, 1979(Slatkin, , 1980. In our experiment, we do not have information on the genetic architecture for body size, nor mutational variance. We did not see a decline in the magnitudes of the rates of change in body size in the experiment, which is consistent with the maintenance of heritability. The fecundity of an individual during the breeding season was a function of its body size and current and past effects of density. Past density influenced fecundity via changes in the physiological conditions of individuals that spent the previous non-breeding season at high densities but survived to breed at low densities (i.e. a carry-over effect that was modelled as the interaction between density at the end of the breeding season in the previous generation and current breeding density (Betini et al., 2013a(Betini et al., , 2014. Current density could also negatively influence fecundity via densitydependent effects. There was a stronger negative effect of density on fecundity for large versus small individuals, although large individuals have higher intrinsic fecundity ( Figure 3 in the main text). During the 'stop selection' phase of the simulations, we decreased carry-over effects (parameters ' 5 and ' 6 ) by 25% given that flies had unlimited amount of food in the 'stop selection' treatments and were likely to be in better physiological condition than those in the 'short-term control' treatments. Parameters used in the fecundity function were chosen such that fecundities as a function of body size and densities were within the normal range observed in D. melanogaster.
The survival of an individual during the non-breeding season was also a function of its body size (when viability selection occurs) and current and past effects of density. As indicated by our previous experiment, larger flies had lower survival and this effect was magnified by increase density (Betini et al., 2014). We included the effects of past density beyond the ones described above, going back two generations. During the 'stop selection' phase of the simulations, survival was not a function of body size nor past density, but was instead only affected by current density, following the logistic equation Su for survivorship and the experimental design. To accomplish this and to ensure proper scaling of survival, we assigned all individuals the same phenotype in the survivorship function. Parameters for the survival function were based on values obtained in the long-term control ( Figure 2a in the main text).
In the context of an individual's breeding and phenotypic values, an individual with breeding value x had an environmental effect y added with a mean effect that was a negative linear function of density and was normally distributed with variance V[e]. Consequently, each breeding value expressed a distribution of phenotypes as a function of the current environment.
In the context of survivorship, since each breeding value expressed a distribution of phenotypes, each breeding value expressed a distribution of survivorship and the model integrated across environmental effects to get the mean survivorship associated with a breeding value. This integration of survivorship across environmental effects gave the value of the projection function for a given breeding value. The projection function in the context of fecundity was modelled similarly; breeding values were transformed into phenotypes via a distribution of environmental effects, which when integrated gave the mean fecundity of individuals with a particular breeding value for body size.
Population size across seasons was the product of the current population size times the average fecundity (breeding season) or average survivorship (non breeding season) at the phenotypic level. Below is a list of variables, parameters and functions used in the model.

Variables
X i -population size at the beginning of the breeding season i generations ago Y i -population size at the beginning of the non-breeding season i generations ago b X -mean breeding value for body size at the beginning of the breeding season b Y -mean breeding value for body size at the beginning of the non-breeding season p X -mean phenotype for body size at the beginning of the breeding season p Y -mean phenotype for body size at the beginning of the non-breeding season z -body size Note, a generation consists of a breeding season and then a non-breeding season. At the beginning of the breeding season the mean breeding value b X ð Þis a function of densities in previous generations, such that i > 0.
Parameters V A ¼ 0:003 -additive genetic variance for body size V E ¼ 0:007 -environmental variance for body size such that h 2 ¼ 0:30 (Prout and Barker, 1989) Fecundity function: Survival function: Integral projection model: ÞNðy; e X ðX 1 Þ; V E Þdydx In the equations above N z; m; v ð Þ is the probability density of a normally distributed random variable with a mean of m and variance v. Functions e X N ð Þ ¼ Àl X Nand e Y N ð Þ ¼ Àl Y Ngive the average environmental effect on the phenotype.

Other parameters and their values
' 1 ¼ 1:5 (intrinsic fecundity for small flies) ' 2 ¼ 2:8 (rate of increase in fecundity with body size) ' 3 ¼ 0:00007 (baseline rate of decline in fecundity with density for small flies) ' 4 ¼ 0:04588 (rate of increase in the magnitude of the decline in fecundity as body size increases) ' 5 ¼ 1 (constant characterizing the effects of population size at the end of the breeding season one generation ago; i.e. carry-over effects) ' 6 ¼ 0:001 (constant characterizing the interaction between current population size and population size at the end of the breeding season one generation ago) 1 ¼ 0:005 (rate of decline in survivorship due to population sizes one generation ago) 2 ¼ 0:0026 (rate of decline in survivorship due to population sizes two generations ago) 3 ¼ 0:7 (constant governing negative effect of the interaction between the current population size and body size) 4 ¼ 0:7 (constant governing the shape of the negative effect of the interaction between current population size and body size on survivorship) 5 ¼ 6 (a second constant governing the shape of the negative effect of the interaction between current population size and body size on survivorship) 6 ¼ 350 (a third constant governing the shape of the negative effect of the interaction between current population size and body size on survivorship) l X ¼ 0:00001 (rate of decline in body size [with density], i.e. the environmental effect of density on body size due to density at the beginning of the breeding season) l Y ¼ 0:00025 (rate of decline in body size [with density], i.e. the environmental effect of density on body size due to density at the beginning of the non-breeding season) All LMM were performed using the lmer function from the lme4 package (Bates, 2010) and p-values were obtained using the lmerTest package (Kuznetsova et al., 2014). All analyses were conducted in R (R Core Team, 2015). Marginal (R 2 LMM(m) ) and conditional (R 2 L<MM(c) ) variance for the LMMs were calculated with MuMIn package according to (Nakagawa and Schielzeth, 2013). R 2 LMM (m) is the variance on the response variable that is explained only by the fixed effects and R 2 LMM(c) is the variance that is explained by both fixed and random effects (Nakagawa and Schielzeth, 2013).