Phenotypic evolution in stochastic environments: The contribution of frequency‐ and density‐dependent selection

Understanding how environmental variation affects phenotypic evolution requires models based on ecologically realistic assumptions that include variation in population size and specific mechanisms by which environmental fluctuations affect selection. Here we generalize quantitative genetic theory for environmentally induced stochastic selection to include general forms of frequency‐ and density‐dependent selection. We show how the relevant fitness measure under stochastic selection relates to Fisher's fundamental theorem of natural selection, and present a general class of models in which density regulation acts through total use of resources rather than just population size. In this model, there is a constant adaptive topography for expected evolution, and the function maximized in the long run is the expected factor restricting population growth. This allows us to generalize several previous results and to explain why apparently “ K ‐selected” species with slow life histories often have low carrying capacities. Our joint analysis of density‐ and frequency‐dependent selection reveals more clearly the relationship between population dynamics and phenotypic evolution, enabling a broader range of eco‐evolutionary analyses of some of the most interesting problems in evolution in the face of environmental variation.

One of the great challenges in evolutionary biology is to understand the reciprocal feedbacks between ecology and evolution to improve our understanding of factors affecting the capacity of organisms to adapt to environmental changes (Govaert et al. 2019). This requires that ecologically realistic assumptions are included in models describing the evolutionary dynamics. Two such general characteristics are important. First, fluctuations in population size reduce mean fitness of the population, and hence may induce density-dependent selection . Second, random fluctuations in the environment affect the growth rate and variability of most natural populations (Lewontin and Cohen 1969;Lande et al. 2003) and may strongly affect the capacity to produce adaptive responses to variation in environmental conditions (Dempster 1955;Gillespie 1973;Felsenstein 1976). Here we develop a general theory for how variation among phenotypes in their sensitivity of fitness to changes in the environment will affect evolution of density regulated populations.
A general framework for analyses of phenotypic evolution was provided by Lande (1976Lande ( , 1979 and Lande and Arnold (1983) who developed the deterministic quantitative genetic theory for evolutionary responses to selection of phenotypes following a multivariate normal distribution. Lande (1976) also considered stochasticity by random genetic drift causing stochastic changes in the mean phenotypesz with variance inversely proportional to population size N. This theoretical framework improved our conceptual understanding of phenotypic evolution (Lande 1979(Lande , 1982 and facilitated the development of statistical methods to analyze selection on quantitative traits (Lande and Arnold 1983) that soon became widely used (Kingsolver et al. 2001;Hendry 2017). However, these models ignored important ecological processes such as density dependence and stochastic fluctuations in the environment.
An important advance in our understanding of how organisms adapt to a fluctuating environment was provided by the results of Dempster (1955) who showed that the final outcome of selection on a single locus was determined by the mean geometric fitness of a genotype (Gillespie 1973(Gillespie , 1974Felsenstein 1976). In continuous time, this corresponds to the long-run growth rate as measure of fitness (Turelli 1981;Tuljapurkar 1982;Orzack and Tuljapurkar 1989;Saether and Engen 2015). An important advance to include ecologically more realistic assumptions into models of phenotypic evolution was provided by Lande (2007) who developed a general model in which the fitnesses of different individuals are differently affected by environmental fluctuations. An implication of this theoretical framework is that it enables the analyses of how environmental stochasticity generates stochastic fluctuations in selection, as well as in the response to selection. Such stochastic variation will be present even in extremely large populations, analogous to environmental stochasticity in population dynamics (Lande et al. 2003). Accordingly, even in large populations, variation in the mean phenotype over long time intervals usually contains considerable stochastic fluctuations (Hendry and Kinnison 1999;Kinnison and Hendry 2001;Hendry et al. 2008). This may be a result of fluctuations in the environment that affects the fitness of individuals differently, generating stochasticity in selection (see references in Hayward et al. 2018;Hunter et al. 2018).
Stochastic fluctuations in the mean phenotype of a population, especially when fitness is frequency dependent, can impact both ecological and evolutionary dynamics. Such frequency-dependent selection in which fitness of a phenotype or genotype depends on its frequency in the population (Ayala and Campbell 1974) is known to strongly affect phenotypic evolution (see references in Svensson and Connallon 2018). Famous examples are Fisher's runaway process of sexual selection (Fisher 1930, Lande 1981 and the evolution of stable sex ratios (Fisher 1930). Frequency-dependent selection affects phenotypic evolution if fitness functions are of the type W (z;z) for a univariate character with phenotype z in a population with mean phenotypes z (Lande 1976; see Abrams 2019 for ecological examples). Because the fitness function depends on the mean phenotype in the population, this model includes frequency-dependent selection. The dependency of the mean fitness on the mean phenotype then occurs in two ways, through thez in the fitness function as well as thez that appears by averaging over the distribution of z when computing the mean fitnessW (z) (Abrams et al. 1993a(Abrams et al. , 1993b. In this case, Lande (1976) showed that selection is not only determined by d lnW (z)/dz but also depends on the frequency dependence introduced by thez inW (z;z).
Another implication of the mean phenotype affecting the individual fitness function is that it has to be considered as a part of the environment. For instance, a large body of theory has been developed to study how the social environment affects the fitness of individuals focusing on how social interactions influence the strength and response to selection (Frank 1997(Frank , 1998. A key type of environmental fluctuation that is an important component of the feedback between ecology and evolution is fluctuations in population density. Furthermore, game theory has been widely applied (see references in Metz and Geritz 2016) to determine evolutionary resting points expressed by genotypes that cannot be invaded by any other types (i.e., the evolutionary stable strategy). In these models, fitness payoffs are often frequency dependent (Maynard-Smith and Price 1973;Abrams et al. 1993a;McGill and Brown 2007). However, these approaches generally ignore environmental stochasticity and rarely consider the role of population dynamics. A key connection between frequency dependence, evolution in fluctuating environments, and population dynamics relates to what Fisher (1930) in general referred to as "deterioration of the environment." Fisher's fundamental theorem of natural selection states that the rate of increase in mean fitness due to natural selection is exactly the additive genetic variance of fitness. Under frequency-dependent selection, the realized change in mean fitness also has a component due to the deterioration of the environment represented by the change in mean phenotype caused by natural selection (Frank and Slatkin 1992). In addition, effects of environmental stochasticity that differ among phenotypes will result in changes in mean phenotype, and under frequency-dependent selection this will affect the expected evolution (Lande 2007). Hence, frequency-dependent selection will not only affect the mean phenotype in the population but also affect the mean fitness, and thus population size.
A key component of the feedbacks between ecology and evolution is fluctuations in population density. If variation in population size N (or population density) affects fitnesses among genotypes or phenotypes differently, selection is said to be density-dependent ). In a seminal paper, MacArthur (1962) showed that in a constant environment evolution maximizes equilibrium population size or the carrying capacity K. In that case, changes in N represent an additional source of environmental deterioration. An important implication of these theoretical results was MacArthur and Wilson's (1967) introduction of the concept of r-and K-selection, based on the assumption that genotypes or phenotypes with high fitness (a high r) at small population densities were inferior to those with lower r at population sizes close to K (for reviews of this concept, see Roughgarden 1979;Boyce 1984;Reznick et al. 2002).
If selection is weak so that the mean phenotype is practically constant over long time intervals, then stochastic fluctuations in N will on average affect selection in a way depending on the mean phenotypez, so that density-dependent selection accumulates to frequency-dependent selection through time, as pointed out by Smouse (1976) and Heino et al. (1998). The response to selection over a short time interval may still depend on N, but over large intervals with an approximately constant mean phenotype due to weak selection, this dependence on density acts only through properties of the stationary distribution of N depending onz and not on N itself. Examples of this approach based on averaging over fluctuations in population size are the models on fluctuating r-and K-selection in a stochastic environment analyzed by Lande et al. (2009) and Engen et al. (2013).
One central concept in evolutionary biology is Wright's adaptive topography that may be defined either in the allele frequency space (Wright 1932) or in the phenotypic space (Simpson 1944). In deterministic quantitative genetics with no densityor frequency-dependent selection, the adaptive topography is the mean fitness as a function of the mean phenotype, or the mean Malthusian fitness in continuous time (Lande 1976;Lande and Arnold 1983). The applicability of this result when there are stochastic fluctuations in the environment has been questioned because there is no constant topography so that the population mean phenotypes constantly pursue a variable optimum (Lande 2008;Chevin et al. 2015). However, as pointed out by Lande (2007), there may be a constant topography that applies to the expected evolution, that is, the evolution averaged over the stochastic environments. Engen et al. (2013) demonstrated, under density dependence and stochastic fluctuations in population size, that this would lead to stationary fluctuations of mean phenotype around the phenotype maximizing the topography for expected evolution.
A central focus in early studies applying the concept of rand K-selection was to identify general characteristics of species that could be classified by fast growth under no density dependence (r-selected) or good performance at high population densities (K-selected). For example, Cody (1966) proposed that population densities of K-selected species should be higher than those of r-selected species because of more intense selection for more efficient resource utilization and successful intraspecific competition at population sizes close to the carrying capacity K. The dichotomy between r-and K-selected species was developed further by Pianka (1970) to associate a number of population dynamical characteristics and life-history attributes with the assumed relative strength of r-versus K-selection. Stearns (1977) and Boyce (1984) criticized such classifications of species, arguing that several other processes than r-and K-selection could generate similar covariation among life-history traits. Although the classifications of species by Pianka (1970Pianka ( , 1972) was shown to be independent of the strength of r-versus K-selection, these early studies still identified some general processes that may affect evolution of specific features as a response to densitydependent selection. Several of the proposed patterns, for example, related to the type of intraspecific competition, are based on assumptions about how the fitness of an individual depends on the mean phenotype in the population. For example, Southwood (1977) suggested that the intraspecific competition in rselected species should be of "scramble" type varying in intensity with fluctuations in the availability of resources, whereas in Kselected species the competition was of "contest" type for access to space or some fixed amount of a limited resource, for example, through territorial defense. Thus, the proposed influence of r-and K-selection on general ecological patterns was based on a combination of assumptions about how the strength and shape of density regulation was affected by the mean phenotype of the population. Accordingly, Gill (1974) suggested the term α-selection to denote selection for competitive ability itself, dependent on the distribution of phenotypes in the population. This provides an early attempt to link density-and frequency-dependent selection.
Three major results for evolution in continuous time were derived by Lande (2007). First, the individual mean fitness under environmental stochasticity was shown to be the deterministic growth rate of the individual minus the covariance between its growth and that of the total population. We will show how this measure of fitness relates to Fisher's fundamental theorem of natural selection. Second, the adaptive topography for the expected selection, without density-or frequency-dependent selection other than that generated by stochastic selection, is the longrun growth rate of the population expressed as function of the mean phenotype. Finally, the expected change in the mean phenotype, the response to selection, is the additive genetic covariance matrix G for z multiplied by the gradient of the long-run growth rate, generalizing the classical deterministic gradient formula of Lande (1979Lande ( , 1982. However, all these results assume the same effect of N on all phenotypes; thus, only ensuring a finite population size but not including any density-dependent selection . Here we extend the model of Lande (2007) to include general forms of density-and frequency-dependent selection. We use this approach to evaluate a rather general class of models, with density dependence depending on the mean phenotype, which gives a constant adaptive topography for the expected evolution. We show that the parameter maximized by evolution is the expected value of the quantity defining the density regulation by limiting the population growth. In nongenetic models this is N itself, but allowing for genetic variability and evolution, this quantity may in general depend onz.
We use this general model to explore how variation among phenotypes in their sensitivity to changes in the environment affects phenotypic evolution in populations subject to fluctuating r-and K-selection. Density regulation is not only allowed to act through population size or density itself but also through other quantities such as functions of total biomass of the individuals in the population. In addition, α-selection may affect the use or access to resources dependent on the distribution of the phenotypes in the population (Gill 1974;Johsi et al. 1998). We show how joint simulations of density N and mean phenotypez can be performed numerically in complex models where analytical expressions for the mean Malthusian fitness are not available.
To illustrate the practical application of this general theory, we show that it can be used to classify species along a continuum ranging from fast pace-of-life species at one end to slow ones at the other end (see Wright et al. 2019). We do this by considering how evolution of body mass is mediated by specific assumptions about the effects of frequency-dependent selection and the form of density regulation. Body mass is a phenotypic character that affects a magnitude of ecological processes (McNab 1971;Peters 1983) and influences the pattern of density regulation in many populations (Lopez-Sepulcre and Kokko 2005). It is also a trait that for more than two centuries has been known to show distinct variation along environmental gradients (Bergmann 1847;Gaston et al. 2008) and to respond rapidly to environmental changes (Ozgul et al. 2009;Sheridan and Bickford 2011). These effects of body mass variation generate interactions between frequencyand density-dependent selection, which may produce relationships between the position of species along the pace-of-life continuum and their abundances expressed by the carrying capacity K. We show that the form of these relationships depends on specific assumptions about how environmental variation affects individuals of different size (Gardner et al. 2011) and how the distribution of body masses in the population influences its dynamics. Lande (2007) performed a general analysis of the stochastic component of selection for a large population in continuous time, essentially a stochastic version of the evolutionary theory for multinormally distributed phenotypes of Lande and Arnold (1983). The population size is assumed to be large enough for genetic drift to be neglected so that all stochasticity is generated by stationary environmental fluctuations. The vital rates have small or moderate fluctuations so that the long-run growth rate gives accurate demographic projections. The theory assumes weak densityindependent selection so that the Malthusian fitness provides an accurate description in continuous time that approximates natural selection in discrete time (Charlesworth 1994). The phenotypes and breeding values are assumed to follow approximately multivariate normal distributions with phenotypic and genetic covariance matrices P and G, respectively. These covariance matrices are assumed to be unaffected by environmental variation and stay approximately constant by mutation and recombination (Lande and Arnold 1983), which requires that selection as well as linkage among loci with epistatic fitness interactions are weak (Bulmer 1971(Bulmer , 1976Lande 1982). Furthermore, most results presented here describe long-term fluctuations of mean phenotype around some central value, so that approximating G and P by their average values over these fluctuations is a reasonable approach.

BASIC DEFINITIONS
Diffusion processes are stochastic processes that are continuous in state and time and can be used as approximations for discrete processes under small or moderate fluctuations. The noise in a diffusion is white, implying that we assume stationary environmental fluctuations with no temporal autocorrelation. Diffusions are characterized by their infinitesimal mean and variance that can be chosen as the mean and variance of the changes in the underlying discrete time model during a time unit (Karlin and Taylor 1981). In multivariate processes, infinitesimal covariances are also required. Let n(z)dz be the number of individuals with a d-dimensional phenotype in the interval (z, z + dz) so that the total population size is N = n(z)dz and the distribution of phenotypes is p(z) = n(z)/N (for a summary of the description of the parameters, see Table 1). The infinitesimal mean of n(z), which is the average growth rate of the population in continuous time, is expressed as [r(z) − g(N )]n(z), where g(N ) is a fixed function determining the density regulation. In this model, there is no density-dependent selection because all phenotypes are equally affected by N, meaning that densities of all phenotypes are reduced by the same fraction by the density regulating term .
In a fluctuating environment the fitness of an individual can be expressed as a reaction norm where both fitness in the average environment and the effect of the environment on fitness will depend on phenotype. To exemplify, assume that r(z) = r 0 (z) + σ(z)ε as illustrated in Figure 1, where ε is a standardized environmental variable with zero mean and unit variance. When σ(z) varies among phenotypes, this is an example of an interaction between genotype and environment (G × E ). Then, considering two phenotypes y and z, it appears that cov[r(y), r(z)] = c(y, z) = σ(y)σ(z) and corr[r(y), r(z)] = 1, if σ(z) has the same sign for all values of z. More generally, the noise for different phenotypes may have some correlation structure ρ(y, z), giving c(y, z) = σ(y)σ(z)ρ(y, z). Another class of models appears when there are several environmental variables affecting the growth rates. If the Malthusian fitness, for example, depends linearly on two different environmental variables, r(z) = r 0 (z) + σ 1 (z)ε 1 + σ 2 (z)ε 2 , then the covariance function becomes more complicated and the correlation will in general be smaller than 1.

r(z)
Malthusian growth rate in the average environment for phenotype z, written as r(z; N,z) if vital rates depend on population size and mean phenotype r 0 (z) Malthusian growth rate at small population densities for phenotype z r(z) Malthusian growth rate of the population in the average environment r(z) Long-run growth rate of the population with mean phenotypez averaged over environmental states m(z,z) Expected fitness of an individual with phenotype z in a population with mean phenotypez W(z) Fitness of an individual with phenotype z , written as W (z,z) if it depends on mean phenotypez W (z) Mean fitness of individuals with mean phenotypezin the average environment z * Value of the phenotype at which the maximum of the adaptive topography is located σ 2 d Demographic variance in the population dynamics due to random variation in individual fitness other than that determined by phenotypes (e.g., binomial variation in survival) σ 2 (z) Sensitivity of individual growth rates to environmental variation σ 2 e (z) Environmental variance in the dynamics of total population size, defined as the variance in per capita growth rate as influenced by the stochastic environment and evaluated for the population with mean phenotypez σ 0 Parameter expressing the variability in the environment c (z,z) Covariance between growth rates of z and that of the total population for a given mean phenotypez c (y, z) Covariance between growth rates (increments) of phenotypes y and z over environmental states, written as c(y, z; N,z) if vital rates depend on population size and mean phenotype μ(z) Vector of covariances cov[σ(z), Variable describing environmental variable i, standardized with mean 0, and unit variance n(z)d(z) The number of individuals with d-dimensional phenotype in the interval (z, z + dz) N Total population size so that N = n(z)dz p(z)

Distribution of phenotypes p(z) = n(z)/N g(N)
A strictly increasing fixed function of N defining the density regulation γ(z) Sensitivity of the individual growth rate of individuals with phenotype z to fluctuations in N α, β Intercept and slope, respectively, describing linear effects of z on r(z) η Coefficient describing how σ(z) depends on z a, b Intercept and slope, respectively, describing linear effects of z on ln γ(z)

R(z)
Function determining how much an individual with phenotype z contributes to within-species competition for resources The theory of Lande (2007) is based on environmental effects being described by a general covariance function c(y, z). Assuming that fluctuations in the environment can be described as a white noise process in continuous time and using the fact that the temporal change in n(z) during time dt is n(z)r(z)dt, the infinitesimal covariance between the processes n(y) and n(z), for any phenotypes y and z, is c(y, z)n(y)n(z).
The infinitesimal mean and variance of N are then mean and variances over phenotypes, [r(z) − g(N )]N and σ 2 e (z)N 2 , respectively, which by the continuous version of the mean and variance of a sum arē The dependence on mean phenotypez = zp(z)dz is due to the fact that the normal distribution p(z) has meanz, which is omitted to simplify the notation. The parameter σ 2 e (z) is the environmental variance in the dynamics of N. In finite populations, there will always also be a demographic variance σ 2 d due to random variation in individual fitness other than that determined by phenotypes, such as binomial survival. Demographic noise gives no contribution to the covariance function c(y, z), but it affects the dynamics of N by adding σ 2 d /N to the infinitesimal variance of ln N and subtracting σ 2 d /(2N ) to the infinitesimal mean. Here we assume that the population is large enough for the demographic noise to be ignored.

EXPECTED FITNESS IN A STOCHASTIC
ENVIRONMENT Lande (2007) showed that the mean selection differential over environmental states is is called the expected fitness of an individual with phenotype z. Except for a trivial additive constant not affecting selection, this is its Malthusian fitness minus the covariance between its growth rate and the growth rate of the total population. More precisely, the quantity subtracted from the deterministic growth rate r(z) in equation (2) can be expressed as which is the infinitesimal covariance between ln n(z) and ln N. A related result based on selection coefficients on alleles in a single locus haploid model was given by Gillespie (1991).

ADAPTIVE TOPOGRAPHY AND GRADIENT FORMULA IN A STOCHASTIC ENVIRONMENT
The long-run growth rate is an important concept in population dynamics and population genetics, expressing the growth rate of log population size over large time intervals in a stochastic environment in the absence of density regulation (Saether and Engen 2015). Using diffusions, this is the infinitesimal mean of ln N which in the present model, by the transformation formula for diffusions, isr(z) =r(z) − 1 2 σ 2 e (z). In this stochastic model, the response to selection depends on the environment, but the longterm evolution is mainly determined by the mean response over environmental fluctuations. Using the expectation symbol E for this mean , a major result of Lande (2007) is the generalized gradient formula where ∇ denotes the operator (∂/∂z 1 , ∂/∂z 2 , . . . , ∂/∂z d ), showing that the long-run growth rate defines an adaptive topography for expected evolution. Note that this topography is different from the one used in deterministic theory, which expresses exactly the response to selection at a given mean phenotype. In addition, it also has a different interpretation. In a stochastic model, the response is determined by the stochastic environment, but there is still a topography determining the expected response. If this topography has a maximum at z * , evolution will on average attempt to reach this value, but due to environmental fluctuations the mean phenotype will end up in stationary fluctuations with z * as a center point, as demonstrated by Engen et al. (2013).
Lande (2007) also derived the infinitesimal covariance matrix for the selection differential vector and the response, as well as the infinitesimal covariance between log population size and the mean phenotype components. These results are given in Appendix A, also showing more details concerning the derivations.

ENVIRONMENTAL NOISE
In general, the covariance function c(y, z) may be complicated, but we obtain interesting results even for the simplified model with linear environmental effects and G × E interactions, as illustrated in Figure 1. For example, in shorebirds such as the Eurasian oystercatcther Haematopus ostralegus cold winters affect all birds negatively (Goss-Custard 1985), but small individuals with smaller intake of food are generally more likely to die than larger birds (Stillman and Goss-Custard 2010). With these assump- Writing v(z) = GP −1 μ(z) and applying results given in Appendix A, the infinitesimal second moments are  From this it appears that the absolute value of the correlation between d ln N and any component dz i is 1, the signs being the same as those of μ i (z). If σ(z) is a decreasing or increasing function of the ith phenotype component z i , then corr[d ln N, dz i ] equals −1 or 1, respectively. The dz i are also fully correlated by Figure 2 illustrates this model, using a single trait and assuming r(z) = α − βz and σ(z) = σ 0 e −ηz . The adaptive topography is thenr(z) = α − βz − 1 2 σ 2 0 e −2ηz+η 2 P with a maximum at z * = ln(ησ 2 0 /β)/(2η) + ηP/2. One plausible interpretation of this model is to consider z as body size, so that the growth rate as well as the effects of environmental noise decrease with increasing size. Lande's (2007) given by equation (2) and scaled by an additive constant to be zero at z = z * (Fig. 2), then depends strongly onz so that frequency-dependent selection here is generated solely by this simple stochastic term.

NATURAL SELECTION
An important implication of the expression for the expected fitness (eq. 2) is that it relates directly to Fisher's (1930) fundamental theorem of natural selection (see derivations in Appendix B). In the stochastic model, the rate of increase in the long-run growth rater(z) due to natural selection is the additive genetic variance of the expected fitnessm given by equation (2). However, if the covariance function c(y, z) is not constant then c(z,z) always depends onz. Therefore, except for the trivial case in which environmental fluctuations affect the growth rate of all phenotypes equally, stochasticity in selection always implies frequency-dependent selection. As a consequence, there is also an additive component to the expected change in the long-run growth rate due to the change inz (expressed by eq. B2), caused by natural selection (for more details, see Appendix B).

GENERAL RESULTS
We can generalize this model for stochastic selection in three ways, explained in some detail in Appendix C. The results will mainly be as above, except for the gradient formula expressed by equation (3), which is not directly applicable as first demonstrated in the one-dimensional deterministic case by Lande (1976). However, it can still be used by appropriate interpretation of derivatives. Assuming, as in Lande (2007), that the environment affects different phenotypes differently, the three generalizations are: (1) introducing density-dependent selection by including a general form of density regulation, (2) introducing frequency-dependent selection by allowing any parameter expressing the vital rates of individuals not only to depend on the individual phenotype but also on the mean phenotype in the population, and (3) introducing density-and frequency dependence in the stochastic effects by allowing the covariances between n(y) and n(z) to depend on both N andz. The generalization (1) was carried out in Appendix A of Engen et al. (2013). We will consider biologically important examples of the second generalization by assuming individual heterogeneity and explicitly modeling how phenotypes mediate density regulation. Because we are relying on Lande's (1976Lande's ( , 1979Lande's ( , 1982Lande's ( , 2007 approach assuming that G and P are approximately constant, the distribution of allele frequencies depends only onz, implying that the model represents a general form of frequency-dependent selection within this framework. The empirical foundation for the third of our generalizations is that environmental variance may change with population size (Gamelon et al. 2017;Hansen et al. 2019), which then must also affect the covariances c(y, z). We thus generalize to the situation in which c(y, z) also depends on population size N and the mean phenotypez.
Because (N,z), according to Lande (2007), is described as a diffusion process, the infinitesimal moments defining it are evaluated for a given value of (N,z) (and actually n(z) for all values of z), which is the state of the diffusion in accordance with standard description of multivariate diffusions (Karlin and Taylor 1981). It is shown in Appendix C that we can consider general models where the growth rates of n(z) for a given (N,z) is r(z; N,z) and the covariances are (y, z; N,z).
All the results derived by Lande (2007), up to the gradient formula, are then valid for this generalized model just by including the way on which the parameters depend on N andz (see Appendix C for details). In particular, the mean Malthusian fitness given by equation (2)  Here the firstz inc appears from the distribution of phenotypes p(y) having meanz, whereas the second represents the assumption that the covariance between dn(y) and dn(z) may depend on the mean phenotype in the population (frequency-dependent stochasticity).
When it comes to the gradient formula, great care must be taken with respect to interpretation. The growth rate of ln N is nowr N, x).
Here, x =z represents how the infinitesimal moments may depend on the mean phenotype, andz in the formula arises from integrations of the distributions p(z) and p(y) that also depend on the meanz. As a consequence, the gradient formula is still valid, provided that the gradient (the derivatives) operate only on z keeping x constant. This can be expressed mathematically as The message from this derivation is that the gradient formula without this conditioning gives incorrect results if the distribution of individual vital rates depend on the mean phenotype in the population. We will demonstrate the application of this formula below, showing examples of more general types of density regulation.
In practice, it may often be impossible to find an analytical expression forr(z; N, x)| x=z . If so, then it is preferable in numerical computations to use the expression ∇ p(z;z) = P −1 (z − z)p(z;z) arising from taking the derivative of the normal distribution with respect to its mean (Lande and Arnold 1983) as in the derivation of equation (3), giving To perform simulations of (z, ln N ), this integration would have to be performed numerically at each discrete step used to approximate the continuous model.

POPULATION SIZE N
A constant adaptive topography is unlikely to result from any general form of density-and frequency-dependent selection. Although equation (4) is expressed as a gradient formula, it is not an adaptive topography due to the dependence of x =z. Under weak selection, evolution during relatively short time intervals may still approximately be defined by equation (4) and considered as an adaptive topography depending onz (and not x), but over the long term the topography changes with the mean phenotype. However, for a special class of models, analyzed by Engen et al. (2013), that includes density dependence covering many ecologically realistic situations, a constant adaptive topography still appears for the expected evolution even under density-and frequency-dependent selection.
Such models have the form r is an increasing function with g(0) = 0 so thatr 0 (z) is the deterministic growth rate under no density regulation. The function γ(z) describes situations where growth rates of different phenotypes z are differently affected by density. Saether et al. (2016) have previously used this model to analyze the effects of density-dependent selection on clutch size evolution in a small temperate passerine, the great tit Parus major. Engen et al. (2013) showed for this model, assuming weak selection and constant environmental variance, that expected evolution will always be in a direction increasing the function Eg(N ) = [r(z) − σ 2 e (z)/2]/γ(z). This is the mean of g(N ) over fluctuations in N for a givenz that are driven by fluctuations in the environment. Accordingly, this is a mean value over environmental fluctuations over a time interval wherez varies little due to the assumption of weak selection. Provided that the presence of trade-offs between population growth rates and effects of density prevents an indefinite increase in fitness, Eg(N ) will take its maximum for some mean phenotype z * . However, this maximum will be reached only in the deterministic case with σ 2 e = 0, while in generalz will fluctuate stationary around the center z * (which is not necessarily the mean value). Hence, Eg(N ) defines an adaptive topography for expected evolution, while in practice evolution is disturbed by environmental stochasticity. Note that environmental stochasticity affects both the center value z * and the degree of evolutionary fluctuations ofz around it, in analogy with well-known results in nongenetic population dynamics where σ 2 e has been known to affect the mean population size as well as the magnitude of its fluctuations. In the deterministic case, z * maximizes g(K ) and hence the carrying capacity K, a result going back to MacArthur (1962).
We can illustrate this model by a linearization using a onedimensional phenotype z. Because r 0 may take any value while γ must be positive to express density regulation, the natural linearized model is the first-order Taylor approximation to r 0 and ln(γ), r 0 (z) = α − βz, ln γ(z) = a − bz. With positive values of β and b, we may think of z as a character that expresses the size of an individual. Then, larger individuals have smaller population growth rates at small densities, but they are better competitors at higher densities by smaller values of γ making them less affected by density. If z is normally distributed with meanz and phenotypic variance P, then the mean parameters arer 0 (z) = α − βz andγ(z) = e a−bz+b 2 P/2 . With a constant environmental variance σ 2 e , meaning that c(y, z) = σ 2 e , the mean phenotype maximizing Q(z) = [r 0 (z) − σ 2 e /2]/γ(z) then takes the value Two general types of such responses have been demonstrated. First, small individuals may be favored at low population sizes but large ones favored close to K. For example, Trinidad guppies (Poecilia reticulata) from high-predation environments are younger and smaller at sexual maturity, invest more resources to offspring production and produce more offspring per brood compared to individuals from low-predation environments (Reznick and Bryga 1996) . Furthermore, analyses of long time series of demographic data, as well as experimental manipulations, indicate that larger guppies subject to low-predation pressures are less affected by high densities than the highly reproductive individuals exposed to high predation pressure (Bassar et al. 2013). This provides strong evidence for genetic differences in responses to fluctuations in population density strongly affecting life-history evolution in this species (Reznick et al. 2019). Accordingly, in this kind of system, our model with positive β and b predicts that increasing densities should favor larger individuals (Fig. 3A). It also appears that selection for body size is frequency dependent because the fitness functions depend strongly on the mean phenotype (Fig. 3B,D), but with no clear directional pattern. Second, an opposite relationship between z and increasing population size (β < 0, b < 0) is also possible (Fig. 3C). In temperate cervides, large body mass is often associated with high reproductive output (Saether and Haagenrud 1983;Albon et al. 1987), which is determined by access to high-quality food during summer (Saether and Heim 1993;Langvatn et al. 1996). However, large individuals are susceptible to high population densities especially during winter because they require greater amounts of food (Andersen and Saether 1992). In this case, density dependence favors smaller individuals (Fig. 3C,D) and generates negative frequency dependence for body size.
An important extension to include more ecological realism is that the environmental variance may depend on the phenotype.
If the sensitivity to environmental fluctuations increases with the value of the phenotype, the environmental variance at equilibrium increases for a given variability in the environment σ 2 0 (Fig. 4A) and strongly affects the z * (Fig. 4B).

INDIVIDUALS' USE OF RESOURCES
In addition to being linear with respect to g(N ), this model is based on the assumption that density regulation acts through the population size N, which is realistic if all individuals have practically the same resource requirement. However, variation in phenotype z may represent variation in body size or other characteristics affecting individual use of resources and hence their impact on density-regulation reducing population growth rate. This may occur if resource use scales allometrically with body mass, as in many species of ungulates, resulting in density effects on population growth rate operating through herbivore biomass (see Fig. 13.4 in Owen-Smith 2002). In standard models, each individual contributes the same quantity 1 to the regulation acting through N (see Travis et al. 2013). More generally, let R(z) be a function determining how much an individual with phenotype z contributes to within-species competition for resources, which underlies density dependence in population growth (MacArthur 1972;Roughgarden 1976;Dieckmann and Doebeli 1999). The density regulation may then act through the sum of these quantities over all individuals,

R(z)n(z)dz = N R(z)p(z)dz = NR(z).
Applying this to the above model then gives a deterministic growth rate of the form r(z; N,z) = r 0 (z) − γ(z)g(NR(z)).
In a stochastic environment, the growth rate of ln N at a given N is theñ so that the expected value of g(NR(z)) must bes(z)/γ(z) = Q(z), wheres(z) = r 0 (z) − σ 2 e (z)/2. Note that the resource use has no effect on the fitness averaged over N which is always zero. Assuming weak selection and applying equation (4), we show in Appendix D that Because G is a positive definite matrix, this time derivative is nonnegative, meaning that the expected evolution ofz is always in a direction of increasing Q(z) and is therefore an adaptive topography for expected evolution. As in the model of Engen et al. (2013), this means that expected evolution is always in a direction approaching the mean phenotype z * maximizing Q(z), while stochastic environments generating fluctuations in N tend to move the mean phenotype away from this value, resulting in stationary fluctuations around z * (see Engen et al. 2013 for details). In conclusion, within the above modeling framework evolution tends to, on average, maximize the expected value of the factor g(NR(z)) regulating the growth rate of the population, resulting in an adaptive topography Q(z) for expected evolution even in the presence of density-and frequency-dependent selection. Note that the solution for z * and the mean value of g(NR(z * )) will be the same for any function g(NR(z)), but the mean of N or the carrying capacity may depend strongly on the choice of that function. If, for example, g is the identity function so that the dynamics are logistic and we still use the assumption of week selection so that the covariance between N andz is practically zero, implying that E[NR(z)] = E(N )R(z) for a given mean phenotype, then Hence, the expected population size is proportional to 1/R(z). At the evolutionary equilibrium the carrying capacity K, defined as the population size at which the deterministic growth rate is zero, is accordinglyr(z * )/[γ(z * )R(z * )], so that the expected population size is somewhat smaller than K, in this model is given by Note that the fitness of individuals with different phenotypes z may still react differently to density expressed through g(NR(z)) due to individual variation in the parameter γ(z).
If the function g is linear, this means that the average total resource requirement ENR(z) is maximized by evolution, and if individual use of resources is independent of z (frequencyindependent selection), then the expected population size is maximized, as in the models of Lande et al. (2009) and Engen et al. (2013). At the deterministic limit (as σ 2 e approaches zero) and when N has a constant equilibrium value, N = K, the constant use KR(z) of resources is maximized, and finally, if in addition R(z) is constant, then K is maximized in accordance with the classical result of MacArthur (1962).

PHENOTYPE-SPECIFIC RESOURCE USE
For a general resource function R(z) and density dependence expressed by g(NR(z)) we then find, using the previous linearized models for r(z) and ln γ(z) and the general results in the previous section, that Figure 5A illustrates how r(z * ) and γ(z * ) may depend on the environmental variance in this model. If the phenotype z here is considered as some measure of log body size of individuals so that the actual size is w = e z , a class of reasonable resource functions is obtained by choosing R(z) proportional to w θ = e θz . In this model, the log 10 of the expected population size, as well as carrying capacity, depend on the environmental variance and how the resource function varies with z. When the R(z) is constant or weakly related to z (θ close to zero), then log 10 of the expected population size as well as the carrying capacity decrease asr(z) increases with increasing environmental variance (Fig. 5B). This agrees with previous results for models with constant R(z) (Lande et al. 2009;Engen et al. 2013), where increased environmental noise leads to r-selection and smaller K (see also Roughgarden 1971). However, for larger values of θ, indicating thatz has a stronger effect on body size and use of resources, population sizes increase together withr(z) as σ 2 e increases. These results for large values of θ may be part of the explanation for why typically r-selected species are usually small and exist in populations with large numbers of individuals, and why K-selected species are large with a small numbers of individuals.

Discussion
An important focus in evolutionary ecology is the inclusion of realistic ecological assumptions in the models describing evolutionary processes (Govaert et al. 2019). As an example of such a close interaction, there is even in age-structured models a clear mathematical link between the demographic variance in population dynamics and random genetic drift expressed by the variance effective population size (Engen et al. 2005a(Engen et al. , 2010Olsson et al. 2013;Myhre et al. 2016;Trask et al. 2017). This is due to the fact that demographic variance and genetic drift arise from stochastic contributions that are independent among individuals. Here we have developed another link between ecology and evolution by showing that environmental stochasticity expressed by the environmental variance, which also is defined for age-structured models (Tuljapurkar 1990;Engen et al. 2005b), is closely related to the environmental stochasticity in selection driven by phenotypes being differently affected by environmental fluctuations (Gillespie 1973(Gillespie , 1974Bull 1987). In the quantitative genetic theory of Lande (2007), this is expressed by the environmental covariance function c(y, z). It has also previously been shown that the theory of fluctuating r-and K-selection developed for nonoverlapping generations (Lande et al. 2009;Engen et al. 2013) is applicable under small and moderate fluctuations in age structure (Lande et al. 2017), indicating that most of the conclusions presented here may be valid under less restrictive assumptions. Thus, we suggest that our extensions of Lande (2007) provide a general theoretical framework for analysing eco-evolutionary dynamics under a wide range of ecologically realistic conditions. This represents an alternative approach to the long tradition of analysing fluctuating selection in population genetic models (e.g., Felsenstein 1976;Gillespie 1991;Ellner and Hairston 1994;Turelli et al. 2001).
A central concept in evolutionary theory is Wright's (1932) adaptive topography (Svensson and Calsbeek 2012). During the last decade, there has been much focus on adaptive topographies that fluctuate in time so that evolution continuously chases a changing optimum, but assuming a constant environmental covariance c(y, z) = σ 2 e (Chevin et al. 2015 and no density-dependent selection. These are examples of environmentally induced stochastic variation in selection due to fluctuations in the fitness function. Lande (2007) considered fluctuations in the environment in a more general way using the covariance function c(y, z) and introduced the concept of adaptive topography for expected evolution rather than evolution realized at any specific time (see also Starrfelt and Kokko 2012). This defines a constant adaptive topography even under stochastic evolution that describes the long-run evolution under stochastic fluctuations. The general result that the long-run growth rater(z) is the relevant adaptive topography in such models (Lande 2007) was also confirmed by Lande (2008) dealing with allele frequencies at one or two loci. The model was exemplified by derivation of the exact stationary distribution of allele frequencies in the Haldane-Jayakar (1963) model in which the three genotypes in a biallelic model have fitnesses fluctuating in time but on average are equal. Here we have used another example, illustrated by Figure 2, showing how frequency-dependent selection may arise from differences in sensitivity to environmental noise through the covariance function c(y, z). This also illustrates Lande's measure of expected fitness given by equation (2), and emphasize the difference between the adaptive topography for expected evolution, and the fitness function that vary in time due variation in mean phenotype.
Plasticity is a basic concept in evolution expressing environmental effects on phenotypes (Tufto 2015). The parameters describing these environmental effects may be genetically determined so that plasticity may be under selection (Lande 2009(Lande , 2019Chevin et al. 2010;Tufto 2015;Chevin and Hoffmann 2017). As in the model of Lande (2007), the present model does not include plasticity in z, but there are phenotype-dependent environmental effects on fitness, as exemplified by our model where r(z) is a linear function of the environment with intercept and slope depending on z. For the effects of plasticity in z, see, for example, Lande (2019). Our generalized model with covariance functions c(y, z; N,z) that are frequency-and densitydependent describes complex environmental effects on fitness. Because these effects determine this covariance function, which changes under evolution ofz, the theory provides also a platform for analyzing evolution of environmental effects on fitness in a rather general way.
One important effect of fluctuating environments on the selection of different phenotypes is density regulation. The environmental effects in the theory of Lande (2007) are expressed as a white noise process with no temporal autocorrelations. Density, when it acts as an environmentally fluctuating effect on selection, may often involve large temporal autocorrelations, consistent with slow return times to the equilibrium. Under weak selection, the fluctuations in population size will be much faster than changes in mean phenotype so that averaging over the fluctuations in N provides an accurate description of the evolutionary process. We have seen that this approach transforms density-dependent selection into frequency-dependent selection (Fig. 2), establishing a formal link between them, as pointed out by Smouse (1976) and Heino et al. (1998). This frequency dependence comes in addition to the frequency dependence generated by the environmental covariance function, as shown by Lande (2007) and in Figure 3. Thus, our theoretical approach provides a general framework to analyze selection that includes phenotypespecific effects of fluctuations in the environment and population size, as well as the influence on selection of the distribution of phenotypes in the population (Figs. 2, 3). Figure 3 also illustrates an important point in relation to our comments on Fisher's fundamental theorem of natural selection under stochastic evolution. When the mean phenotype fluctuates, the fitness function defined by Lande (2007) may change considerably, and there will be corresponding fluctuations in the additive genetic variance of fitness. We have seen that this variance at any time is exactly the rate of change in the long-run growth rate due to natural selection. However, as Fisher (1930) pointed out, responses to selection, changes in population size, or external environmental changes, can all affect the fitness function and hence the goal of evolution. Such changes, which are examples of what Fisher called "deterioration of the environment," will al-ways introduce an additional component to the rate of change in the long-run growth rate under stochastic evolution (Frank and Slatkin 1992), here expressed by equation (B2).
Although selection favors increased fitness at each time under density-dependent selection, long-term fluctuations in population size will always make changes in the mean phenotype accumulate in some complex way , depending on the characteristics of the density regulation. One key question is whether there exists an adaptive topography that describes the final outcome of this long-term evolution. There are two classical results in deterministic theory relevant for this. MacArthur (1962), using a logistic type of model, showed that evolution always maximizes the carrying capacity K. A related result is that if population growth in an age-structured model is regulated by the density of individuals in one particular age class, then the density of this class is maximized (Charlesworth 1994). This age class was named "the critical age class." Engen and Saether (2016) showed that if all density regulation of vital rates in an age-structured model acts through the same function of the vector of population densities across ages or stages, then the value of this function is maximized by evolution under no stochasticity. Engen et al. (2013) and Saether (2016, 2017) also included density-dependent selection in a stochastic environment, assuming that the Malthusian fitness depends linearly on some function of population size. In this case, the topography for the expected evolution turned out to be the expected value of this function, a result closely related to the previous conclusions of MacArthur (1962), Charlesworth (1994), and Engen and Saether (2016). Here we have shown, using the same type of model, but with density-regulation acting through resource requirement by the function g(NR(z)) rather than g(N ), that this more general density-regulating function, also depending on mean phenotype, is an adaptive topography for expected evolution. However, in the most general case when mean vital rates are general functions of the mean phenotype, as well as population density, as for example in Soay sheep Ovies aries (Kentie et al. 2020) and Trinidad guppies (Renzick et al. 2019), we have seen that it is only possible to find the conditional gradient formula given by equation (4). This implies that there is no constant adaptive topography even for the expected evolution.
The presence of r-and K-selection involves the basic tradeoff that phenotypes with the highest fitness under small densities are those most strongly affected by population sizes close to the carrying capacity (Boyce 1984;Mueller 1997;Travis et al. 2013). Empirical evidence now exists both for natural populations and under experimental conditions, that such densitydependent selection may be important in affecting the capacity for adaptive changes to fluctuations in the environment (Mueller 1997;Travis et al. 2013;Saether et al. 2016;Kortessis and Chesson 2019), with adaptive topographies varying through time. Our generalized model with frequency-dependent selection based on Lande (2007) shows that if density regulation acts through total resource requirement depending onz, there is still a constant adaptive topography with important implications for such populations. In particular, such an effect of total biomass on the population dynamics seems to be a general feature in species and ecosystems (e.g., Caughley and Krebs 1983;de Roos and Persson 2013;Bassar et al. 2015). For this type of density regulation, this topography determines the expected evolution and defines a function of mean phenotype that is maximized through longterm evolution and not the stochastic realizations of response to selection at each time. By also taking into account the parameter γ(z) that is a factor in the term defining the density regulation, this makes it possible to evaluate how specific assumptions, concerning how different phenotypes are affected by the density regulation, will affect the expected evolution ofz. In particular, the way average body size affects density-regulation determines the direction of the expected evolution of body size (Fig. 3). This indicates that how N as well asz affect resource limitation can strongly influence the pattern of phenotypic evolution through density-dependent selection in natural populations (Abrams 2019) and distribute species along a slow-fast pace-oflife continuum (Wright et al. 2019(Wright et al. , 2020. Early studies suggested that the pattern of covariation among a wide array of species-specific characteristics should be associated with the relative strength of r-versus K-selection acting on the species. For example, Pianka (1970) proposed that Kselection should lead to large body size and high efficiency in resource use, whereas r-selection should produce small individuals favoring large output of new individuals, but simultaneously result in smaller densities of r-selected than of K-selected species at equilibrium population size (Cody 1966;Pianka 1972). Although such patterns are based on a large number of assumptions that make general comparisons across species or systems extremely difficult (Boyce 1984;Mueller 1997;Saether and Engen 2015), our results indicate that there exists a potential for establishing general patterns between phenotypic characteristics and patterns of covariation of life-history traits, dependent on how density regulation through resource limitation affects the population dynamics (e.g., Schoener 1976) as well as response toz (Abrams 2019).
Here we suggest that r-selection in general will result in higher abundances compared to K-selection (Fig. 5) if larger individuals possess a stronger effect on density regulation than small individuals. This corresponds to one of the general differences between r-and K-selected species suggested by Southwood (1977). Hence, population density tends to increase with decreasing body mass in a wide array of taxa (Damuth 1981(Damuth , 1987Peters 1983). Furthermore, if large individuals are more sensitive to environmental fluctuations this will result in evolution of smaller body size (Fig. 4B) and larger environmental stochasticity in the pop-ulation dynamics (Fig. 4A). Accordingly, Wright et al. (2019) argued that stochastic variation in population size may be an important factor in generating within-population variation among phenotypes in the tempo of individual life histories, which ultimately will lead to the co-evolution of life-history traits that distribute species along a slow-fast pace-of-life continuum.
Density-dependent selection represents the simplest form of feedback between the ecological effect of an organism's own making and the changes it causes in the selective environment. Evolutionary models of density-dependent selection in stochastic environments have allowed the study of feedbacks between ecology and evolution in more realistic scenarios (Abrams 2019;Govaert et al. 2019). Here we integrate density-dependent selection with frequency-dependent selection in stochastic environments, enabling us to further study how phenotypic variation interacts with population regulation in stochastic environments. Further integrating ecological and evolutionary process will increase our ability to predict the capacity of different organisms to adapt to increasing levels of stochastic variation in the environment.

AUTHOR CONTRIBUTIONS
Engen developed the theoretical results initiated by discussions with Wright, Araya-Ajoy and Saether on ecological correlates of r-and Kselection. Saether, Wright and Araya-Ajoy reviewed relevant biological literature and contributed substantially to the introduction and discussion. All 4 authors continuously discussed the development of the manuscript.

ACKNOWLEDGMENTS
The study was funded by the Research Council of Norway (SFF-III 223257/F50 and FRIKLIM 267511) and the Norwegian University of Science and Technology (NTNU). We are extremely grateful for the excellent comments by two reviewers on a previous version of the article.

DATA ARCHIVING
No data are included in the present study.

Environment
Using the transformation formulas for diffusions (Karlin and Taylor 1981) the infinitesimal moments of the log population size for the model of Lande (2007) given in the main text are the so-called long-run growth rate During time increment dt, the number n(z)dz, as well as the fraction p(z)dz of individuals in (z, z + dz), changes by selection.
Writing d p(z) for the change in p(z), the mean phenotypez after selection has changed by an amount dS = zd p(z)dz, known as the selection differential (denoted S(t ) by Lande (2007), where t is time, but we prefer dS(t ) because it is an infinitesimal increment of a vector). The general theory of Lande (1976Lande ( , 1979 and Lande and Arnold (1983) then gives the response to selection, also taking into account the transmission from parents to offspring, as dz = GP −1 dS(t ). Now, considering the two variables p(z) = n(z)/N and p(y) = n(y)/N, Lande (2007) used the transformation formulas for diffusions and the infinitesimal moments of the threedimensional diffusion (Karlin and Taylor 1981) n(z), n(y) and N to find their infinitesimal moments as and where dependence onz is not shown, meaning thatc(z) = c(z,z) = c(y, z)p(y)dy. Inserting the infinitesimal mean into the expression for dS(t ) and using the fact that r(z) −r + σ 2 e − c(z,z) has zero mean then yields the expected selection differential as wherem (z,z) = r(z) −c(z,z) ( A 6 ) is the expected fitness of an individual with phenotype z that except for a trivial additive constant not affecting selection, is its Malthusian fitness minus the covariance between its growth rate and the growth rate of the total population. More precisely, the quantity subtracted from the deterministic growth rate r(z) in equation (7) can be expressed as which is the infinitesimal covariance between ln n(z) and ln N. From this, Lande (2007) derived the stochastic generalization of the gradient formula of Lande and Arnold (1983) for the response to selection, showing that the mean Malthusian fitness in that formula simply can be replaced by the longrun growth rater(z) given by equation (2). To see this, we use the formula ∇ p(z;z) = P −1 (z −z)p(z;z), valid for the normal distribution p(z;z). To find ∇σ 2 e (z) we need to evaluate ∇[p(y;z)p(z;z)] in the expression for σ 2 e (z) in equation (1) which is ∇ p(y;z)p(z;z) + p(y;z)∇ p(z;z) giving By the same substitution we find ∇r(z) = P −1 cov[r(z), z] giving ∇r(z) = P −1 cov[m(z,z), z] = P −1 EdS. Furthermore, using E(dz|z) = GP −1 EdS, this yields the generalized gradient formula showing that the long-run growth rate defines an adaptive topography for expected evolution. Using the infinitesimal moments of p(z) and p(y) given by equation (5), Lande (2007) also showed that the infinitesimal covariance matrix for the selection differential vector is where the covariance refers to the distribution p(y)p(z) of (y, z) for two randomly chosen individuals from the population. From this the infinitesimal covariance matrix for the response dz is Finally, using the above relationship Lande (2007) With a k-dimensional phenotypic vector equations (2), (3), (8), (10), and (11) give all the infinitesimal moments for the (k + 1)-dimensional diffusion (z, ln N ).

Natural Selection
Considering equation (A7) in for the expected selection differential and ignoring the stochasticity in dS, the model is within the framework of Fisher's derivation of his fundamental theorem for natural selection in continuous time (Fisher 1930). We can then write Lande's fitness function given by equation (A6) asm(z; x), where x =z is kept constant when considering only natural selection. In general, this is a model with frequency-dependent selection because the fitness function, which is a function of z, has z as parameter. The mean phenotype must then be considered as a component of the environment affecting the whole population and the total change in mean phenotype will be affected by the environmental change inz caused by the natural selection, which is an example of what Fisher called "deterioration of the environment." His fundamental theorem only deals with natural selection and not the effects of "deterioration." Writing h(z, x) = c(z, x)p(z)dz, it appears that σ 2 e (z) = h(z,z). The mean fitness is then m(z, x) = m(z, x)p(z)dz =r(z) − h(z, x).
It follows from Fisher's theorem that the expected change in mean fitness due to natural selection is where V A [m(z,z)] is the additive genetic variance of the fitnessm(z,z). Because the covariance function c(y, z) is symmetric in y and z, we also see that h(z, x) = h(x,z), so that (d/dt )h(z, x)| x=z = 1 2 (d/dt )σ 2 e (z). Inserting this into the middle part of the above equation finally leads to the interpretation of Fisher's fundamental theorem for natural selection in a stochastic environment as that is, the expected rate of change in the long-run growth rate due to natural selection equals the additive genetic variance of the fitness as it was defined by Lande (2007). However, if the covariance function c(y, z) is not constant, the fitness function of z given by equation (A6) always depends onz, which must then be considered as an environmental variable, meaning that nontrivial stochasticity in selection always implies frequency-dependent selection. As a consequence, there is also an additive component dr(z) det in the expected change in the long-run growth rate due to the deterioration of the environment dz nat caused by natural selection. Utilizing the fact that ∇ x h(z, x)| x=z = 1 2 ∇σ 2 e (z), then gives If σ 2 e is constant then this component is zero.

The Diffusion Approach
The derivation of Lande (2007) was based on using the infinitesimal moments of n(y), n(z) and N to derive the moments of p(y) = n(y)/N and p(z) = n(z)/N. Because these calculations are valid for all pairs (y, z), their moments determine uniquely the moments ofz = zp(z)dz given by equations (A3) and (A4).
In that derivation, the moments of n(y), n(z) and N did not depend on N orz. Engen et al. (2013) showed that the moments could more generally also depend on population size N, implying that the parameters occurring in equations (A3) and (A4) may, in general, also be functions of N, so that the formulas derived by Lande (2007) can be used to analyze models with density dependence.
In the present model that also includes frequency dependence, the moments also depend onz. Formally, we then have to consider X 1 = n(y), X 2 = n(z), X 3 = N and X 4 =z as a four-dimensional diffusion, but still use the approach of Lande (2007) of calculating the moments of p(y) = g 1 (X ) = X 1 /X 3 and p(z) = g 2 (X ) = X 2 /X 3 , where X = (X 1 , X 2 , X 3 , X 4 ) T . If we do this for all pairs (y, z), then their infinitesimal moments determine uniquely the moments of the target variable X 4 .
In this model, we then use the general transformation formulas for multivariate diffusions (Karlin and Taylor 1981) to find the moments of g 1 (X ) and g 2 (X ). These moments are linear combinations of derivatives ∂g i /∂X j and ∂ 2 g i /(∂X j ∂X k ), where all coefficients are moments of the diffusion X . Because the functions g i (X ) do not include X 4 , all derivatives with respect to X 4 are zero, so that this component only occurs in the coefficients of the other terms. The result of the transformation is therefore the same as if X 4 was some constant parameter and we only considered the transformation of the three-dimensional diffusion (X 1 , X 2 , X 3 ) with X 4 as a parameter occurring in the moments.
As a consequence, the fourth variablez may be a parameter in the moments of n(y), n(z), and N, and it can simply be plugged in as such in the derivation of the moments of p(y) and p(z) that finally determine the infinitesimal moments ofz.

Frequency-Dependent Resource Requirement
Let R(z) be a function determining how much an individual with phenotype z contributes to within-species competition for resources. Then, under the model described in the main text the growth rate of ln N at a given N is r(z; N,z) =r 0 (z) − 1 2 σ 2 e (z) −γ(z)g(NR(z)).
Assuming weak selection so that fluctuations in population size are much faster than those of the mean phenotype, this growth rate has mean of zero under the almost stationary fluctuations of N so that Eg(NR(z)) =s(z)/γ(z) = Q(z), wheres(z) =r 0 (z) − 1 2 σ 2 e (z). According to equation (4) in the main text, the expected evolution conditioned on N andz is now