## Abstract

Understanding the variation in selection pressure on key life-history traits is crucial in our rapidly changing world. Density is rarely considered as a selective agent. To study its importance, we partition phenotypic selection in fluctuating environments into components representing the population growth rate at low densities and the strength of density dependence, using a new stochastic modelling framework. We analysed the number of eggs laid per season in a small song-bird, the great tit, and found balancing selection favouring large clutch sizes at small population densities and smaller clutches in years with large populations. A significant interaction between clutch size and population size in the regression for the Malthusian fitness reveals that those females producing large clutch sizes at small population sizes also are those that show the strongest reduction in fitness when population size is increased. This provides empirical support for ongoing *r*- and *K*-selection in this population, favouring phenotypes with large growth rates *r* at small population sizes and phenotypes with high competitive skills when populations are close to the carrying capacity *K*. This selection causes long-term fluctuations around a stable mean clutch size caused by variation in population size, implying that *r*- and *K*-selection is an important mechanism influencing phenotypic evolution in fluctuating environments. This provides a general link between ecological dynamics and evolutionary processes, operating through a joint influence of density dependence and environmental stochasticity on fluctuations in population size.

## 1. Introduction

One of the central issues in ecology has been the influence of density dependence upon fluctuations in population size, which lead to controversies lasting for several decades [1]. An important contribution was provided by May and co-workers who showed that even simple population models could produce extremely complicated population dynamics [2] and that most populations showed patterns in population fluctuations that clearly indicated regulatory effects of density dependence [3]. However, even though ecologists have for a long time realized the importance of density dependence, estimating the form of density regulation and strength of density dependence in natural populations have proved difficult [4,5]. This occurs because a large sample covering the full range of variation in population size is necessary for obtaining reliable estimates of how a change in population size affects the population size the following year(s) [6]. This requires long time series, preferably based on precise population estimates [7].

In population genetics, classical models assume a constant population size. In an attempt to include more ecological realism in population genetic models MacArthur [8] showed that in a stable environment under density-dependent selection, evolution will maximize equilibrium population size or the carrying capacity of the population. This important result formed the foundation for the introduction by MacArthur & Wilson [9] of *r*- and *K*-selection as an important selective agent in natural populations. Based on the assumption of a trade-off between traits affecting selection at small population sizes and those favoured at larger population sizes, they argued that high population growth rates should be favoured when the number of individuals was small (*r*-selection). By contrast, at large population sizes selection should favour highly competitive phenotypes, resulting in an increased carrying capacity (*K*-selection).

This introduces a complication different from that in ecological analyses of density dependence; density-dependent selection requires that changes in population size affect the fitness of different phenotypes differently. For example, the occurrence of *r*- and *K*-selection in a population requires that the fitness of phenotypes favoured at small population sizes are more strongly influenced by a change in population size than the fitness of those phenotypes with high mean individual fitness close to the carrying capacity *K.* Thus, analyses of density-dependent selection must include time series of how fitnesses of different phenotypes are affected by fluctuations in population size, which requires estimates of second-order interaction effects between population size and phenotypic variation on fitness that may be small compared with the first-order effects. In addition, as in ecology, the effects of density dependence should be separated from the influence of environmental variation on the distributions of phenotypes. Consequently, this makes data requirements for statistical analyses of density-dependent phenotypic selection in natural populations extensive because long time series of individual phenotype-specific differences in fitness are needed.

A central focus in evolutionary biology has been to understand factors affecting the rate of phenotypic evolution [10]. A large number of studies covering a wide variety of taxa document that phenotypic selection acting on quantitative characters occurs regularly in natural populations [11–13], yet processes explaining variation in the strength of this selection in time and space are poorly known [11,13–15]. Because populations of even closely related species may show large differences in the magnitude of population variability [16,17], density-dependent selection may provide a general mechanism with important implications for how fast phenotypes will change over time. For instance, if there is no density dependence, the long-run growth rate of the population is maximized [18] and stasis will only occur when there is balancing selection around an intermediate optimal phenotype (see Chevin *et al*. [19] and references therein for analyses of the effects of stabilizing selection on phenotypic evolution in a fluctuating environment). Basically, this corresponds in the terminology of MacArthur & Wilson [9] to a form of *r*-selection. By contrast, selection may also favour individuals with phenotypes that are little affected by population density and thus produces *K*-selection. Although these cases represent *r*- and *K*-selection, this does not necessarily imply that there are fluctuations within a population between these two types of selection. However, if there is genetic variation in growth rate as well as strength of density regulation with a trade-off between phenotypes favoured at small population sizes and those that are less affected by population density, fluctuating *r*- and *K*-selection around an intermediate phenotype, dependent upon the environmental stochasticity may occur [20–23]. Thus, the influence of fluctuations in population size on phenotypic selection may be crucial for understanding basic features of phenotypic evolution [24].

Our aim with the present study is to use a unique long-term demographic dataset of great tits in the Netherlands to analyse how phenotypic selection depends upon population size. We choose the total number of eggs produced per female during the breeding season as the phenotypic trait studied. Previous studies have shown stabilizing selection around an intermediate clutch size in great tits [25,26], which is affected by variation in population size [27]. We use the model of Engen *et al*. [21] to analyse the joint process for mean character and log population size in a fluctuating environment. A general result obtained from this model is that with reasonable trade-offs preventing unlimited growth of fitness, evolution approaches stationary fluctuations of the mean character in the population, with the expected response always in the direction that increases average population size. Hence, the mean phenotype during a long time interval represents a compromise between phenotypes with large growth rates at small population sizes and those with weak density regulation [21,28]. In addition, stochastic effects make the mean fluctuate continuously, generating fluctuations in population size. This provides a tool for empirical analyses of density-dependent selection if the strength of density dependence *γ*(*z*) is not constant among phenotypes. For *r*- and *K*-selection to occur, there must exist a trade-off between the growth rate at low densities *r* and *γ* (or between *r* and *K*) because otherwise the evolution of the mean phenotype would always increase the mean growth rate and decrease the mean strength of density regulation resulting in correlated selection for large and *K*. Thus, the requirements of such phenotypic trade-offs between major life-history characters means that empirical analyses of evolution by means of *r*- and *K*-selection are challenging because obtaining reliable estimates of density-dependent selection as well as genetic covariances among phenotypic characters are in general difficult in natural populations [29]. Here we provide evidence that fluctuations in population size can act as an important selective agent of mean clutch size in great tits. We show, using a novel approach based on a parametrizing of a model for *r*- and *K*- selection including realistic ecological dynamics [21], that there is phenotypic variation in the strength of density regulation as well as in the population growth rate at small densities. More importantly, we show that there is an interaction effect between population size *N* and mean phenotype on the Malthusian fitness. Consequently, when the population size fluctuates, this generates *r-* and *K-*selection as predicted by MacArthur & Wilson [9], resulting in evolution towards an intermediate clutch size that varies little through time. The evidence for such density-dependent selection is important because it provides empirical support for a general mechanism for a reciprocal coupling between ecological and evolutionary processes [30].

## 2. Material and methods

### (a) Study population and field procedures

We analyse selection on the number of eggs produced per season based on data from a long-term (1955–present) population study of great tits at the Hoge Veluwe National Park, the Netherlands. The area consists of mixed pine-deciduous woodland on poor sandy soils. To ensure that availability of artificial nest sites did not limit population size, we provided a surplus of nest-boxes [31]. The adult survival rate of great tits in this population was around 0.40 [32], strongly dependent on the beech crop cycle. The population dynamics of great tits is dependent on a combination of density-dependent effects as well as environmental fluctuations especially caused by variation in beech crops and winter temperatures [33–35]. In our study populations, immigrants constitute an important part of the population and strongly affect population dynamics [31].

We checked nest-boxes for nest building, number of eggs or number of nestlings at least once a week from April onwards. The parents were caught when the nestlings were 7–10 days old. When already ringed they were identified, if unringed they were given a metal ring with a unique number. Nestlings were ringed on day 7. Female great tits can produce a second brood within the same season (i.e. raising a new brood after successfully fledging the first brood), although the frequency of this double-brooding has declined in recent decades in this population [36]. Unringed females were not included in the analyses, as their survival to future breeding seasons could not be determined. The probability of recapturing a female was very high (average = 98.7%) [37]. Manipulated broods (both parents and offspring) were excluded from all analyses.

### (b) Model

Our analysis is based on females' individual fitness *W* = *I* + *B*/2, which is an individual's contribution to the next year's population. Here *I* is an indicator of female survival with value 1 for survival and otherwise 0, while *B* is its number of recruits of both sexes entering the breeding population the following year or later. Individual fitness is a stochastic variable [38] with some distribution expressing stochastic variation among individuals within years, representing demographic stochasticity, depending on phenotype *z*, population size *N* as well as the environment defined by some vector **ε**. Writing the expected fitness of an individual of phenotype *z* at population size *N* in environment *ɛ* as the function *M*(*z*, *N*, **ε**) is then the Malthusian parameter of a hypothetical population of individuals with phenotype *z* and population size *N* under constant environmental conditions *ɛ*. Averaging *W* over demographic as well as environmental stochasticity using the notation **E**_{ɛ} for the expectation, we define where *r*(*z*, *N*) is the deterministic growth rate of a hypothetical population of individuals with phenotype *z* at density *N*. Averaging the conditional Malthusian parameter over the fluctuating environment *ɛ* we obtain the mean Malthusian parameter for given *z* and *N*, which is the long-run growth rate of a hypothetical population with parameters *z* and *N*. The first-order approximation to the mean Malthusian parameter is
where the last term, which is generated by the common environmental noise, will depend slightly on *z* and *N*, but is here approximated by its mean value through time.

An important general result [21] is that expected evolution of the mean phenotype subject to density-dependent selection always increases the function
2.2Here is the long-run growth rate of the population in the absence of density regulation so that where is the deterministic growth rate of the mean phenotype in the average environment, is the environmental variance, is the average strength of density dependence in the population growth rate and *g*(*N*) is an increasing function of *N* describing the density regulation (for theta-logistic density regulation *g*(*N*) = *N ^{θ}* and for the logistic model

*g*(

*N*) =

*N*). The parameter is then the carrying capacity in the corresponding deterministic model so that The expected evolution of the mean phenotype under density-dependent selection will be towards the value

*z*

_{opt}maximizing which is the mean of

*g*(

*N*) for a given [21]. The expected evolution of thus corresponds to an increase of but environmental stochasticity always perturbs it back to values slightly below the maximum so that has temporarily stationary fluctuations around the optimal phenotypic value

*z*

_{opt}.

We adopt the model of Engen *et al*. [21] writing defining a logistic dynamic model [17] for a hypothetical population of individuals with phenotype *z*. Here *r*_{0}(*z*) is the deterministic growth rate at small population sizes, while *γ*(*z*) defines the strength of density regulation. Because large variation in the deterministic growth rates at small densities with variation in clutch size is known to occur [25,27], it is required to use a second-degree approximation also allowing for variation in the slopes. The deterministic growth rate then has a maximum if *β*_{3} < 0. There is, on the other hand, much less knowledge about the phenotypic variation in the strength of density regulation. We, therefore, believe that the best choice for achieving our main goal of detecting a decrease in *γ*(*z*) with increasing phenotype, which is required for density-dependent selection to occur, is to avoid over-parametrization by using a simple linear relationship to describe how the strength of density dependence depends on phenotype *z*.

### (c) Estimation procedures

We perform the statistical analysis on the individual fitness of females *W* using a mixed log-linear model with 2*W* = 2*I* + *B* as dependent variables taking integer values 0, 1, 2,…. Then where and the temporal variance of *ε* is the environmental variance Writing **W** for the vector of observed individual fitness in a given year with population size *N* and where superscript T denotes matrix transposition, the model takes the form , where **E** denotes the expectation. The left side is here the column vector with components **E**ln(2*W*) defined for each observed female with observed dependent variable 2*W* a given year, the rows of **X** are the individual vectors of these individuals and *U* is a normal environmental noise variable with zero mean and variance assumed to be independent among years.

### (d) Ecological and evolutionary dynamics

Engen *et al*. [21] showed that the dynamics of total population size *N* in this model is of the logistic type determined by the mean value of the parameters across individuals. Let *p*(*z*) be the normal phenotypic distribution of *z* with evolving mean and variance *P* assumed to be constant. The mean parameters are then and Assuming that selection is weak so that *N* changes much faster than the infinitesimal mean and variance of the diffusion approximation for ln *N* for a given is the mean Malthusian fitness and respectively, where is the long-run growth rate for densities close to zero. Then the selection differential is according to Lande & Arnold [39] and Lande *et al*. [20] where *P* is the phenotypic variance. With appropriately defined trade-offs preventing unlimited growth of Engen *et al*. [21] showed that will have stationary fluctuations around the mean phenotype *z**, maximizing the function , where appearing in in general also may depend on the mean phenotype (see Engen *et al*. [21] for details). Writing for the carrying capacity in the corresponding deterministic model, the function that is maximized (equation (2.2)) can alternatively be written as which has the same form as the function found to be maximized in the haploid model of Lande *et al*. [20].

So far, we have ignored migration, which constitutes an important component of the dynamics of this population [31]. This includes treating emigrants as dead individuals, which makes on average negative. To obtain a stationary model for using the results of Engen *et al*. [21], we must assume an immigration rate *μ* replacing by in the expression for The average growth rate at the log scale can be estimated as giving the estimate of the immigration rate so that on average equals zero. We assume that immigration is independent of phenotype and that the mean Malthusian parameter (equation (2.1)) does not differ between immigrants and locally born females. Accordingly, the total clutch size did not differ significantly between locally born females and immigrants (generalized linear mixed model, *P* = 0.223 including year as a random effect).

In the present parametric model, we find
2.3*a*and
2.3*b*which determines the function and the centre phenotype *z**, maximizing . The individual fitness will show a similar pattern of variation with *z* as the mean Malthusian fitness, just differing with a constant *β*_{3}*P*. The estimate of the variance *P* was

### (e) Estimation of heritability

The heritability was estimated following procedures in Husby *et al*. [40] using a residual maximum-likelihood mixed model (animal model). Animal models are based on pairwise relatedness between all individuals and partition phenotypic variance into its additive genetic and environmental variance components [41].

### (f) Simulations of phenotypic evolution

The ecological and evolutionary process can be simulated taking selection and random genetic drift as well as the estimated heritability and stochastic environmental variation into account. The dynamics of *N* are simulated by the mixed log-linear model for individual fitness that was fitted to the data (equation (2.2)), using the actual parameter estimates. The response to selection is whereas the change in breeding value due to random genetic drift is where is the effective population size [42] and *V* is a random variable with expectation 0 and unit variance [38,43]. Following procedures in Sæther *et al*. [44], the estimate of the demographic variance was

The simulations of the temporal variation in breeding values assume that environmental stochasticity affects all phenotypes similarly, that there is no phenotypic plasticity and that the immigrants have the same as natives.

## 3. Results

Our analyses reveal that mean Malthusian fitness (equation (2.1)) as well as the selection differential (the gradient of with regard to ) changed with phenotypic mean in a way that depends on population size (figure 1 and table 1). The mean fitness decreased considerably with population size (figure 1*a*), confirming strong density regulation. The maximum mean fitness fluctuated considerably with *N* (figure 1*a*), producing selection for larger and smaller clutch sizes at small and large population densities, respectively (figure 1*b*). As revealed by the gradient of the strength of density regulation (−*β*_{5} in equation (2.3*b*)), females producing a larger number of eggs at small population sizes were those more strongly affected by variation in population density. The significant negative sign of *β*_{5} indicates that when the mean clutch size was small, selection will most of the time be towards large clutch sizes but in the opposite direction at large clutch sizes. When clutches were on average of intermediate sizes the direction of selection may fluctuate considerably.

The optimal clutch size for a given population size, i.e. the number of eggs produced that maximizes the fitness of a female, decreased with increasing population size (figure 1*a*, see also Both *et al*. [27]). This relationship was in particular owing to low reproductive success of second clutches in our population [27,36]. Accordingly, none of the *β _{i}* for the first clutches were significantly different from zero, whereas significant stabilizing selection was found for the second clutches (table 1).

Our estimates (table 1) imply that the strength of density regulation increases significantly with increasing clutch size. Accordingly, at small population sizes the Malthusian fitness was highest among females with large clutch sizes, but the fitness of these females decreased rapidly with increasing population size and was very low at large population sizes (figure 1*c*). By contrast, the Malthusian fitness of females laying few eggs showed less variation with changes in population size.

This pattern of density-dependent selection demonstrated selection for high reproductive output at smaller population sizes, but the fitness of these high-fecundity females was at the same time more susceptible to fluctuations in population size as revealed by the significant increase in the strength of density dependence with increasing total clutch size (table 1). This trade-off between *r*- and *K*-selection resulted in stabilizing, stochastic phenotypic selection towards an intermediate clutch size *z** = 15.97 in this population (figure 2), maximizing the function (equation (2.1)) where we have inserted the average immigration rate ensuring that that fluctuates around zero. The maximum value of at *z** is of an order of magnitude similar to the estimated carrying capacity (using a logistic model) in this population. Furthermore, the mean value of *z* during the study period was 11.20, which is close to *z** considering the large stochasticity in the process as well as the uncertainty in the function mainly owing to large standard errors in *β*_{3} and *β*_{5} (table 1). We also note that the mean Malthusian fitness at *z** is still quite small (figure 1*a*). This illustrates that a large immigration rate *μ* = 0.34 is necessary for long-term persistence of this population [31].

In this population, the heritability *h*^{2} of the total number of eggs produced per season is significantly larger than zero (). As a consequence, density-dependent selection results in the evolution of the mean breeding value that will fluctuate slowly around *z** (figure 3). Thus, *r*- and *K*- selection can produce stabilizing selection around an intermediate phenotype even though there is directional selection acting in a single year, environmental stochasticity is substantial (table 1) and demographic stochasticty induces genetic drift that causes random variation in [38].

## 4. Discussion

Our analyses reveal that in great tits the total number of eggs produced per female during a breeding season is subject to fluctuating *r*- and *K*-selection. This occurs because the selection is density dependent (figure 1). When mean clutch size is small, selection will most of the time be towards large clutch sizes and in the opposite direction for large clutch sizes. The significant negative interaction between total clutch size and population size (table 1) indicates that the fitness of high-fecundity females is most strongly affected by large population sizes, causing the breeding value of total clutch size to fluctuate around an intermediate value (figure 3). Thus, such density-dependent selection provides a reciprocal link between current ecological dynamics and evolutionary processes operating at longer time scales [30].

The Malthusian fitness, which is the growth rate on the logarithmic scale of a hypothetical population of individuals with a given phenotype and population size (equation (2.1)) showed a maximum at intermediate total clutch sizes that fluctuates with the population size (figure 1*a*). This is supported by evidence from large-scale brood size manipulation experiments in years which differed in density. They suggest that there exists an intermediate clutch size that maximizes the offspring production of a breeding female [25,27] and that this optimal clutch size varies with density [27]. Furthermore, analyses of the dynamics of populations of great tits have shown that changes in population size between years are determined by a balance between the number of new recruiting offspring produced, which is affected by density-independent factors such as weather and beech mast [32], and density-dependent losses of individuals during the non-breeding season [45]. In particular, variation in juvenile survival strongly affects the fluctuations in population size [33] as well as the fitness of individual females (i.e. the number of offspring produced that recruit to the next generation) [46]. These effects operate through a joint effect of environmental variation on the demography of both local residents and immigrants [31,35].

These analyses are based on several simplifying assumptions. First, we assume that the effects of fluctuations in the environment are independent of the phenotype so that there are no stochastic components in the selection process except what is governed by fluctuations in *N* (see Lande [47], Engen & Sæther [38] and Chevin *et al*. [19] for other stochastic components). Second, there is no phenotypic plasticity with regard to population size in total number of eggs produced [27], which may affect the estimated response to selection (figure 3). Third, we disregard age-dependency which is known to cause transient fluctuations in selection differentials even in short-lived passerines such as the great tit [48]. Fourth, we consider only one phenotypic characteristic, ignoring the influence on total clutch size of correlated selection on other characters. For instance, in great tits variation among females in reproductive success is also affected by the timing of breeding [37], which influences the number of broods that can be produced per season [36]. Accordingly, the density-dependent selection on total clutch size occurs mainly as an effect of population size on the breeding success via second clutches (table 1), which has been shown to strongly affect the reproductive rate of this population [36]. We also assume that the immigrants have the same as the native population. Still, in spite of these assumptions, our model captures essential characteristics of the dynamics of this population (figure 2), enabling us to demonstrate that those females producing a large number of clutches at small population sizes are also those that experience the stronger density regulation as revealed by the significant interaction between the total number of eggs produced and population size (table 1).

As is the case for many reproductive traits of tits [49], there was significant heritable variation in the total number of eggs produced per breeding season, which is a precondition for evolutionary responses to selection. In this case, density-dependent selection to maximize mean population size (equation (2.2)) will result in a stationary distribution of the breeding values of *z*. Simulations of sample paths from this stationary distribution of breeding values show that the temporal variation in mean phenotype will be small (figure 3) because of the trade-off involving higher susceptibility of individuals with high growth rates at small densities to changes in population size (table 1). However, the recorded values of was somewhat smaller than the predicted intermediate value *z** (figure 3), which is expected from the large uncertainties in parameter estimates (table 1) as well as by our ignoring of the effects of phenotypic plasticity and fluctuations in age-distributions. Accordingly, the simulations (figure 3) can only be used to characterize the general properties of the stationary fluctuations in the mean breeding values.

A general problem in evolutionary biology has been to explain the long periods typically seen in the fossil record of small changes in morphology [10,50,51]. In general, this is explained by stabilizing selection around an optimal phenotype [10,52], which assumes stable fitness peaks, corresponding to niches stable over time [50]. Our results indicate that ongoing *r*- and *K*-selection may provide a general mechanism for such stabilizing selection in environments that fluctuate over time (figure 3) because it produces stabilizing selection around an intermediate phenotype *z** (figures 1*a* and 2). This implies that characters may show variation over longer time periods, still with no temporal trend in spite of directional selection on heritable characters [14,50,53].

Evolutionary biology is currently experiencing a period of great synthesis where several important contributions have been made in our understanding of the evolutionary process [10,54,55]. Our results suggest that the inclusion of specific assumptions about the ecological dynamics involved can produce precise predictions about expected patterns of phenotypic evolution. We have recently experienced several advances in our understanding of eco-evolutionary dynamics typically involving feedback processes in which the organism modifies the effects of its environment through *s*trong trophic interactions [56–58]. Here we have shown that feedback loops are not a prerequisite for generating eco-evolutionary dynamics, only density-dependent selection is necessary. Such selection will then alter the evolutionary response to changes in the environment [22,30].

## Ethics

This study was carried out with the approval of the Animal Experimentation Committee of the Royal Dutch Academy of Sciences under licence NIOO 10-07.

## Data accessibility

The data on which the analysis in this paper are based are available from the data repository of the Netherlands Institute of Ecology. This is a data repository which is managed by the Institute and is available via http://data.nioo.knaw.nl/index.php.

## Authors' contributions

S.E. and B.-E.S. developed the model; S.E., B.-E.S. and V.G. performed the statistical analyses; M.E.V. provided the data and contributed to the biological interpretation of the results and all authors contributed to the writing of the paper.

## Competing interests

We declare we have no competing interests.

## Funding

B.-E.S., S.E. and V.G. were supported by the European Research Council (ERC-2010-AdG 268562) and the Research Council of Norway (SFF-III 223257/F50) and MEV by the European Research Council (ERC-2013-AdG 339092).

## Acknowledgements

We are grateful to A. Husby for providing heritability estimates. We thank T. Reed for comments to a previous version of the paper.

- Received October 8, 2015.
- Accepted April 4, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.