Ecological dynamics and large scale phenotypic differentiation in density-dependent populations

Spatial differentiation of phenotypes is assumed to be determined by a combination of fluctuating selection producing adaptations to the local environment and a homogenizing effect of migration. We present a model with density regulation and a density-dependent fitness function affected by spatiotemporal variability in population size driven by spatially correlated fluctuations in the environment causing fluctuating rand K -selection on a set of traits. We derive the variance in local mean phenotypes and show how the spatial scales of the correlations between the components of the mean phenotype depend on ecological parameters. The degree of spatial differentiation of phenotypes is strongly influenced by parameters affecting ecological dynamics. In the case of a one-dimensional character the geographical scale of variation in the mean phenotype has simply an additive term corresponding to the Moran effect in population dynamics as well as a term determined by dispersal and strength of local selection. The degree of phenotypic differentiation increases with decreasing strength of local density dependence and decreasing strength of local selection. These results imply that the form of the spatial autocorrelation function can reveal important information about ecological and evolutionary processes causing phenotypic differentiation in space. © 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
One striking biological phenomenon is the decrease in similarity in many characteristics of populations with increasing geographical distance. For example, a large number of studies covering a wide array of taxa have revealed that genetic differentiation increases with increasing distance between localities (e.g. Rousset, 1997;Templeton, 2006;Aguillon et al., 2017). Similarly, mean phenotypes including both morphological and life history traits of individuals from populations located close to each other in space are more similar than when individuals sampled from distantly located populations are compared (Gould and Johnston, 1972;Zink and Remsen, 1986;Linhart and Grant, 1996;Conover et al., 2006). Finally, another general pattern that has emerged is that the degree of temporal synchrony in population fluctuations decreases with increasing distance (Myers et al., 1997;Ranta et al., 2006;Saether et al., 2007).
Despite these general patterns of isolation by distance in a wide range of population characteristics, the underlying * Correspondence to: Centre for Biodiversity Dynamics, Department of Mathematical Sciences, Realfagsbygget, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway E-mail address: Bernt-Erik.Sather@ntnu.no (B.-E. Saether).
processes affecting the degree of differentiation with distance between the localities are poorly known (Zink and Remsen, 1986;Linhart and Grant, 1996). This occurs even though many classical models in population genetics focus on how geographical structure affects spatial differentiation of genotypes or phenotypes. Haldane (1930), considering a single population receiving immigrants from a large neighboring population with an allele that was deleterious in homozygotes, showed that for a sufficiently high dispersal rate the locally adapted allele would be lost. Wright (1932) also formulated his shifting balance theory in a spatial context, proposing that natural selection will, together with random perturbations in gene frequencies due to random genetic drift, produce adaptive changes (see Svensson and Calsbeek, 2012). Gene flow also constitutes a critical element of this theory because under natural selection alone the population would be trapped at one of the adaptive peaks (Slatkin, 1987;Goodnight, 2012). These early models encapsulated the basic idea that spatial differentiation is the outcome of two opposing processes: natural selection causing adaptations to local environmental conditions and gene flow reducing the mean fitness of the population (Felsenstein, 1976;Lenormand, 2002). Several results are known for small scale genetic differentiation by analysis of models with no selection. An important extension of the earliest theories in population genetics was https://doi.org/10.1016/j.tpb.2019.04.005 0040-5809/© 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). provided by Wright (1943), who proposed that spatial similarity in genetic composition between populations should increase with increasing migration. Wright derived this relationship by a model in which all populations (demes) had constant size and equal probability of exchanging individuals with each other, but considered also continuous populations. Kimura (1953) introduced a more realistic spatial structure than the island model, assuming that populations exchanged migrants only with their immediate neighbors. In this stepping stone model the expected divergence between alleles sampled from different populations also increases with distance between the populations (Slatkin, 1993;Barton et al., 2013). In general, differences between populations increase rather fast at short distances for then to become almost independent of distance. Such a pattern of isolation by distance in genetic divergence has also received considerable empirical support (Rousset, 1997;Templeton, 2006;Aguillon et al., 2017).
An important advance in our understanding of processes affecting phenotypic differentiation in space was provided by Lande's (1991) model of isolation by distance in a quantitative trait. He showed that at equilibrium the geographic pattern of phenotypic variation depended on genetic drift, amount of genetic variance, strength of stabilizing selection and the spatial scale of dispersal, but was independent of the form of the dispersal function. He was also able to derive the covariance between expected phenotypes of individuals as a function of distance for both one-and two-dimensional habitats. However, this approach was based on several simplifying assumptions (e.g. a constant optimal phenotype, the same strength of selection anywhere in space and constant population size and environment). As in the early theories in population genetics (Slatkin, 1985(Slatkin, , 1987, important ecological processes such as the effects of density dependence and stochastic fluctuations in the environment were ignored. In spatial ecology, Moran (1953) derived one of the few quantitative predictions in population ecology. He showed that the temporal correlation in fluctuations of the size of two populations at different locations was equal to the correlation in the common environmental noise in population growth rates if the populations were subject to the same linear form of density regulation. In this model only environmental noise affecting all individual equally is considered which is reasonable for large population (Lande et al., 2003). For smaller populations demographic noise affecting individuals independently (Engen et al., 1998) must be included and will make correlations smaller. In a different approach, Lande et al. (1999) included density dependence and dispersal to derive a very simple formula for the spatial scale of population synchrony, measured as the standard deviation of the spatial autocorrelation function for population density when this is scaled to become a distribution. The squared scale has a minimum occurring under no migration equal to the corresponding squared scale of the environmental noise. With migration a term ml 2 m /γ is added, where m is the dispersal rate, l 2 m is the variance of the dispersal distance, and γ is the local strength of density regulation. The Moran (1953)-effect then becomes a special case (m = 0) of this more general model (see also Engen et al. (2002) and ). These theoretical results have in turn provided a framework for empirical analyses of factors affecting the degree of population synchrony over larger distances Grøtan et al., 2005;Saether et al., 2006Saether et al., , 2007Saether et al., , 2008a.
Here we use the ecological model of Lande et al. (1999) as a basis for analyzing the variance and the spatial scale of the autocorrelation of mean phenotype of a multinormally distributed character in a fluctuating environment. This provides an extension of the approach by Engen and Saether (2016a) and Hadfield (2016), who derived expressions for the degree of spatial differentiation under the influence of dispersal of phenotypes and fluctuations in fitness optima, but assuming constant population densities and no density-dependent selection. In the present model we also include density dependence and fluctuating population densities affecting fitness, generating continued fluctuating r-and K -selection (Lande et al., 2009;Engen et al., 2013), meaning that there is selection for fast growth at small densities and for competitive ability at large densities. In our model the populations are continuously tending to increase their mean fitness or the Malthusian parameter (Saether and Engen, 2015), relative to the given density, according to Fisher's (1930) fundamental theorem of natural selection. However, variation in population density due to a fluctuating environment creates changes in the vital rates and the Malthusian parameter generating phenotypic change in addition to that given by Fisher's theorem (Frank and Slatkin, 1992;Edwards, 1994;Frank, 1997;Kirkpatrick, 2009). Under weak local selection this will create stationary but slow fluctuations in mean phenotypes and hence affect the magnitude of adaptive divergence recorded in space.
The aim of the present model is to analyze how fluctuations in local population size, local density dependence, local selection and dispersal affect the spatial distribution of mean phenotype. We derive scaling results for mean phenotypes valid over large distances, ignoring minor local fluctuations in mean phenotypes due to genetic drift. We argue, based on previous results (Engen and Saether, 2016a,b) summarized in Appendix A, that demographic stochasticity in the population dynamics and local genetic drift can be ignored when population densities exceed certain values determined by the spatial scale of environmental fluctuations.
Although the model is based on several simplifying assumptions, it still includes a number of basic ecological and evolutionary processes and can be considered as a theoretical null model for how different parameters are expected to affect the large scale pattern of spatial differentiation in mean phenotypes in stationary fluctuating populations spread out over a large area. This also provides a theoretical foundation for analyzing how different ecological processes will affect the form of the spatial autocorrelation of phenotypes, which has a long history in analyses of phenotypic and genetic differentiation in a variety of taxa (e.g. Slatkin and Arter, 1991;Sokal and Jacquez, 1991;Sokal and Oden, 1991;Sokal et al., 1998).

Outline of the modeling approach
We shall consider population density N(x, t) at location x = (x 1 , x 2 ) at time t as well as the vector of mean phenotypes z(x, t) as continuous spatio-temporal fields, neglecting demographic stochasticity and random genetic drift. In Appendix A we summarize the main results of two previous papers (Engen and Saether, 2016a,b) on demographic noise in continuous spatial population models, defining two crucial parameters. Writing ρ e (h) for the spatial autocorrelation of the environmental noise in growth rates at locations at separation h, the characteristic area A 0 = ∫ ρ e (h)dh can be considered as an area in which individuals are similarly affected by the environment. The spatial demographic effect s d = σ 2 d /(KA 0 σ 2 e ) is an expression of the importance of demographic noise relative to environmental noise.
Here σ 2 d and σ 2 e are the demographic and environmental variances in population growth, respectively, and K is the carrying capacity at the density scale. The main conclusion is that if s d is small compared to 1, then demographic noise and the random genetic drift it generates can be neglected. Typical values of σ 2 d and σ 2 e are of order 1 and 0.01 (Lande et al., 2003;Saether et al., 2013), respectively, giving an s d of order 100/(KA 0 ), where KA 0 is the mean number of individuals in A 0 , so that large density K or large spatial scale of the environmental noise is required to give a small s d . In the following we assume that the carrying capacity and the environmental spatial scale are large enough for using a continuous spatial population model with environmental noise only.
Even if drift operating over relatively small distances is ignored, the problem is extremely difficult to analyze in general. There are two spatial fields in two dimensions, one for population densities and one for multivariate mean phenotypes that interact by temporally correlated environmental noise and dispersal, and undergo stochastic fluctuations in time. However, we avoid using time in the formulas writing N(x) andz(x) for the spatial fields to simplify the notation when this is not required. So, if there are q different phenotypes to be studied jointly, there is altogether a (q + 1)-dimensional field in two spatial dimensions. In order to find any transparent results that can give some general insight into this complex process it is crucial that some assumptions of weak selection are made. If selection is sufficiently weak, the temporal and spatial variation in mean phenotypes will be slow and often not observable in ecological time, but faster and observable at an evolutionary time scale over thousands of years. The population densities, however, fluctuate very fast in evolutionary time, commonly with a return time to equilibrium from some months up to some few decades, although there is less variation among species if time is measured in generations. The coupling between the fields can then be dealt with by assuming that the genetics is fixed in many mathematical arguments concerning fluctuations in ecological time. In this way, the fluctuations in the density field can be considered as stochastic noise in the evolutionary process, and can even be approximated by a noise with no temporal but a spatial correlation at the evolutionary time scale, determined by the spatial population dynamics governed by local density regulation, dispersal and spatially correlated environmental noise.

Ecological dynamics
The logistic model for population density N in continuous time with environmental noise and no spatial or genetic components takes the form dN = N(r − ηN)dt + σ e NdB(t), where σ 2 e is the environmental variance and informally dB(t) are increments of standard Brownian motions during time dt with zero mean and variance dt. Writing γ = K η, where K = r/η is the mean population size in the corresponding deterministic model (σ 2 e = 0), yields dN = N(r − γ N/K )dt + σ e NdB(t). In this initial formulation with no genetics it appears that γ = r, but as we later want to study genetic variation in r and γ (introducing r(z) and γ (z), where z is the phenotype vector) it is required that we keep the intrinsic growth rate r and the strength of density regulation γ as distinct parameters. The characteristic return time to equilibrium is then 1/γ (May, 1974). Expressing the model in terms of n = N/K −1 then yields dn = −n(n+1)γ dt +σ e (n+1)dB e . Linearizing at n = 0 (N = K ) and approximating the stochastic term by its average σ e dB e (approximating N/K by 1), yields an accurate approximation for small and moderate population fluctuations (Engen, 2017) as dn = −γ ndt + σ e dB(t).
Following Lande et al. (1999) we write N(x, t) for the population density at location x = (x 1 , x 2 ) at time t, assuming that mean densities are large enough for demographic noise to be ignored. Under small fluctuations in local population density and density independent dispersal at rate m the linearized dynamics of the proportional deviation from carrying capacity K , n(x, t) = N(x, t)/K − 1, are given by where f (h) is the distribution of dispersal distance and the last term is the noise. The factor dB e (x, t)  Based on this model Lande et al. (1999) derived a remarkably simple formula for the spatial scale of population synchrony. Scaling autocorrelation or autocovariance functions to become distributions, they defined the spatial scale along a given direction in space as the standard deviation of these distributions in that direction. Writing this spatial scale for the autocorrelation of population density as l n (n and N have the same spatial scale), for the noise as l e , and the standard deviation of dispersal distance as l m , this yields This scaling result based on linearization of the dynamics is formally only valid in the limit that σ 2 e approaches zero, but Engen (2017) analyzed its accuracy and showed that it is very accurate under moderate and even rather large population fluctuations. The Fourier transform of the spatio-temporal covariance Here, we analyze a rather similar but a bit more complex spatial model for mean phenotypes. In this model the process n(x, t) appears as noise and we show in Appendix C that it can be ap- The Fourier transform of this autocovariance function, which is what is required to perform numerical calculations, is given by Eq. (C.2) in Appendix C. We also show in this appendix that the spatial scale of this noise term is l 2 that this is larger than the spatial scale of n(x, t) is due to the temporal autocorrelations in n(x, t).

Quantitative genetic model
Let z = (z 1 , z 2 , . . . , z q ) T , where T denotes matrix transposition, be a phenotypic vector with a multinormal distribution p(z; x) among individuals at location x with meanz(x) varying in space and time and phenotypic covariance matrix P and additive genetic matrix G, which are both assumed to be constant in space and time. The time variable has been omitted in the mean phenotype to simplify the notation. The local change in mean phenotype will have two additive where the first is due to dispersal and the second is the response to local selection. It is shown in Appendix D that with an additional stochastic term (when conditioned on the field of mean phenotypes) that is small compared to the stochastic term in the response to be defined by Eq. (4) below. Hence, in the following we neglect the stochastic term in the migration component.
We now consider the quantitative genetic model where r(z) and η(z) vary among individuals at a given location x according to the multivariate normal density p(z; x). For simplicity we assume that each individual is equally affected by the environment so that the environmental variance σ 2 e is constant. The mean deterministic growth rate for small densities is thenr(z) = ∫ r(z)p(z; x)dz and similarly forη(z), while the corresponding long-run growth rate, the expected growth rate on the log scale, is (Lande, 2007). These are functions of the mean phenotypez at location x since p(z; x) has meanz not shown in the notation.
Introducing γ (z) = K η(z) for individuals with phenotype z as we did in the non-evolutionary Eq.
(1), we show in Appendix E that the response to selection takes the form where n = n(x, t) is the spatio-temporal process given by Eq.
Actually, to obtain a realistic model this function must have a maximum in the genetic model for fluctuating r-and K -selection because expected population size would otherwise evolve towards infinity (see Appendix E and Engen et al., 2013 for details).
There is no simple transparent formula for the autocovariance function for the noise, but in Appendix C an approximation for its Fourier transform is given by Eq. (C.2). This enables numerical calculation of the autocovariance function by the backward transformation (see Appendix B). From Eqs. (3) and (4) the total change in mean phenotype dz = d mz + d Rz has three terms. The first term in Eq. (4) is deterministic, representing expected response to selection towards z * by the adaptive topography Q (z). Next, there is second deterministic term expresses the effect of migration by Eq. (3), while the third appearing in Eq. (4) is a stochastic noise term with zero mean generated by fluctuations in population density that are fast compared to those ofz. Assuming weak selection we can approximate Q to the second order at its maximum writing where the H −1 ij = −∂ 2 Q (z)/(∂z i ∂z j ) evaluated atz = z * are the elements of the Hessian matrix H −1 for the function −Q . Using this approximation we find ∇Q (z) = −H −1 (z − z * ). Inserted into the equation for the response to selection which is averaged over the stochastic fluctuations in n (and are fast compared to those of z and considered as noise in the genetic model) and conditioned on the fieldz(x) can be written as where Γ =γ GH −1 is a matrix expressing the strength of local stabilizing selection. The total change in mean phenotype is now the sum of the dispersal and response component, that is, The dynamic equation (6) has striking similarities with Eq. (1) (with n replaced byz) althoughz is in general a vector while n is a scalar. First there is a linear term sending the mean phenotype towards z * , then a linear dispersal term of the same form as in (1), and finally a noise term which as in Eq. (1) can be approximated by white noise with spatial autocorrelation that has spatial scale l 2 e + 2ml 2 m /γ (see Appendix C).

The scale of spatial autocorrelations of mean phenotypes
In Appendix F we give the general expression for the spatial scale l ij of all covariances/variances betweenz i andz j , that is, the spatial scale of the autocovariance function cov[z i (x, t),z j (x + h, t)], in any direction. The expression is complex and not transparent, so here we focus on the analysis of a single trait z with additive genetic variance G. Provided that trade-offs between rapid growth expressed by large r, and ability to compete under large densities, expressed by a small γ , has evolved so that individuals with large r tend also to have large γ , Q (z) will have a maximum at z * with a second derivative −H −1 at z * .
Eq. (F.4) in Appendix F for a one-dimensional trait then yields only one eigenvalue λ 1 = G/H. Omitting the subscript 1 the squared spatial scale of cov[z(x),z(x + h)] takes the form where Γ =γ GH −1 is the strength of stabilizing selection in the sense that the local response is d Rz two terms at the right side of Eq. (6) represent the direct effect of fluctuations in population size N, which has a spatial scale larger than that of N itself (exactly by a factor 2) given by Eq.
(2) due to the temporal autocovariances in N as outlined in Appendix C.
In other words, these two terms correspond to the Moran-effect (Moran, 1953) in population dynamics as it is the squared spatial scale of the noise in the evolutionary process. The last term is the component generated by dispersal through local selection, and is the analogy of the term ml 2 m /γ in the expression for the scale of N(x, t) given by Eq.
(2). This shows that dispersal may have a large effect on the spatial scale of mean phenotype if the local selection is weak.

The magnitude of fluctuations in mean phenotypes
In Appendix G we show how to compute the spatial autocovariance function by Fourier analysis, first deriving the Fourier transform for autocovariance functions such as where F e (u) and F m (u) are the Fourier transforms of ρ e (v) and the distribution f (v) of dispersal distance, respectively, in the isotropic case (see Appendix B for definitions).
If there is no dispersal, m = 0, so that the correlated noise is the only effect that combine locations, using ∇γ /γ = ∇ lnγ and the backward transformation (2π Because G is a factor in Γ this variance is proportional to G as well as to σ 2 e . The factor (∇ lnγ ) 2 represents density dependent selection by expressing how stronglyγ (z) depends onz atz = z * . However, both density independent and density dependent selection has additional effects by determining the factor −Q ′′ (z * ) of Γ . Notice that Eq. (8) for the variance in mean phenotype refers to the simple model with no dispersal. By inspection of Eq. (7) we see that increased dispersal makes the integral in the equation smaller and hence induce reduction in the variance corresponding to a smoothing effect of dispersal that tends to reduce the magnitude of spatial variation (Hadfield, 2016).

A univariate illustrating example
We now illustrate the application of this theory by a simplified univariate linear model. First notice that if γ (z) represents a local density regulation in the sense that increased local density yields a decrease in population growth rate. Then γ (z) > 0 so that ln γ (z) may take any real value just like the intrinsic growth rate r(z). Accordingly, as an illustration of the univariate model we consider the linearizations of r(z) = r 0 + αz giving s(z) = s 0 + αz, and ln γ (z) = ln γ 0 + βz at the equilibrium z * as a first order approximation to more general models valid for small fluctuations ofz in time and space. The results are independent of K if we scale the density by choosing γ 0 so thatγ (z * ) =s(z * ) giving Q (z * ) = 1. We also choose to scale z by a factor so that the additive genetic variance G = 1. Then, from (8) In this model, the standard deviation SD(z) = √ var(z), relative to its maximum at m = 0 given by Eq. (9), decreases as the dispersal rate m increases for different values of the local strength of selection Γ (Fig. 1). This illustrates that the smoothing effects on the mean phenotype of dispersal are dependent upon the strength of selection. Similarly, the spatial autocorrelation for the mean phenotype is also dependent on the strength of local selection (Fig. 2). Strong local selection decreases the scaling of the spatial autocorrelation function. However, the actual value of the autocorrelation, shown in Fig. 2, decreases with the scale shown as vertical lines. This occurs because the form of the autocorrelations deviates more from Gaussian curves towards functions with heavier tails under very weak strength of selection due to the larger effect of dispersal.

Discussion
The major conclusion of our model is that processes affecting the magnitude of spatial synchrony in population dynamics such as environmental covariation in space, migration and strength of local density dependence (Lande et al., 1999) will influence the variance and spatial scale of covariances for mean phenotypes. This occurs even though the model is homogeneous in space and therefore does not account for adaptations to permanent local environments varying in space (e.g. Hereford, 2009). Hence, the results should be viewed as a theoretical null model in the same way as for example the spatial model of Lande et al. (1999). However, as such, it gives considerable insight into how different parameters affecting the population dynamics influence the variance of the mean phenotype in space as well as the spatial scale of its autocorrelation function. Thus, this provides a theoretical framework for interpretation of differences in patterns of variation in spatial autocorrelation functions that are often used in analyses of phenotypic differentiation (Slatkin and Arter, 1991;Sokal et al., 1996).
The present model is an analysis of variation and large scale synchrony in covariances of mean phenotypes. Based on previous results we have argued that demographic stochasticity and random genetic drift can be ignored in the calculations provided  that the mean number of individuals in characteristic areas defined by the spatial scale of environmental fluctuations, is large. For example, the spatial scale of environmental fluctuations in population dynamics of butterflies in the tropics is of order 2 km (Lande et al., 2003), giving a characteristic area of A 0 = 2π l 2 e ≈ 25 square kilometers. Although there is a number of rare species, most species will have a very large number of individuals within A 0 so that the demographic coefficient s d is very small compared to 1. Similarly, for species with very large spatial scales of covariation in the population dynamics, such as marine fishes (Myers et al., 1997), a typical spatial scale may be as large as 500 km, corresponding to A 0 close to a million square kilometers that also, for most species, contains a huge number of individuals making s d small. Birds (Saether et al., 2007) and mammals (Grøtan et al., 2005(Grøtan et al., , 2008 represent species with typical intermediate spatial scale of covariation in the population dynamics, at the order of 10 km, giving A 0 at about 600 square kilometers. Many species will then have a number of individuals within such an area that will make s d sufficiently small to ignore demographic noise, although some rare species will not. In accordance, based on extensive comparative analyses of geographical variation in quantitative genetic (Q ST and neutral F ST ) markers Leinonen et al. (2008) concluded that spatial divergence due to natural selection and local adaptation was mainly the norm. Our evolutionary model is also based on several other simplifying assumptions. First, we assume that the phenotype is normally distributed among individuals with a constant covariance matrix (Lande, 1976) as in the model for r-and K -selection of Engen et al. (2013). There are many factors, including dispersal, that will affect this distribution and its covariance matrix, but we assume that the deviations from normality and fluctuations in covariances are small. Second, the local dynamics of population size N is assumed to have small or moderate fluctuations so that the dynamics can be approximated by a linear model in n = N/K − 1. Third, we assume weak local selection so that changes in mean phenotype are much slower than changes in population size. Comparative analyses indicate that directional selection coefficient on phenotypic traits in natural populations often are significant but quite small (Kingsolver et al., 2001(Kingsolver et al., , 2012Hendry, 2017) and consistent over time (Morrissey and Hadfield, 2012). Fourth, we assume that the environmental variance σ 2 e is constant, independent of the mean phenotype. With these approximations we argue that the stochastic term in the increment ofz(x) due to dispersal can be ignored. Fifth, we approximate temporal fluctuations in local population size, which has a temporal autocorrelation depending on the local density regulation as well as the dispersal, by a noise process with no temporal correlations considered as noise in the evolutionary process. The spatial autocorrelation of this noise has squared spatial scale l 2 e + 2ml 2 m /γ , while the spatial scale of the autocorrelations in population size N is l 2 e + ml 2 m /γ (Lande et al., 1999). Sixth, we assume no phenotypic plasticity in the environmental response, which can influence the spatial scaling of the mean phenotype (Hadfield, 2016). In spite of these simplifying assumptions, our model provides an extension towards increased ecological realism of previous models that have analyzed how spatial differentiation of phenotypes depends on distance (Lande, 1991;Engen and Saether, 2016a;Hadfield, 2016) because it includes density dependent selection as well as stochastic fluctuations in the environment generating stationary fluctuations in mean phenotype. For example, in the model of Engen and Saether (2016a) population densities are constant (like in Lande, 1991) but the optima of the fitness function θ(x) fluctuate both in time and space. In that model we obtained some simple scaling result for spatial variation in the mean phenotype: fluctuations in θ have spatial scale l θ and the temporal process for θ has an average return time to equilibrium of 1/γ θ so that γ θ represents the strength of the force driving θ towards its spatial and temporal equilibrium value. Writing Γ for the strength of local selection in that model, we found that the squared spatial scale of the mean phenotype was l 2 θ + ml 2 m /(Γ + γ θ ). Here it appears that a short return time to equilibrium of θ reduces the scaling effect of local selection through Γ , since selection then continuously is disturbed by the fast fluctuations in θ. Still, there are some links between the results of the two models: dispersal increases the scale while increased strength of local selection decreases it. Including all ecological effects in a single model would be interesting, but challenging, and is unlikely to provide general transparent results.
Although we have given the general solution for a multivariate character, we have chosen to focus the discussion on the results for a single character only because the multivariate results are rather complex with little transparency. However, even if the parameters r(z) and γ (z) determining the Malthusian fitness (Saether and Engen, 2015) in this model are functions of a multivariate phenotype, their variation may often be determined by some function of phenotypes, for example a linear combination (corresponding to a first order Taylor expansion), so that the univariate model with this single phenotype used as the only single trait affecting fitness may often be realistic. Hence, the more transparent univariate results illustrate well the links between the spatio-temporal variation in population dynamics and the evolutionary processes.
The univariate analysis yields two major observations, summarized in the expression for the autocovariance function c z (v) in the isotropic model given in Appendix F. Plugging in distance v = 0 this yields an expression for the spatial variance inz as expressed by Eq. (8) Even under no dispersal we can draw some interesting conclusions. First, in order to obtain a variable mean phenotype there must be some environmental variance σ 2 e generating fluctuations in local population size. Second, in order to have any local response to selection the trait must have some additive genetic variance G so that it actually has some heritability. Third, there must be some gradient in the local strength of density regulationγ (z) (and its log) at the phenotype z * maximizing Q (z) =s(z)/γ (z). This means that the effect of density on vital rates of individuals must depend on the phenotype (see Saether et al., 2016a for an empirical example), which is a requirement for continued fluctuating balance between r-and K -selection. Finally, the (co)variance is inversely proportional to the local strength of stabilizing selection Γ so that large spatial fluctuations in mean phenotype may indicate weak local selection. This effect of Γ is intuitive from simple results in first order autoregressive models in discrete time or the Ornstein-Uhlenbeck model in continuous time. If such a process ξ has infinitesimal mean and variance a − bξ and σ 2 , the stationary variance is σ 2 /(2b), where b corresponds to our strength of selection Γ , always making the mean phenotype approach its temporal mean (see also Chevin et al., 2015Chevin et al., , 2017. Our model yields a particularly simple result for the spatial scale of the autocorrelation of mean phenotype given by Eq. (7). As in the ecological model of Lande et al. (1999), if weak density regulation is weak, dispersal can generate a very large spatial scale for the autocorrelations of population size Eq. (2). Our formula provides a similar result with respect to local strength of stabilizing selection. If the strength Γ is small, dispersal will have a large effect on the spatial scale of the autocorrelation in mean phenotype, in addition to the general effects of dispersal as homogenizing agent for phenotypic variation in space as expressed by Eq. (8) and illustrated in Fig. 2 (see references in Engen and Saether, 2016a;Hadfield, 2016).
The univariate model with linearization of s(z) = s 0 + αz and ln γ (z) = ln γ 0 + βz represents a simple decomposition of the local selection into a density independent component expressed by α and density dependent selection expressed by β. When z is scaled so that the mean is at zero and the additive genetic variance is 1, the local strength of selection Γ is simply the product αβ. The intrinsic growth rate s * at the equilibrium is α/β and the variance of the mean phenotype in space and time under no dispersal is 1 2 σ 2 e β/α = 1 2 σ 2 e /s * . Thus, the variance will increase with increasing density dependent selection, which is actually the link between the fluctuations inz and those in population size N. When α gets very small so that there is only density dependent selection, there will in the limit (α = 0) be no genetic equilibrium. Comparative studies have revealed large spatial variation in directional selection in natural populations (Siepielski et al., 2013), but the mechanisms involved in generating the geographical patterns in selection are still poorly known (Kingsolver et al., 2012).
To illustrate further the relationship between phenotypic differentiation and population dynamics, it may be illustrative to consider some typical parameter values. Environmental stochasticity is often referred to based on analyses of vertebrate population dynamics Lande et al., 2003;Saether et al., 2013Saether et al., , 2016b as being moderate when σ 2 e = 0.01 reflecting a standard deviation of 0.1 in the intrinsic growth rate, while σ 2 e = 0.04 with standard deviation 0.2 represents rather large fluctuations. The value of s * , the intrinsic growth rate of the species, that is, the growth rate in the absence of density regulation, varies considerably even among closely related species Saether et al., 2008bSaether et al., , 2013Saether et al., , 2016b. Similarly, the deterministic growth rate is often correlated with the generation time Bjørkvoll et al., 2012), expressing the slow-fast continuum of life history variation (Saether and Bakke, 2000). For example, s * = 0.5 may represent a small fast organism while s * = 0.02 is in the typical range of large species with long generation times. The largest variance in mean phenotypes under no dispersal is expected for such species in highly fluctuating environments, exemplified by 1 2 σ 2 e /s * = 1 2 · 0.04/0.02 = 1 so that the variance equals the additive genetic variance of the character under no dispersal. In general, dispersal decreases the variance (Fig. 1). A fast growing species in a moderately fluctuating environment will at the other end have variance 1 2 · 0.01/0.5 = 0.01 so that the standard deviation of the mean phenotype is only 10% of the standard deviation of the breeding values under no dispersal.
To summarize, our model indicates that there should be a relationship between the pattern of spatio-temporal variation in population dynamics and degree of spatial differentiation in phenotypes (Kirkpatrick and Barton, 1997). These patterns should be analyzed using spatial variation of individual phenotypes z rather than mean phenotypesz. Because the phenotype z of an individual equals the mean phenotypez plus a deviation from this, which is independent among individuals, the spatial autocovariance inz must be the same as the spatial autocovariance in the individual phenotypes z and can therefore be estimated by standard methods in spatial statistic (Sokal et al., 1998), e.g. using variograms based on individual observations (Cressie, 1993). drift We describe population density N(x, t) at location x = (x 1 , x 2 ) at time t as well as mean phenotypesz(x, t) as continuous spatiotemporal fields, neglecting demographic stochasticity affecting changes in population size in small areas with finite number of interacting individuals, as well as random genetic drift in such subpopulations. Therefore, our theory can only be realistic for description of fluctuations in density and mean phenotypes over areas large enough for ignoring such independent stochastic variation among individuals, and our scaling results are valid only for spatial covariances between mean phenotypes at large distances. Hence, we first need to find the characteristics of the spatiotemporal distribution of population density of different organisms that determine minimum area and distance for application of these results.
Demographic noise in continuous spatio-temporal processes is difficult to model accurately, requiring that we deal with each individual separately and their position in space. We have previously developed a theory for analyzing population dynamics including demographic and environmental noise as well as dispersal . We have also examined selection and genetic drift under fluctuations in optimal phenotype in a spatio-temporal model ignoring density-dependence . In both cases there is density independent dispersal of individuals at a given rate m (a fraction mdt disperse during time dt) and a specified distribution of dispersal distance. The present model has considerably higher degree of complexity because it uses two fields, N(x, t) andz(x, t), interacting through a common dispersal of individuals, while in our genetic model , the field of optimal phenotypes is not affected by dispersal. This makes it impossible to derive any transparent analytical results in the present model in the presence of demographic noise, but still the above papers contain some results that can be used to address how demographic stochasticity may affect the spatial scaling of phenotypic variation . In simple population models with no spatial structure the noise in the dynamics of population size N, defined as the variance of the yearly relative change in population size ∆N/N, is σ 2 e +σ 2 d /N. Here σ 2 e is the environmental variance affecting all individuals in a similar way, whereas σ 2 d is the demographic noise representing independent variation in vital rates among individuals (Lande et al., 2003). In spatial models the environmental effects at two locations at spatial distance h will have a spatial correlation ρ e (h) that usually decreases with increasing distance.  showed that the importance of demographic noise in spatio-temporal population dynamics can be expressed by a parameter s d = σ 2 d /(KA 0 σ 2 e ) called the spatial demographic effect determined by mean population density K and the characteristic area A 0 = ∫ ρ e (h)dh, provided that the environmental correlation vanishes at infinity. If the environmental autocorrelation is isotropic with a Gaussian form with variance l 2 e in any direction, Roughly we can say that individuals within a circle with area A 0 are typically affected in approximately the same way by the environment. The demographic noise can be ignored in computations of the spatial scale of population density if s d is small compared to 1 . The effect of genetic drift in the genetic model is more complex (Eq. (14) in Engen and Saether, 2016b), dependent upon several parameters. A characteristic area in this model for a single trait is defined by the spatial autocorrelation ρ θ (x) for the optimal phenotype θ, as A 0,θ = ∫ ρ θ (x)dx. The effect of drift depends on an effective characteristic population size defined as N e = N 0,θ /(σ 2 d +m), where N 0,θ is the (constant) number of individuals in A 0,θ . From Eq. (4) in Engen and Saether (2016b) we find that, under weak local selection, drift can be ignored in the spatial scale of mean phenotypes if Here τ is the width of the local fitness function with Gaussian form, G is the additive genetic variance of the trait, and V θ is the spatial variance of θ. As a numerical illustration let τ 2 /G = 10 (weak selection), τ 2 /V θ = 10 (so that V θ = G), m = 1 (each individual disperses on average once during a generation) and σ 2 d = 1, gives that N 0,θ must be large compared to 100.
From these two models it appears that demographic noise and genetic drift can be ignored in the evaluation of spatial scale of phenotypic variation when the number of individuals in the characteristic areas are sufficiently large. Although the present model is more complex, the results for these two more simple, but closely related models, give a strong indication that demographic noise and genetic drift have only a minor influence on the spatial scale of phenotypic correlations, provided that the number of individuals affected similarly by the environment is sufficiently large.

Appendix B. The temporally white noise approximation for fluctuations in population size
The mathematical problem is most easily solved by using the Fourier transform of functions f (x 1 , x 2 ) that may be an autocorrelation function, an autocovariance function or a distribution, where (x 1 , x 2 ) = x are two-dimensional spatial coordinates. The Fourier transform of f is then where ω = (ω 1 , ω 2 ) and the integration is taken over the whole plane R 2 . The backward transformation is given by For simplicity of illustrations we shall consider isotropic models where all functions f can be expressed as functions of the For such models the Fourier transforms are functions of u = √ ω 2 1 + ω 2 2 . We then use the notation f (v) for the functions to be transformed and F (u) for their Fourier transforms. The inverse transformation is then given by the univariate real integration is the Bessel function of the first kind of order zero. Lande et al. (1999)  Here F e (ω) and F m (ω) are the Fourier transforms of the spatial autocorrelation of the ρ e (h) noise and the dispersal separation, respectively.
Appendix C. The white noise approximation to the spatiotemporal process n(x, t) The last term in Eq. (6) determined by ndt = [g(N) − Q (z)]dt is a noise term in the evolutionary processz. In the spatial evolutionary process n(x, t) represents a temporally and spatially autocorrelated noise. By the assumption of weak selection the mean phenotype will be a process with much slower fluctuations than the n(x, t) so that a long interval ∆T relative to fluctuations in population size is a short interval in the evolutionary process.
The effects of n(x, t) during successive intervals ∆T 1 , ∆T 2 , . . . can then be considered as independent so that the noise in the evolutionary process can be approximated by a temporally white, but spatially autocorrelated noise. To find this approximation consider an interval ∆T during which the effects of the noise at locations x and x + h are a(x) and a(x where ρ w (0) = 1 and subscript w indicates the use of the white noise approximation. Hence, under weak selection we can approximate cov[a(x), a(x + h)] by σ 2 w ρ w (h)∆T so that the noise n(x, t) can be replaced by a noise σ w dB w (x, t) with no temporal autocorrelation, but a spatial correlation given by E Using the approach of Engen (2001) the spatial scale of this autocovariance function along the first axis is most easily computed as

Appendix D. The dispersal component
Now consider the smoothing effect of dispersal on the field z(x) by first considering the increment of the sum of phenotypes T (x) = N(x)z(x) over all individuals at location x. We assume that the probability that an individual disperse during time dt is mdt, not depending on its phenotype. Dispersal out of location x during time dt yields a reduction mdtT (x) of the sum of phenotypes at x, while dispersal in from the surrounding locations yields an addition mdt Here the time variable has been omitted to simplify the notation. Similarly, N(x) is reduced by mdtN(x) and increased by mdt Now,z(x) = T (x)/N(x) and the increment due to migration during dt is .
Inserting the above expressions for dT (x) and dN(x) then yields splits the increment in mean phenotype conditioned on the field of mean phenotypes into a deterministic part and a stochastic component This stochastic component will generally be very small compared to the deterministic component given by (D.1) for several reasons. First we have based the analysis of a population density field using the linearization given by Eq. (1) in the main text, which is valid only under small or moderate fluctuations in n(x), say with standard deviations smaller than about 0.3. From Lande et al. (1999) and the results in Appendix C it appears that n(x), when studying its effect over long intervals in ecological time, can be approximated by a field with no temporal autocorrelations but a spatial correlation with scale √ l 2 e + 2ml 2 m /γ . If either the density regulation is weak (small γ ) or spatial scale of environmental noise is large (large l e ), then this is much larger than the scale l m of the dispersal. Accordingly, over the area contributing to dz(x) by dispersal the spatial correlations for the scaled densities are approximately 1 so that the density is approximately constant

Appendix E. Quantitative genetic theory
Now consider the local selection at some given location x and time t. In order to simplify the notation the symbols x and t are omitted when not required. The selection will vary in space, but only through the value ofz at each spatial location x. Engen et al. (2013) analyzed a quantitative genetic model for r-and K -selection in a fluctuating environment using continuous time and ignoring spatial effects. Under logistic dynamics the deterministic growth rate of a hypothetical sub-population N sub (z) with individuals of phenotype z in a population in which density regulation is governed by the total population size N, was assumed to have the form [r(z) − η(z)N]N sub (z). Environmental stochasticity is in general described by the density-independent infinitesimal covariance between the changes in sizes of subpopulations of type z and w due to selection only as c(z, w)N sub (z) N sub (w) (see Lande, 2007 for details). The infinitesimal mean and variance for the total population size are then [r(z)−η(z)N]N and σ 2 e (z)N 2 , wherē with a similar definition forη(z), are the population means. Hence, the dynamics of N is determined by the mean values of r(z) and η(z) as well as c(z, w) across the individuals in the population, more precisely after transforming to log scale wheres(z) =r(z) − σ 2 e (z)/2. Assuming perfect correlations between the environmental stochastic effects on phenotypes by assuming a constant c(z, w) = σ 2 e , Engen et al. (2013) showed that the response to selection by this model for a given population size N is wheres(z) =r(z) − σ 2 e /2 is the long-run growth rate of the population in the absence of density regulation. Since we assume weak selection the parameters determining the dynamics of N in space will show very little variation, although N itself may fluctuate considerably. Remembering that the stochastic field of population densities only acts as noise in the genetic model, it can be approximated by its average properties, which are those at the equilibrium phenotype. This is exactly the spatio-temporal model for fluctuations in density analyzed by Lande et al. (1999) with linearized dynamics as expressed by Eq. (1) in the main text.
Writing n = N/K − 1 and defining γ (z) = K η(z), where K is the equilibrium density as in Eq. (1) that can be considered as a scaling factor for densities with no genetic variation, Eq. (E.1) takes the form d R (z|z, N) = G[∇s(z) − (n + 1)∇γ (z)]dt.

(E.2)
Defining the function Q (z) =s(z)/γ (z) the expected density is EN =s(z)/η(z) = KQ (z) which is approximately K under small environmental fluctuations so that n fluctuates around zero.

(E.3)
Hence, the expected evolution, obtained by ignoring the last stochastic term, will always tend to increase Q (z) approaching an equilibrium value z * maximizing Q (z), while the noise term proportional to n = N/K − 1 generates fluctuations in mean phenotype around z * (see Engen et al., 2013 for details).

Appendix F. Solution for the spatial scale of fluctuations in mean phenotypes
In order to solve the full spatial model for a multivariate character we assume that GH −1 can be diagonalized writing GH where the ∇γ (z) are approximated by their values at z * denoted ∇γ .
Notice that Eq. (1) for the scaled population size and Eq. (F.2) for the centered mean phenotype has the same form, although Eqs. (F.2) and (1) yield a vector and a scalar, respectively. To find the spatial covariance functions we utilize the assumption that parameters describing the system has values that ensure stationarity in population size as well as mean phenotypes. For example, there must be trade-offs between growth rates and density regulations that prevent fitness to increase indefinitely. Inserting dū(x) from Eq. (F.2) then yields, after dividing throughout by dt, where γ isγ (z) evaluated at z * . Utilizing that Λ is a diagonal matrix this yields an equation for the (ij)-element C uuij (h) of C uu (h) as where b ij are the elements of the matrix b = M G∇γ ∇γ T GM T .
Let U be the Fourier transform of C uu (h) in the sense that each element in the matrix is transformed separately so that the (i, j)-element of U is U ij (ω) = 1 2π ∫ ∫ e i(h 1 ω 1 +h 2 ω 2 ) C uuij (h)dh.
Similarly we write F and F w for the transforms of the distribution of the dispersal distance f (h) and the spatial correlation ρ w (h) of the noise generated by fluctuations in population sizes.
Utilizing that the integral in Eq. (F.2) is a convolution we then find [2m + γ (λ i + λ j )]U ij = 2mU ij F + b ij F w with the solution This is the Fourier transform of the autocovariance function C uuij (h) = cov[ū i (x),ū j (x + h)]. Using thatū = M (z − z * ) the autocovariance function forz is accordingly M −1 C uu (h)(M −1 ) T , and the matrix of its Fourier transforms with elements Z ij is Z = M −1 U (M −1 ) T . From this we may compute the spatial autocovariance of the components ofz in isotropic models using the backward transformation given by the univariate integration shown in Appendix B.
The squared spatial scale of C zz (h) = cov[z i (x),z j (x + h)] with Fourier transform Z ij (ω) in the direction of the first coordinate, by the way we have defined this scale, is